Bcor loss perturbs myeloid differentiation and promotes leukaemogenesis.

The BCL6 Corepressor (BCOR) is a component of a variant Polycomb repressive complex 1 (PRC1) that is essential for normal development. Recurrent mutations in the BCOR gene have been identified in acute myeloid leukaemia and myelodysplastic syndrome among other cancers; however, its function remains poorly understood. Here we examine the role of BCOR in haematopoiesis in vivo using a conditional mouse model that mimics the mutations observed in haematological malignancies. Inactivation of Bcor in haematopoietic stem cells (HSCs) results in expansion of myeloid progenitors and co-operates with oncogenic KrasG12D in the initiation of an aggressive and fully transplantable acute leukaemia. Gene expression analysis and chromatin immunoprecipitation sequencing reveals differential regulation of a subset of PRC1-target genes including HSC-associated transcription factors such as Hoxa7/9. This study provides mechanistic understanding of how BCOR regulates cell fate decisions and how loss of function contributes to the development of leukaemia.

Introduction

The BCL6 Co-repressor (BCOR) gene is located on chromosome X and encodes a ~180 kDa nuclear protein that is ubiquitously expressed across adult human tissues. BCOR is essential for normal embryonic development; germline hemizygous BCOR loss-of-function mutations are responsible for male lethality while heterozygous mutations in females result in oculocardiofaciodental (OFCD) syndrome, a rare genetic condition characterised by craniofacial, ocular and cardiac abnormalities[1]. Next-generation sequencing studies have demonstrated that somatic BCOR mutations occur in a range of conditions that affect the myeloid and erythroid haematopoietic lineages including acute myeloid leukaemia (AML), myelodysplastic syndrome (MDS), chronic myelomonocytic leukaemia (CMML) and aplastic anaemia[2-6]. The mutations in these diseases almost always result in a premature stop codon and nonsense-mediated decay or protein truncation, strongly suggestive of a tumour suppressor role for BCOR in these contexts. Interestingly however, BCOR can mediate pro-oncogenic functions in some cell types, while in other contexts it behaves as a tumour suppressor[7-11]. Hence, the role of BCOR is highly tissue-specific and should be functionally analysed in a physiologically relevant setting that is appropriate for the disease being studied.

Polycomb group (PcG) proteins are evolutionary conserved chromatin modifiers that regulate a broad array of genes in mammals and play major roles in development and cancer[12]. PcG proteins are present in multi-protein complexes that can be classified into two types, Polycomb repressive complex 1 and −2 (PRC1 and −2). PRC1 and PRC2 possess distinct enzymatic activities: PRC1 ubiquitinates histone 2A at lysine 119 (H2AK119ub) whereas PRC2 di- and trimethylates lysine 27 of histone 3 (H3K27me2/3). The canonical pathway of PRC transcriptional repression, established predominantly from studies in Drosophila and mammalian embryonic stem cells, involves the sequential recruitment and activity of PRC2 followed by PRC1[12,13]. However, PRC1 complex components are known to be highly variable and at least six distinct complexes have been described. Each complex is composed of a catalytic ubiquitin-ligase core containing RING1/RNF2 (also known as RING1A/B) and a Polycomb Group Ring Finger (PCGF) paralogue, bound to different accessory proteins[12,13]. The mode of recruitment of these distinct entities to chromatin, their effects on the regulation of gene expression and their distinct biological functions remain under investigation.

The BCOR protein was initially identified in a yeast two-hybrid screen as an interaction partner of the transcription factor (TF) B-cell Lymphoma 6 (BCL-6)[14] and was subsequently shown to be a member of a non-canonical PRC1 complex. In multiple different cell types BCOR co-purifies with RING1/RNF2, PCGF1, RYBP, SKP1 and the histone demethylase KDM2B, a complex commonly referred to as PRC1.1[15-18]. Interestingly, in germinal centre B cells, BCOR can assemble into an alternate CBX8 containing Polycomb complex[7] underscoring the context-dependent nature of PRC1 complexes. The C-terminal PCGF Ub-like fold discriminator (PUFD) domain of BCOR is necessary and sufficient for its interaction with PRC1 and can bind to PCGF1 or −3[19]. BCOR has also been shown to directly bind AF9, although whether this occurs outside of a PRC1 context is unknown[20].

Herein, a conditional mouse model of Bcor inactivation was developed to explore its function in myeloid differentiation and transformation. Using small numbers of haematopoietic cells isolated ex vivo, comprehensive analyses of the transcriptional and epigenetic consequences of Bcor loss were conducted. We demonstrate that Bcor has a pivotal role in the regulation of haematopoietic stem cell (HSC) associated transcriptional networks. Loss of Bcor results in expansion of myeloid progenitor cells, and in the context of oncogenic Kras, the initiation of cancer.

Results

Bcor mutation leads to expansion of myeloid progenitors

To analyse the function of Bcor in haematopoiesis in vivo a novel conditional mouse model was generated that mimics the truncating BCOR mutations observed in AML. We made use of a previously developed Bcor allele that has exons 9 and 10 flanked by LoxP sites (Fig. 1a) (Hamline et al., in prep.). Cre recombinase-mediated excision of exons 9 and 10 causes a frameshift and the introduction of a premature in-frame stop codon in the Bcor open reading frame. The recombined Bcor locus encodes a truncated Bcor protein that lacks the PUFD domain required for interaction with the PRC1 complex[19] (Fig. 1a).

Fig. 1

Bcor mutation increases repopulating capacity of myeloid progenitors. a Diagram of Bcor wild-type and conditional alleles (middle, exons 9–11 labelled) and the encoded truncated protein resulting from Cre-mediated recombination (bottom, protein domains labelled; B6BD, BCL-6 binding domain; Af9B, Af9 binding domain (present only in isoforms a and b) Ank Ankyrin repeat, PUFD PCGF ub-like fold discriminator, black box denotes 11 out-of-frame amino acids before stop codon). b PCR genotyping of sorted KL (Kit+Sca1−Lin−) cells performed 20 weeks post tamoxifen administration. PCR products corresponding to Bcor (red arrow), Bcor (green arrow), and Bcor (orange arrow) are shown. c Schematic overview of analysis of Bcor and Bcor mice. d Representative flow cytometry plot of bone marrow (BM) from Bcor and Bcor mice at 20 weeks. The strategy for gating KL (Kit+Sca1−Lin−) and KSL (Kit+Sca1+Lin−) cells within lineage negative (CD5−, B220−, CD11b−, Gr-1−, 7–4− and Ter-119−) BM cells is indicated. e Absolute numbers of indicated cell populations per femur quantified by flow cytometry (n = 3–4 mice/genotype; KL p = 0.0049, KSL p = not significant (N.S.)). f Schematic overview of transplantation experiments using total BM mononuclear cells or isolated BM KL cells. g Percentage of donor-derived cells (CD45.2+) within the indicated populations in the BM of Ptprca recipients transplanted with total BM mononuclear cells quantified by flow cytometry 20 weeks post-transplant (n = 8 recipients, KL: p = 0.00012; myeloid: p = 0.0075). h Absolute numbers and percentage chimerism of donor-derived cells in the peripheral blood and BM of Ptprca recipients transplanted with sorted KL cells quantified by flow cytometry (n = 9 recipients; PB week3, p = 2.084E-12, week 4, p = 3.84E-09, BM chimerism p = 4.0E-4, KL p = 1.2E-3, KSL p = 3.7E-3, Myeloid CD11b+only p = Myeloid CD11b+Gr1+ p = 2.3E-2). Values from individual animals are shown, bars indicate mean ±standard deviation, means compared with unpaired t-test: *p < 0.05, ** p <  0.005, *** p < 0.0005)

