BET inhibition prevents aberrant RUNX1 and ERG transcription in STAG2 mutant leukaemia cells

Mutations in the subunits of the cohesin complex, particularly in the STAG2 subunit, have been identified in a range of myeloid malignancies, but it is unclear how these mutations progress leukaemia. Here, we created isogenic K562 erythromyeloid leukaemia cells with and without the known leukemic STAG2 null mutation, R614*. STAG2 null cells acquired stem cell and extracellular matrix gene expression signatures that accompanied an adherent phenotype. Chromatin accessibility was dramatically altered in STAG2 null K562 cells, consistent with gene expression changes. Enhanced chromatin accessibility was observed at genes encoding hematopoietic transcription factors, ERG and RUNX1. Upon phorbol 12-myristate 13-acetate (PMA)-induced megakaryocytic differentiation, STAG2-null cells showed precocious spike in RUNX1 transcription from its P2 promoter. A similar precocious spike was observed in transcription of ERG. Interestingly, spikes in RUNX1-P2 and ERG only occurred as immediate early response to differentiation induction. Treatment of STAG2 null cells with enhancer-blocking BET inhibitor, JQ1, dampened precocious RUNX1 P2 expression and led to a complete loss of RUNX1 P1 and ERG transcription during PMA stimulation in both parental and STAG2 null K562 cells. These results suggest that precocious RUNX1 and ERG expression in STAG2 null cells is enhancer-driven. Furthermore, JQ1 treatment reduced stem cell-associated KIT expression in STAG2 null cells. We conclude that STAG2 depletion in leukemic cells amplifies an enhancer-driven transcriptional response to differentiation signals, and this characteristic is dampened by BET inhibition. The results have relevance to the development of therapeutic strategies for myeloid leukaemia

Read counts were retrieved by exon and summarized by gene using featureCount version v1.5.3. Differentially expressed genes in the STAG2-null lines versus parental (WT) were identified using DESeq2 (Love et al., 2014). P-values were adjusted for multi test using Independent Hypothesis Weighting (Ignatiadis et al., 2016). Euclidean distance matrix correlation was computed in R, based on clustering of biological replicates for each cell type.
GSEA software (Subramanian et al., 2005) (settings; enrichment statistics: weighted; max size: 5000) was used to carry out enrichment analyses on ranked gene lists based on Wald statistics from DESeq2 on MgSigDB gene sets. The ranked lists contained genes significant at adjusted p-values ≤0.05 from the DESeq2 analyses.

ATAC sequencing and analyses
50,000 cells were transposed as described previously (Buenrostro et al., 2015;Litzenburger et al., 2017). Libraries from three biological replicates of parental K562 (WT) and STAG2-nullA were prepared using Illumina® Nextera DNA library preparation (FC-121-1030) and indexing (FC-121-1011) kits, and sequenced at 100 bp single end reads on HiSeq 2500 at the Otago Genomics Facility (New Zealand). ATAC-sequencing reads were aligned to the human genome GRCh37 (hg19) using bowtie2. Peaks were called using Macs2 without model and shift size to accommodate transposase integration (--nomodel --shift -37 --extsize 73). Peak summits were extended ±250 bp. Peaks from all biological replicates were merged and normalized. Only peaks overlapping in all 3 replicates were considered. Biological replicates correlation matrix was computed using Euclidean distance.

ATAC sequencing differential analyses
Differential peaks were identified using edgeR (McCarthy et al., 2009). Briefly, peak count data were retrieved using featureCounts (Smyth et al., 2013). Differential events were estimated using Generalized Linear Models. After likelihood test ratio, p-values were adjusted for multiple testing using Benjamini-Hochberg FDR <0.05.

ATAC sequencing peak annotation
Peaks were annotated to genomic sites and genes using HOMER with Hg19 annotation (Heinz et al., 2010). Enrichment of differential ATAC peaks at super enhancers (dbsuper database; (Khan and Zhang, 2016)) was determined using the 'intersectBED' function from BED tools utilities (Quinlan and Hall, 2010).

ATAC sequencing motif analyses
Enrichment of transcription factor binding motifs were analysed using HOMER for known motifs (Heinz et al., 2010).

Array CGH
Array comparative genomic hybridization was performed using the Agilent SurePrint G3 ISCA CGH+SNP 4x180k(hg19) array platform at the Wellington Genetic Services Facility (New Zealand). The parental K562 was used as the reference for array CGH.

Flow cytometry
Cells treated with 500 nM JQ1 or DMSO were stained with PE anti-CD117 (KIT) (340529, BD Biosciences) and APC anti-CD15 (551376, BD Biosciences), flow cytometry was carried out on the Navios EX flow cytometer at Southern Community Laboratories (New Zealand) and analysed using Kaluza Software (Beckman Coulter). Cell cycle progression was determined following double thymidine block and propidium iodide staining as previously described (Antony et al., 2015;Dasgupta et al., 2016) and analysed using Flow-JO software.

Statistical analyses
Unless otherwise stated all statistical analyses were carried out using Prism version 7 software (GraphPad).
Supplementary Data 1. Differentially regulated transcripts (p≤0.05) in STAG2-null lines compared to parental wild type K562 cells.

Supplementary Data 2.
Array CGH analyses of STAG2-null lines. The tables list the regions with gain or loss in STAG2-null lines compared to parental wild type K562 cells.

Supplementary Data 3. Differentially accessible ATAC peaks (p≤0.05) identified in
STAG2-null lines compared to parental wild type K562 cells.

Supplementary Data 4.
List of super enhancers (SE) (described in K562, CD14+ and CD34+ cells) and associated transcripts that have differential accessibility (ATAC) and or differential expression in STAG2-nullA compared to parental K562.

A.
B.