Bcor mice were crossed with HSC-SCL-Cre-ER mice to facilitate tamoxifen inducible Bcor deletion specifically in HSCs[21]. Administration of tamoxifen to a cohort of Bcor;HSC-SCL-Cre-ER animals resulted in efficient recombination at the Bcor locus (Fig. 1b, Supplementary Fig. 1A-B) (hereafter referred to as Bcor). Bcor mice were compared to tamoxifen treated single transgenic Bcor or HSC-SCL-Cre-ER littermates which expressed wild-type Bcor (hereafter referred to as Bcor). Using immunobotting we confirmed that Bcor results in low level expression of a truncated protein (Supplementary Fig. 1C), consistent with previously published findings of the same allele[8].

Bcor mice did not display overt haematopoietic dysfunction over a twenty-week observation period and their peripheral blood counts remained within normal limits (Supplementary Fig. 1D). To identify more subtle phenotypes in Bcor mice, a comprehensive immuno-phenotypic analysis of bone marrow populations was performed (Fig. 1c). There was a significant increase in the absolute number of cKit+Sca1−Lin− myeloid progenitor cells (KL cells) in the bone marrow Bcor mice, whilst there was no difference in the number of cells in the cKit+Sca1+Lin− HSC compartment (KSL cells) (Fig. 1d, e and Supplementary Fig. 1E, J). Further analysis of the sub-populations that make up the KSL compartment using a previously published gating strategy[22] (Supplementary Fig. 1F) demonstrated no significant difference in long term haematopoietic stem cells (LT-HSC), multi potent progenitor cells (MPP) or haematopoietic progenitor cells (HPC1/HPC2) (Supplementary Fig. 1G). The absence of an effect in KSL cells was not due to inefficient gene knockout, as comparable recombination of the Bcor locus was observed in both KL and KSL cells at 12 weeks (Supplementary Fig. 1A, B) post tamoxifen administration.

To further investigate the effect of Bcor mutation on haematopoiesis, competitive transplantation experiments were performed in congenic Ptprca mice (Fig. 1f). First, bone marrow cells from tamoxifen treated Bcor or Bcor donor mice were transplanted into sub-lethally irradiated Ptprca recipients and the level of leukocyte chimerism in the bone marrow was assessed after 20 weeks. Bcor cells were capable of tri-lineage differentiation, but also demonstrated an expansion of KL cells and an increase in mature myeloid cells compared with Bcor cells (Fig.1g, Supplementary Fig. 1H). As there was no difference between Bcor and Bcor KSL cells either at steady state or in the transplant setting, we hypothesised that the increased KL compartment in Bcor animals was due to an enhanced self-renewal or proliferative capacity of Bcor KL cells. To evaluate this, KL cells were isolated from Bcor or Bcor donor mice and transplanted into sub-lethally irradiated Ptprca recipients (Fig. 1f). The absolute number of donor-derived cells in the peripheral blood was tracked over 4 weeks after which point the bone marrow was assessed for donor-derived populations. Bcor KL cells displayed enhanced repopulating capacity, expanding to more than three times the size of the Bcor population in the peripheral blood of recipients over the 4-week period (Fig. 1h, Supplementary Fig. 1I). Further investigation revealed that the transplanted Bcor KL cells remained in the bone marrow after 4 weeks, whereas Bcor cells were largely depleted (Fig. 1h). Additionally, recipients transplanted with Bcor KL cells had higher numbers of donor-derived mature myeloid cells compared to those recipients that received Bcor KL cells (Fig. 1h, Supplementary Fig. 1I). Taken together, these data demonstrate that Bcor mutation in HSCs results in increased self-renewal and an expansion of myeloid progenitor cells in vivo.

Bcor regulates an HSC-associated transcriptional program

In order to discern the molecular mechanisms that underpin the aberrant phenotype of Bcor myeloid progenitors, gene expression analysis using 3′ RNA sequencing (Quant-seq) was performed on KL cells isolated from Bcor and Bcor mice. We identified 1497 significantly differentially expressed genes (DEG; see Supplementary Data 1) including many genes involved in regulation of haematopoiesis, leukocyte differentiation and activation (Fig. 2a and Supplementary Fig. 2A and B). Gene set enrichment analysis (GSEA) revealed that Bcor KL cells expressed higher levels of transcripts that are typically enriched in HSCs, and conversely expressed lower levels of transcripts that are upregulated during normal myeloid development (Fig. 2a, b). Quant-seq was also performed on donor-derived Bcor and Bcor KL cells isolated from Ptprca recipient mice. Comparison of the two independently derived data sets demonstrated a high degree of correlation between transcriptional changes in the transplant and primary transgenic settings (Supplementary Fig. 2C). These data suggest that Bcor is partially responsible for reprogramming the transcriptome during myeloid commitment and that Bcor mutation results in aberrant expression of a stem cell-like transcriptional program in myeloid progenitor cells.

Fig. 2

Bcor regulates master HSC TFs. a, b Quant-seq was performed on KL cells sorted from the BM of Bcor and Bcor mice at 20 weeks post tamoxifen administration. a Heatmap of normalized log counts per million (CPM) of the 500 most significantly differentially expressed genes. Genes associated with myeloid or HSC GSEA gene sets are indicated. b GSEA enrichment plots comparing Quant-seq from Bcor and Bcor KL cells against myeloid cell development and HSC gene sets. FDR False discovery rate, NES Normalized enrichment score. c–e Super-enhancer analysis based on H3K27ac ChIP-seq data in Bcor KL cells. c All enhancers identified in one representative Bcor sample ranked by H3K27ac tag count. Red dashed line indicates a super-enhancer slope of 1. Genes associated with super-enhancers identified in 2/3 Bcor samples and differentially expressed in Bcor KL cells compared with Bcor KL cells have been highlighted. d Venn diagram showing overlap of genes significantly differentially expressed in Bcor compared with Bcor KL cells (p < 0.05, absolute logFC >0.5), genes associated with SEs identified in 2/3 Bcor samples, and mouse transcription factors from RIKEN transcription factor database[68] (DEG Differentially expressed genes, TFs transcription factors, SE super-enhancer). e Heatmap of row-scaled Quant-seq counts for six differentially expressed, super-enhancer-associated genes encoding transcription factors. f H3K27ac ChIP-seq read density for the Hoxa7, Hoxa9 and Zfp36 loci in the indicated haematopoietic cell populations[25]. g Line graph showing normalized RNA expression for super-enhancer-associated Bcor-regulated transcription factor genes across in the indicated haematopoietic cell populations[25]

To identify genes that may mediate these effects, we focused on differentially expressed TFs, as these proteins play a key role in controlling cell-type-specific gene expression. In total 89 TFs were differentially expressed between Bcor and Bcor KL cells (Fig. 2b and Supplementary Table 4). These included Runx2, Myc and multiple members of the Gata, Cebp and Hoxa families with established roles in myeloid differentiation.

Recent studies have demonstrated that an interconnected network of TFs the expression of which is regulated by super-enhancers (SEs) plays a key role in the establishment of cell identify[23,24]. We reasoned that defining the SE network in Bcor KL cells would aid in the identification of key downstream mediators of Bcor function. To that end we generated a ChIP-seq dataset for H3K27ac, a mark associated with active chromatin that is widely used to determine the location of enhancers and SEs[24]. A total of 168 SEs were identified in at least two out of the three Bcor samples analysed (Fig. 2c, Supplementary 2E). In agreement with the existing literature[24], SE-associated genes were highly expressed (Supplementary Fig. 2D) and enriched for master haematopoietic TFs such as Tal1 and Fos. We identified six SE-associated TFs that were differentially expressed and all were upregulated in Bcor KL cells compared with Bcor KL cells (Fig. 2d, e). Analysis of publicly available H3K27ac ChIP-Seq and expression data from isolated mouse haematopoietic cell populations[25] demonstrated that each of the six genes were also associated with both SEs and high expression in HSCs. Moreover, their expression and SE-status was lost during differentiation (Fig. 2f, g and Supplementary Fig. 2E), suggesting that they are key drivers of HSC identity and that their downregulation is important for lineage commitment.

H2AK119ub is regulated during myeloid differentiation

Bcor is a component of PRC1.115 and inactivation of Bcor has been shown to impact on the activity of PRC1 and PRC2 at key loci[8,17]. This led us to initially characterise the genome-wide distribution of the PRC1-catalysed histone modification H2AK119ub and the PRC2-associated mark H3K27me3 in wild-type KL cells (Fig. 3a). To understand how these marks relate to other chromatin modifications and to chromatin accessibility, ChIP-seq for H3K4me3 (a modification associated with active transcription start sites) and Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq) were also performed. Analysis of the genomic distribution of H2AK119ub-enriched regions using a window-based method optimised for broad histone marks[26] demonstrated enrichment around promoter proximal regions of genes, with a lesser degree of occupancy at genomic loci distal to transcribed genomic regions (Supplementary Fig. 3A, D). Therefore, subsequent analysis was focussed on regions surrounding transcription start sites (TSS). K-means clustering was used to define five clusters based on ChIP and ATAC-seq results (Fig. 3a) and the transcriptional state of each of these clusters was assessed (Fig. 3b). H2AK119ub and H3K27me3 signal was highly enriched in clusters 1 and 2 and genomic loci enriched for these two marks were significantly overlapping when assessed by permutation test (p < 0.001; see Supplementary Table 5). These clusters also had low H3K4me3 signal and low chromatin accessibility as measured by ATAC-seq indicating a ‘closed’ conformation. In contrast clusters 3 and 4 had high levels of H3K4me3, ‘open’ chromatin and low levels of H2AK119ub and H3K27me3. While PRC1.1 has been associated with CpGs[27], no significant enrichment of CpG islands in H2AK119ub-containing clusters 1 or 2 compared to other clusters was identified. Furthermore, we found no evidence of H2AK119ub signal independent of H3K27me3 as has been previously suggested to be characteristic of PRC1.1 mediated H2AK119ub[10], highlighting that in these cells H2AK119ub appears to be predominantly co-occurring with H3K27me3. Gene expression analysis confirmed that genes in clusters 3 and 4 were highly expressed, whereas those in clusters 1 and 2 had lower expression (Fig. 3b), consistent with the active and repressive nature of their promoters, respectively.

Fig. 3

Bcor regulates H2AK119ub at a subset of PRC1 target genes. a Heatmap representing H2AK119ub, H3K27me3 and H3K4me3 occupancy, GpG island density (CpGi) and DNA accessibility (ATAC) at all annotated gene promoters (5 kb flanking TSSs of Refseq genes regardless of direction; scales indicated) in Bcor KL cells, separated into 5 groups by K-means clustering. b Tukey boxplot of RNA expression in Bcor KL cells for each of the 5 clusters identified in a. c Gene ontology analysis of genes associated with differentially ubiquitinated regions in Bcor KL cells compared with Bcor KSL cells (up in KL, log2FC >0.5; down in KL, log2FC <−0.5). d H2AK119ub ChIP-seq read density at Tal1 and Hoxa9 loci in Bcor KL and KSL cells. Differentially ubiquitinated regions are highlighted in red. e Tukey boxplot of normalized H2AK119ub ChIP-seq signal from top 100 differentially ubiquitinated regions in Bcor and Bcor KL cells (means compared with unpaired T-test, p = 2.03E-6.) f Venn diagram overlaying all genes associated with H2AK119ub-enriched regions in KL cells and genes associated with differentially ubiquitinated regions in Bcor KL cells compared with Bcor KL cells. g Hoxa9 and Bahcc1 gene loci showing H2AK119ub ChIP-seq read density in Bcor and Bcor KL cells. Differentially ubiquitinated region is highlighted in red. h Ranked plot comparing RNA expression between Bcor and Bcor KL cells. Genes that are significantly differentially expressed (log2FC >0.5, p < 0.05) and associated with decreased H2AK119ub (log2FC <−0.3, p < 0.01) in Bcor KL cells compared with Bcor KL cells are highlighted. i Venn diagram showing genes associated with H3K27me3, H2AK119ub or H3K4me3 differentially enriched regions (logFC >0.3 or <−0.3, p < 0.05) in Bcor KL cells compared with Bcor KL cells. Overlaps were assessed for significance by hypergeometric test[69] (H3K27me3 and H2AK119ub p = 3.4E-22, H3K27me3 and H3K4me3 p = 1.8E-37, H2AK119ub and H3K4me3 p = 0.00015, see Supplementary Table 5). j Dmwd and Six5 gene loci showing H2AK119ub, H3K27me3 and H3K4me3 ChIP-seq read density in Bcor and Bcor KL cells. One representative sample per genotype is shown. For Tukey boxplots centre line denotes median, whiskers extend to 1.5xIQR

To investigate whether deposition of the H2AK119ub mark is dynamically regulated during differentiation we compared ChIP-seq data derived from wild-type KSL and KL cells. Analysis of H2AK119ub ChIP-seq data identified 2118 H2AK119 ubiquitinated regions that were significantly altered during differentiation from KSL to KL cells. Association of the differentially ubiquitinated regions with the nearest gene identified 561 genes that were enriched for differentiation-related gene ontology (GO) terms such as development, cell fate commitment and regulation of haemopoiesis (Fig. 3c). Importantly, we observed a gain of H2AK119ub signal in KL cells at key TF loci Tal1, Hoxa7 and Hoxa9 (Fig. 3d). Taken together these results demonstrate that H2AK119ub is associated with gene repression and is redistributed during HSC differentiation including at multiple cell identity-defining TFs.

Bcor regulates H2AK119ub at a subset of PRC1 target genes

To understand how loss of Bcor impacts histone ubiquitination, the H2AK119ub distribution in Bcor and Bcor KL cells was compared. Analysis identified 9216 regions enriched for H2AK119ub 3-fold above input threshold in Bcor KL cells. Consistent with the hypothesis that PRC1.1 accounts for a fraction of overall PRC1 activity, we observed locus-specific changes rather than a global loss of H2AK119ub in Bcor KL cells. Analysis revealed 661 regions were significantly differentially ubiquitinated, predominantly within promoter regions of genes (Supplementary Fig. 3A). Notably, for the top 100 most differentially ubiquitinated loci, the majority were decreased in Bcor cells compared with controls (Fig. 3e and Supplementary Fig. 3B).

Next, we assigned each of the gene-proximal ubiquitinated regions in KL cells to the nearest gene. This analysis identified a total of 2757 H2AK119ub-associated genes, of which 115 were significantly altered in Bcor cells (Fig. 3f, Supplementary Data 2), including Hoxa7, Hoxa9 and Bahcc1 (Fig. 3g). Despite the small number of genes associated with differential ubiquitination, GO analysis demonstrated that they were highly enriched for hematopoietic processes including blood cell differentiation and cell fate commitment (Supplementary Fig. 3C). Loss of H2AK119ub was significantly associated with increased gene expression as assessed by hypergeometric testing (p =  0.012, Supplementary Table 5) and we identified 10 putative direct targets of PRC1.1., where loss of H2AK119ub was accompanied by upregulation of transcription (Fig. 3h). Of these identified target genes, six have known roles in transcriptional regulation, signalling or proliferation in cancer (Hoxa9[28], Hoxa7[29], Prex1[30], Clic4[31], Zfp36l1[32] and Calm2[33]), and are therefore potentially contributing to the transcriptional reprogramming and expansion of Bcor myeloid progenitor cells.

H2AK119ub deposited by non-canonical PRC1 complexes has previously been shown to drive recruitment and activity of PRC2[34]. In order to investigate whether Bcor mutation was associated with alterations in H3K27me3 or H3K4me3, further ChIP-seq in Bcor and Bcor KL cells was performed. Our analysis revealed no global change in H3K27me3 or H3K4me3 signal (Supplementary Fig. 3D), although as expected we did identify locus-specific differences (Fig. 3i). Interestingly, we observed a significant decrease in H3K27me3 signal in the top 100 differentially methylated regions, but this was not the case for regions where H2AK119ub was decreased in Bcor cells (Supplementary Fig. 3E). There was however, a significant overlap between genes associated with H2AK119ub changes and those associated with a change in H3K27me3 alterations as assessed by hypergeometric test (p = 3.4E-22, Fig. 3i, Supplementary Table 5). The changes in H3K4me3 were balanced (i.e. increased and decreased in equal measures) and there were no changes in signal either in the top 100 differentially methylated regions or in those regions where H2AK119ub was decreased in Bcor cells (Supplementary Fig. 3E). We did identify seven loci, including Hoxa7, Hoxa9, Dmwd and Six5, associated with alterations in all three marks (centre, Fig. 3i). These genes uniformly exhibited significant loss of both H2AK119ub and H3K27me3, and a gain of H3K4me3 (Fig. 3j), indicative of epigenetic activation. Alterations in H3K27me3 and H3K4me3 signal were significantly correlated with changes in mRNA expression (Supplementary Fig. 3F), suggesting these epigenetic changes may be indicative of consequential transcriptional changes downstream of Bcor.

BcorΔE9-10 and KrasG12D co-operate to initiate leukaemia

To determine whether the phenotypic and molecular changes driven by Bcor inactivation can contribute to malignant transformation we aged a cohort of Bcor and control mice for over 12 months following tamoxifen injection. During this period, the animals did not develop leukaemia suggesting that inactivation of Bcor alone was not sufficient for cancer initiation (Fig. 4a). Mutations in epigenetic modifiers that alter haematopoietic differentiation have been shown to co-operate with dysregulated growth factor signalling to drive oncogenic transformation[22,35] and BCOR mutations co-occur with RAS pathway mutations in both AML and MDS patients[2,36,37]. Thus, Bcor mice were crossed with Lox-stop-lox-Kras mice that express oncogenic Kras following excision of the transcriptional stop cassette by Cre (subsequently referred to as Kras)[38]. After somatic recombination of both alleles in HSCs, Kras and BcorKras mice developed a lethal haematopoietic disease characterised by leucocytosis, splenomegaly and increased leukaemic blasts in the peripheral blood and bone marrow (Fig. 4a, b, Supplementary Fig. 4A). Notably however, the survival of Bcor Kras mice was greatly reduced compared with Kras controls (Fig. 4a). Bone marrow cells from both moribund Kras and BcorKras animals could recapitulate a lethal disease in primary and secondary recipients upon retransplantation (Supplementary Fig. 4B). Flow cytometry analysis revealed a significant expansion of the KL cell compartment in the bone marrow and spleen of BcorKras compared to Kras animals (Fig. 4b), whilst mature lineage populations in both spleen and bone marrow were similar (Supplementary Fig. 4C). Our data formally demonstrates that Bcor is a tumour suppressor in the myeloid lineage and can co-operate with additional events to initiate and accelerate leukaemia in vivo.

Fig. 4

Bcor cooperates with Kras to initiate leukaemia. a Kaplan–Meier survival curve of primary transgenic mice following tamoxifen administration. (median survival for BcorKras mice was 58 days compared with 113 days for Kras mice; survival compared with logrank Mantel–Cox test, p < 0.0001; median survival not reached for Bcor and Bcor groups). b Percentage of KL and KSL populations in the bone marrow and spleen of moribund BcorKras (n = 5) and Kras only (n =5) mice quantified by flow cytometry (means were compared using unpaired t-test BM KL p = 0.02, Spleen KSL p = 0.027 KL p = 0.024). c Peripheral white blood cell (WBC) counts and spleen weights of moribund transgenic mice (WBC: BcorKras n = 9, Kras n = 4; spleen: BcorKras n = 12, Kras n = 5). d Heatmap showing correlation in gene expression changes regulated by Bcor in the non-malignant and malignant contexts. The top 250 most significantly differentially expressed genes (comparing BcorKras and Kras KL cells) are included. Correlation analysis was performed with Spearman’s rank correlation (ρ = 0.69, p-value <2.2E-16). e GSEA enrichment plots comparing Quant-seq from BcorKras and Kras murine leukaemia against gene expression signatures derived from BCOR-mutant AML patients (see Methods). f Tukey boxplots of normalized expression of selected target genes in AML patient samples. g Tukey boxplot of normalized H2AK119ub ChIP-seq tag counts from the top 100 differentially ubiquitinated regions comparing BcorKras and Kras KL cells (means compared with unpaired t-test, p = 3.9E-05). h Zfp36l1 and Hoxa9 gene loci showing H2AK119ub ChIP-seq read density in BcorKras and Kras KL cells. Differentially ubiquitinated regions are highlighted in red. i Ranked plot comparing RNA expression between BcorKras and Kras KL cells. Genes that are significantly differentially expressed (log2FC >0.5 or <−0.5, p < 0.05) and associated with decreased H2AK119ub (log2FC <−0.3, p < 0.05) in BcorKras KL cells compared with Kras cells are highlighted. Values from individual animals are shown, bars indicate mean ± standard deviation. For Tukey boxplots centre line denotes median, whiskers extend to 1.5xIQR

Bcor regulates H2AK119ub and gene expression in leukaemia

To establish whether Bcor regulates expression of the same genes in a leukaemic context expression of Hoxa5, 7 and 9 was assessed by qPCR in bone marrow from moribund animals. Upregulation of Hoxa5, 7 and 9 was confirmed in BcorKras compound mice compared with Kras controls (Supplementary Fig. 4D) suggesting that the transcriptional consequences of Bcor may be conserved in the malignant setting.

To expand on these findings CD45.2+ bone marrow cells from donor BcorKras and Kras mice that had comparable disease burden were transplanted into irradiated Ptprca recipients. Donor-derived KL cells were isolated and analysed by Quant-seq. Differential gene expression analysis identified 434 significantly altered genes between Kras and BcorKras cells that were highly enriched for differentiation and development-associated GO terms (Supplementary Fig. 4F, see Supplementary Data 3). GSEA analysis revealed that a pattern specification gene set strongly associated with Hox genes was upregulated in BcorKras cells (Supplementary Fig. 4E). Changes in gene expression following Bcor mutation in the non-transformed and malignant settings were significantly correlated (Fig. 4d, Supplementary Fig. 4G), with key TFs Hoxa7, Hoxa9 and Tal1 also upregulated in the context of leukaemia (Fig. 4d).

To assess the  relevance of the BcorKras murine model to human AML, we utilised two previously published mRNA expression data sets[39-41] and analysed differential expression between BCOR-mutant (BCOR) and BCOR wild-type (BCOR) leukaemia patients. Patient data sets #1 and #2 contained 518 and 437 AML patient samples with 13 and 19 BCOR patients, respectively (BCOR frequency ~3%). GSEA analysis revealed that BCOR signatures derived independently from the two data sets were strongly enriched in BcorKras cells compared with Kras cells (Fig. 4d). A number of relevant genes were identified as upregulated by BCOR mutation in both the mouse model and human patient data including NAV1 and H1F0 (Fig. 4e). HOXA family genes were not upregulated in BCOR patients, however, these genes are known to be upregulated by multiple independent mechanisms in AML[28,42], thus introducing a confounding factor to our analysis. Overall, these results suggest that Bcor is a faithful model of BCOR AML.

We next assessed the distribution of H2AK119ub in leukaemic KL cells using ChIP-seq. We identified 9221 H2AK119ub-enriched regions, 1103 of which were differentially ubiquitinated between BcorKras and Kras cells (See Supplementary Data 4). Consistent with our earlier findings, for the top 100 most differential regions there was a significant loss of H2KA119ub signal in BcorKras cells (Fig. 4g). We assigned genes to H2AK119ub-enriched regions and identified 285 genes associated with differential ubiquitination (Fig. 4h, Supplementary Fig. 4H), including 17 genes associated with reduced ubiquitination in both the normal and leukaemic contexts (Supplementary Data 2 and 4). GO analysis of genes associated with differential ubiquitination revealed terms linked to embryonic morphogenesis, differentiation and proliferation (Supplementary Fig. 4I). H2AK119ub loss was significantly associated with increased gene expression when assessed by both hypergeometric testing (p = 9E-11) and GSEA (p = 0.0165; Supplementary Fig. 4J and Supplementary Table 5). We identified twenty-three putative PRC1.1 target genes that were associated with both a loss of H2AK119ub and upregulation of transcription in BcorKras cells (Fig. 4i). Importantly, this included Hoxa7, Hoxa9 and Zfp36l1 that were also identified as putative PRC1.1 targets in the non-leukaemic context. Altogether, these data indicate that Bcor-dependent loss of histone ubiquitination and transcriptional regulation is largely conserved in both normal and malignant haematopoiesis.

Hoxa9 is essential for proliferation of BcorΔE9-10KrasG12D cells

To determine the impact of Bcor target genes on the proliferation of BcorKras leukaemic cells a CRISPR drop-out screen was performed. We focussed on a subset of Bcor-regulated genes at which H2AK119ub loss was accompanied by increased transcription in BcorKras cells. These genes were selected because they represent putative direct targets of PRC1.1. Primary BcorKras tumour cells were extracted from the spleen of a moribund mouse and transduced with Cas9 followed by our custom sgRNA library (see Supplementary Data 5). Samples were collected post transduction and then after 24 and 42 days of in vitro culture, and guide representation determined by sequencing (Fig. 5a). Statistical analysis of sgRNAs that were significantly depleted at day 24 compared to day 0 revealed that guides targeting Hoxa9 were significantly depleted (p = 0.04737, Fig. 5b, c, Supplementary Data 6). Whilst not statistically significant, the next highest scoring genes were Rhoq, Hoxa7 and Hoxb5 (Fig. 5b, c, Supplementary Data 6), indicating these genes may also be important for BcorKras leukaemic cell survival. Hox family genes have been implicated in leukaemia previously[43] and are therefore likely to be key drivers of Bcor phenotype. Guides targeting Dnmt3a were enriched at the end time point, indicating the upregulation of this gene is not essential to the Bcor phenotype and may even be antagonistic to some degree (Supplementary Fig. 5B, C). Accordingly, previous studies have demonstrated loss of Dnmt3a enhances self-renewal in HSCs[44]. Furthermore, these results were recapitulated when cells were harvested at 42 days (Supplementary Fig. 5D). Through these additional analyses, we can confirm that Hoxa9 upregulation is essential for the leukaemic phenotype of BcorKras tumours, and Hoxa7 and Hoxb5 may also be important.

Fig. 5

Hoxa9 is essential for the proliferation of BcorKras leukaemic cells. a Schematic overview of the CRISPR drop-out screen. b Negative enrichment of all genes in the CRISPR screen calculated by MAGeCK[65] (negative enrichment for Hoxa9 = 0.014268, p = 0.04737). c Selected sgRNA counts at T0, T24 and T42. d Model for Bcor function in haematopoietic differentiation and the initiation of leukaemia

Discussion

BCOR is recurrently mutated in AML and a recent study implicated Bcor in the proliferation of myeloid cells in vitro[5]. The data presented here demonstrates an important role for Bcor in regulating myeloid differentiation and leukaemogenesis in vivo. Loss of Bcor function resulted in expansion of myeloid progenitor cells in Bcor mice and competitive transplantation experiments confirmed that Bcor-mutant myeloid progenitors outcompeted their wild-type counterparts. Notably, isolated Bcor KL cells exhibited enhanced repopulating capacity in vivo, suggesting that Bcor may function to restrict self-renewal in lineage committed myeloid progenitors.

Consistent with the observed phenotype, RNA sequencing revealed that ex vivo sorted Bcor KL cells were characterised by an HSC-like transcriptional signature. Loss of Bcor function was accompanied by broad gene expression changes and ~1500 Bcor-regulated genes were identified. Of particular note, we observed upregulation of six SE-associated TFs: Hoxa7, Hoxa9, Lmo2, Fos, Tal1 and Zfp36. Each of these genes have well-established roles in both normal haematopoiesis and the initiation and maintenance of leukaemia[28,29,45-51]. However, we cannot exclude the contribution of other DEGs to the Bcor phenotype.

Bcor is a known component of the PRC1.1 complex[15,16,52] and has been suggested to modulate H2AK119ub at key loci in myeloid cells[8]. This led us to examine the impact of Bcor loss on the distribution of H2AK119ub in KL cells. The catalytic components of PRC1, RING1 and RNF2, are essential for normal and leukaemic cells, implying that H2AK119ub is required for cell survival[53]. In contrast, altered expression of other PRC1 members, such as upregulation of BMI1 and CBX7 and downregulation of MEL18, are known to play a role in malignant progression[54-58], suggesting that altered PRC1 function may promote leukaemogenesis. Consistent with the proposed role of BCOR as a member of the PRC1.1 complex[15,19], Bcor KL cells had reduced levels of H2AK119ub at a small subset of ubiquitin-associated genes. Loss of H2AK119ub was significantly associated with increased transcription in malignant and non-transformed Bcor KL cells and we identified putative direct targets of PRC1.1 including Hoxa7, Hoxa9 and Zfp36l1. The contribution of Hoxa7 and 9 to cancer are well established[29,59] and Zfp36l1, which encodes butyrate response factor 1, has been implicated in myeloid differentiation and negative regulation of the anti-apoptotic protein BCL2[32,60]. It is tempting to speculate that these genes are the initial drivers of transcriptional reprogramming in Bcor-mutant KL cells and indeed we found that Hoxa9 is indispensable for the proliferation of BcorKras leukaemic cells. However, Bcor-driven differential gene expression cannot solely be explained by PCR1.1 activity and we also identified key haematopoietic factors such as Tal1, Lmo2 and Fos that were not associated with H2AK119ub changes. Altered expression of these genes in Bcor cells may be the result of Bcor transcriptional regulation that is independent of the ubiquitin-ligase activity of PRC1.1 or the downstream effects of primary Bcor targets.

The interaction between PRC1 and 2 and the requirement for catalytic activity to mediate transcriptional repression at target genes is not well defined. Consistent with previous studies[10] we found a robust genome-wide association between H2AK119ub and H3K27me3 at proximal promoters of target genes suggesting that the two complexes predominantly act in tandem. There was a significant overlap between genes associated with changes in H2AK119ub and those with changes in H3K27me3. However, actual genomic regions at which H2A ubiquitination was altered in Bcor cells, did not exhibit a significant difference in H3K27 methylation. Interestingly, there were also many loci where H3K27me3 loss was not accompanied by a reduction in H2AK119ub signal. Whether this is directly or indirectly mediated by Bcor is unclear, but one intriguing possibility is that the PRC1.1 complex has an H2AK119ub-independent impact on PRC2 recruitment or activity[17].

Functional studies across multiple systems suggest that BCOR is capable of both oncogenic and tumour suppressive functions, even in closely related cell types[7-10,61]. Here we generated a clinically-relevant model and showed that Bcor inactivation resulted in functional co-operation with oncogenic Kras to initiate leukaemia in vivo. Interestingly, a recent study that used an internal deletion mutant of Bcor found that removal of the BCL6-binding domain also has an impact on haematopoietic differentiation and the development of T-cell leukaemia[61]. While specific deletion or inactivation of the BCL6-binding domain has not been described in human cancer to date, we cannot exclude the possibility that Bcor has additional important functions independent of histone ubiquitination or polycomb binding. As protein expression from the Bcor allele is severely reduced, these additional functions are likely to be abrogated in our model as well, and thus also contribute to the phenotype we observe.

In summary, our work identifies BCOR as an important regulator of normal and malignant haematopoiesis. BCOR controls the expression of key HSC TFs, in part through ubiquitin-mediated silencing during myeloid differentiation. Accordingly, loss of Bcor function results in enhanced self-renewal of myeloid progenitors and in the context of other oncogenes, promotes malignant progression.

Methods

Animals

The Peter MacCallum Cancer Centre Animal Ethics Committee approved all in vivo procedures in this study. Complete details of the generation of the Bcor allele will be described elsewhere (Hamline et al. in prep). In short, loxP sites flanking exons 9 and 10 of the Bcor gene were inserted into the mouse genome. Cre-mediated excision of exons 9 and 10 results in a frameshift and the introduction of a premature stop codon in the Bcor open reading frame. HSC-SCL-Cre-ER mice (also known as Tg(Tal1-cre/ERT)42–056Jrg) and Lox-stop-lox-Kras mice (also known as Kras) were obtained from Jackson Labs. To activate Cre recombinase 3–6-month-old mice were dosed with a total of 5 mg tamoxifen suspended in corn oil (1% tamoxifen, 10% ethanol, 89% corn oil) via intra-peritoneal (IP) injection over five consecutive days. Both male and female animals were used for immunophenotyping, whereas only male animals were used for molecular analysis (including Quant-seq, ChIP-seq and ATAC-seq). Congenic C57BL/6.SJL-Ptprca (referred to throughout the text as Ptprca, a congenic C57BL/6 strain that expresses the differential pan leukocyte marker CD45.1) were purchased from the Animal Resources Centre of Peter MacCallum Cancer Centre or the Walter and Eliza Hall Institute of Medical Research (Melbourne, Australia). All experimental mice were maintained under specific pathogen-free conditions.

Genotyping

Mice were initially genotyped at three weeks of age and additional analyses performed at endpoin. Genomic DNA was extracted from ear clippings, tail clippings bone marrow or spleen cells using DNeasy Blood & Tissue Kit (Qiagen). Genotyping was performed with gene specific primers (Supplementary Table 3) and Phusion DNA polymerase (Thermo Fisher) according to manufacturer’s instructions. Recombination of the Bcor locus was quantified by qPCR using the SYBRgreen High-ROX mix (Agilent) on a LightCycler LC480 (Roche) with primers internal to the deletion site (Bcor exon 9); P primers targeting Bcor exon 6 were used for normalization (Supplementary Table 3).

Flow cytometry analysis

Single-cell suspensions of whole blood, bone marrow and spleen were incubated in ACK red cell lysis buffer (150mM NH4Cl, 10mM KHCO3, 0.1mM EDTA) to lyse erythrocytes and isolate nucleated cells. Cells were resuspended in FACS buffer (2% foetal calf serum (FCS; [Thermo Fisher] in phosphate buffered saline[PBS]) and stained with fluorophore-conjugated antibodies targeted against cell surface markers (Supplementary Table 1) on ice for 15 min. Cells were then washed with FACS buffer and immediately analysed on BD LSR Fortessa X-20 or BD LSR II flow cytometers (BD).

Isolation of myeloid progenitor cells

Bone marrow was harvested from male transgenic mice 20 weeks post tamoxifen administration. Bone marrow cells were incubated with ACK red cell lysis buffer to remove erythrocytes and depletion of lineage positive cells was performed using the MACS Lineage Cell Depletion kit for mouse (Miltenyi Biotech, Cat#130–090–858) and magnetic LS columns (Miltenyi Biotech, Cat#130–042–401) according to the manufacturer’s instructions. Lineage depleted cells were stained with fluorophore-conjugated antibodies targeted against cell surface markers (Supplementary Table 1) as above and 4′,6-Diamidino-2-Phenylindole (DAPI) was used as a viability marker. Viable KL (cKit+Sca1−Lin−) were sorted on BD FACSAria Fusion 5, BD FACSAria II and BD FACSAria Fusion 3 (BD Biosciences) and used for transplantation and molecular analyses.

Cell culture and lentiviral transduction

Sorted KL cells or single-cell suspensions from the spleens of tumour bearing BcorKras mice (referred to as BK cells) were cultured in Dulbecco's Modified Eagle Medium (DMEM) supplemented with penicillin/streptomycin, 20% FCS, L-Asparagine (0.1mM), IL-3 (2ng/ml), IL-6 (2ng/ml), and SCF (10ng/ml). Cells were maintained at 37 °C and 10% CO2. Cells were genotyped as above to confirm their identity and were tested for Mycoplasma contamination. For transduction of BK cells, lentiviral supernatants were harvested from packaging 293T cells transfected with the custom CRISPR library (see below). Lentiviral particles were centrifuged onto RetroNectin-coated (Takara) non-treated tissue-culture plates and virus infection performed according the manufacturer’s protocol.

In vivo competitive transplantation assays

Female Ptprca recipients were sub-lethally irradiated (5.5Gy administered as a single dose) using an X-RAD 320 irradiator (Precision X-Ray Incorporated) and sorted KL cells (5 × 104) or bone marrow mononuclear cells (5 × 105) from donor animals were inoculated via tail vein injection. Chimerism in peripheral blood and bone marrow was assessed by flow cytometry using the markers CD45.1 (recipient) and CD45.2 (donor).

Western blotting

Cells were lysed in Laemmli buffer containing β-mercaptoethanol, subjected to SDS–PAGE on 4–20% Mini-PROTEAN TGX precast polyacrylamide gels (Bio-Rad Cat#4561093) and transferred to PVDF membrane (Millipore). Membranes were incubated with primary Bcor or β-actin antibodies (see Supplementary Table 2 for details) followed by anti-rabbit or anti-mouse IgG HRP-conjugated antibodies (Dako) Blots were developed with enhanced chemo-luminescence (ECL) reagents (GE Healthcare Cat#RPN2106) and exposed to film.

In vitro CRISPR screen

Target genes were identified by integrating Quant-seq and H2AK119ub ChIP-seq data using the following selection criteria: increased mRNA expression (log2FC >0.5 and p < 0.05) and decreased H2AK119ub in the proximal promoter (log2FC <−0.3 and p < 0.05) in BcorKras cells compared with Kras cells. We designed a custom CRISPR library targeting 25 genes (4 guides per gene) with 10 non-targeting control guides (guide sequences are provided in Supplementary Table 5). Guide sequences were PCR amplified from a CustomArray Inc oligo pool and cloned into the lentiGuide-Crimson lentiviral vector (Addgene plasmid #70683) using Golden Gate cloning[62]. The representation of all guides in the pool was validated by sequencing on a MiSeq (Illumina) with single-end 50 bp sequencing. BK cells were transduced using FUCas9Cherry (Addgene plasmid #70182) and Cherry+ cells were sorted on a BD FACSAria Fusion 5 (BD Biosciences). 8 × 106 BKCas9 cells were transduced with a multiplicity of infection of <0.1 to achieve single sgRNA integration per cell at an average 500-fold representation. Five days post viral infection Crimson+ cells were sorted on a BD FACSAria Fusion 5 (BD Biosciences) and returned to culture. Samples were then collected at four days (T0), 28 days (T24) and 46 days (T42) post sorting and genomic DNA extracted using DNeasy Blood & Tissue Kit (Qiagen). Sequencing libraries were prepared by nested PCR method as described previously[63] and were sequenced to a depth of 106 reads with single-end 75 bp sequencing on a NextSeq (Illumina). Sequencing reads were demultiplexed using bcl2fastq (v2.17.1.14) and low-quality reads Q < 30 were removed. The reads were trimmed using cutadapt (v1.14)[64] to extract the 20 bp targeting sequence and sgRNAs that were depleted relative to T0 determined using MAGeCK (v0.5.7)[65]. ggplot2 (v2.2.1) was used for figure generation.

Quantitative RT-PCR

RNA was extracted from whole bone marrow using Trizol reagent and Directzol RNA miniprep kit (Zymo research) according to manufacturer’s instructions. cDNA was synthesised using High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher) as per manufacturer’s instructions. qPCR was performed with SYBRgreen High-ROX mix (Agilent) on a LightCycler LC480 (Roche). Actin was used as a housekeeping control gene. Primer sequences are listed in Supplementary Table 3.

Quant-seq

1–5 × 105 sorted KL cells were lysed in Trizol reagent and RNA was extracted using Directzol RNA miniprep kit (Zymo research) according to manufacturer’s instructions. The Quant-seq 3′ mRNA-seq Library Prep Kit for Illumina (Lexogen) was used to generate libraries as per the manufacturer’s instructions. Libraries were pooled and sequenced with 75 bp single-end sequencing to a depth of 6–10 × 106 reads on a NextSeq (Illumina). Sequencing reads were demultiplexed using bcl2fastq (v2.17.1.14) and low-quality reads Q < 30 were removed. The RNA sequencing reads were trimmed at the 5′ end using cutadapt (v1.14)[66] to remove bias introduced by random primers used in the library preparation and 3′ end trimming was performed to eliminate poly-A-tail derived reads. Reads were mapped to the reference genome (mm10) using HISAT2. Reads were counted using subread software package (v1.5.0-p3)[63]. Differential gene expression analysis was performed using R package LIMMA (v3.32.4). R packages pheatmap (v1.0.8) and ggplot2 (v2.2.1) were used for figure generation. GSEA (v3.0) was used for analysing enrichment of gene sets and GO analysis was performed using metascape [http://metascape.org].

ChIP-seq

KL and KSL cells were sorted from Bcor (n = 3), Bcor (n = 2), BcorKras (n = 3) and Kras (n = 3) mice. 1–5 × 105 cells were crosslinked with 1% formaldehyde and 5% FCS in PBS at room temperature for 10 min, then the reactions were quenched with 1.25 M glycine for 5 min. Crosslinked cells were washed with ice cold dilution buffer (20 mM Tris pH 8, 150 mM NaCl,  2mM EDTA, 1%Triton-X) then resuspended in nuclear extraction buffer (50 mM TRIS pH7.5, 150 mM NaCl, 5 mM EDTA, 0.5%NP-40, 1% TritonX-100). Cells were washed in dilution buffer again before being lysed in ChIP lysis buffer (20 mM TRIS pH7.5, 150 mM NaCl, 0.5M EDTA, 1%NP-40, 0.3% SDS). Chromatin was fragmented to ~150–300 bp by sonicating lysates in a Covaris Ultrasonicator using 100 μl microtubules (Covaris). Immunoprecipitation was performed overnight at 4 °C in 50% dilution and 50% lysis buffer with magnetic protein A and G Dynabeads (Life technologies) pre-coupled to 1μg antibody per sample (antibodies listed in Supplementary Table 2). Beads were then washed with wash buffer 1 (20 mM TRIS pH 8, 500 mM NaCl, 2 mM EDTA, 1% Triton-X, 0.1% SDS) followed by wash buffer 2 (220 mM TRIS pH 8, 250 mM LiCl, 2 mM EDTA, 0.5% NP-40, 0.5% deoxycholate) and then TE buffer (10 mM TRIS pH7.5, 1 mM EDTA) on ice. Beads were then resuspended in reverse crosslinking buffer (1% SDS, 100 mM NaHCO3, 200 mM NaCL, 300 μg/mL proteinase K) and incubated at 55 °C for 30 min followed by 65 °C for 30 min. DNA was then purified using Zymo ChIP clean and concentrate kit according to manufacturer’s instructions (Zymo Research) and quantified using Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific). Indexed libraries were prepared using KAPA Hyper Prep Kit for Illumina platforms (Kapa Biosystems) and the SeqCap Adapter Kit (Roche) following vendor’s instructions. Library QC and quantification was performed using D1000 high-sensitivity screen tape with 4200 TapeStation Instrument (Agilent Technologies) and size selected for between 200 bp and 500 bp using a Pippin Prep system (Sage Science). Libraries were pooled and sequenced with 75 bp single-end sequencing to a depth of 15–20 × 106 reads per sample on a NextSeq (Illumina). bcl2fastq (v2.17.1.14) was used for de-multiplexing. The Fastq files generated by sequencing were aligned to the mouse reference genome (GRCm38/mm10) using bowtie2 (v2.3.3). Samtools (v1.4.1) was used for manipulation of SAM and BAM files. MACS (v2.1.1) was used for traditional peak-calling. R package Csaw (v1.10.0) was used to quantify regions enriched for histone marks, which were then associated with the closest TSS using Bedtools (v2.26). Superenhnacer analysis was performed using Homer (v4.8). Browser viewable TDF files were generated using IGVTools (v2.3.95) and ChIP-Seq tracks were visualized using IGV (v2.3.55). Graphics were generated using deeptools (v2.5.3) and R packages pheatmap (v1.0.8), plotly (v4.7.1) and ggplot2 (v2.2.1).

ATAC-seq

ATAC-seq was performed using a previously published protocol[67]. In brief, 5 × 104 sorted KL cells were lysed and nuclei washed, followed by chromatin tagmentation using TD1 enzyme and TDE buffer supplied in Illumina Nextera DNA prep kit (cat#FC-121–1030). Chromatin was tagmented at 37 °C for 30 min, then purified using the Minielute kit (Qiagen). A DNA library was prepared for sequencing by ligating adapters and amplifying using KAPA hot start polymerase (KAPA biosystems) according to manufacturer’s instructions then purifying using Minielute kit (Qiagen). QC was performed using D1000 high-sensitivity screen tape with 4200 TapeStation Instrument (Agilent Technologies) before size selection for fragments between 150 bp and 700 bp using a Pippin Prep System (Sage Science). Libraries were pooled and sequenced with 75bp single-end sequencing to a depth of 15–20 × 106 reads per sample on a NextSeq (Illumina). bcl2fastq (v2.17.1.14) was used for de-multiplexing. The Fastq files generated by sequencing were aligned to the mouse reference genome (GRCm38/mm10) using bowtie2 (v2.3.3). Samtools (v1.4.1) was used for manipulation of SAM and BAM files. Data were visualised with deeptools (v2.5.3) alongside ChIP-seq data.

Analysis of publically available data sets

H3K27ac and RNAseq BigWig files (From GSE60103[25]) were downloaded using SRA toolkit (v2.5.1), visualised in IgV (v2.4.5) and analysed alongside our ChIP-seq and Quant-seq data as above. GpG island data were downloaded from UCSC table browser. Data were converted using bedtools (v2.26) and ucsc-utils (v20161019) then analysed and visualised alongside our data with deeptools (v2.5.3) and IgV (v2.4.5).

For patient dataset #1 the series matrix file from Verhaak et al. (2009; GEO accession GSE6891[41]) was downloaded using the GEOquery R package (v2.42.0) and probe expression values were collapsed to median values per gene. Samples that had not been exome-sequenced were then removed from the analysis. For patient dataset #2 raw RNA counts and series matrix files were downloaded for the Leucegene project[39] (GEO accessions GSE49642, GSE52656, GSE62190, GSE66917 and GSE67039). Counts were normalized with R using limma-voom (v3.32.10) and batch effects were removed using ComBat with the sva R package (v3.24.4). Differential expression analysis was performed using limma and results visualised using ggplot2 (v2.2.1). GSEA(v3.0) was used to compare to our RNA sequencing data above.