Glucocorticoid response leads to pervasive effects on the ALL chromatin landscape
The cellular response to GCs was mapped in two human B-ALL cell lines (697 and Nalm6) using diverse functional genomic assays over a time-course spanning a 24-hour window (0hr, 2hr, 6hr, 12hr and 24hr timepoints; schematic in Supplementary Fig. 1) following treatment with prednisolone13. Importantly, cell viability was not affected during the time-course.
We assessed for GC-responsive epigenomic changes using ATAC-seq and ChIP-seq for the histone H3 lysine 27 acetylation (H3K27ac) post-translational modification, both of which mark active cis-regulatory elements34–36. In total, 12597 and 16984 changes to chromatin accessibility, as well as 29073 and 28854 alterations in H3K27ac enrichment were uncovered in 697 and Nalm6 cells, respectively across all time points (Fig. 1A). Intriguingly, these two chromatin features demonstrated opposing temporal patterns, with most changes to acetylation occurring early within 6hrs of GC treatment, whereas most accessibility changes occurred after 6hrs (Supplementary Fig. 2).
Because GC treatment leads to changes in chromatin state, we wanted to further determine if GCs promote the formation of super-enhancers (SEs). Rank Ordering of Super-Enhancers (ROSE)37, 38 using H3K27ac data identified 278 (697) and 254 (Nalm6) reproducible GC-responsive SEs following GC treatment present at two or more timepoints. The vast majority were de novo SEs that were identified after GC treatment, whereas a small subset was pre-established SEs that show > 2-fold H3K27ac enrichment following GC treatment. (Fig. 1B). Notably, the top GC-responsive SE in each cell line was located within the ZBTB16 gene locus, a GC response gene identified in vivo39 (Fig. 1C,D).
Glucocorticoid receptor (GR) ChIP-seq over the time-course identified that 54.9% (697) and 58.1% (Nalm6) GR occupancy sites harbored a discernable GR motif, consistent with previous studies40 (Fig. 1E, Supplementary Fig. 3). Notably, most GC-responsive accessible chromatin (61%) and H3K27ac (52%) sites harbored GR occupancy (Supplementary Fig. 4), including all GC-responsive SEs. To identify high-confidence sites of GR occupancy mapping to cis-regulatory elements, we intersected GR sites with ATAC-seq and H3K27ac peaks at each time point (GR + ATAC-seq + H3K27ac sites; henceforth referred to as HGRs) and uncovered 26096 and 29058 HGRs in 697 and Nalm6 cells, respectively. Motif analyses uncovered that HGRs with canonical GRE motifs preferentially harbor enhanced chromatin accessibility and H3K27ac enrichment following GC treatment compared to HGR without GREs (Fisher’s Exact p < 2.2x10− 16; Fig. 1F), consistent with a role in transcriptional activation.
RNA-seq (0hr, 6hr and 24hr) identified 3414 and 2767 differentially expressed genes (DEGs, FDR < 0.01, absolute fold change > 2; combined DEGs from 6hr and 24hr) in 697 and Nalm6 cells, respectively. Supporting a role for GR in regulating the expression of many of these GC-responsive genes, HGRs with enhanced H3K27ac following GC treatment were enriched near up-regulated genes (K-S test p < 2.2x10− 16; Supplementary Fig. 5,6). By contrast, HGRs with reduced H3K27ac following GC treatment were enriched near down-regulated genes (K-S test p < 7.9x10− 10). To more directly link HGRs and GC-responsive chromatin alterations with transcriptional effects we performed H3K27ac HiChIP41 to map long-range, three-dimensional interactions between distal cis-regulatory elements and promoters (0hr, 6hr and 24hr). Most loops (FitHiChIP q < 0.01) were identified after 24hr of GC treatment, and an average of 37.8% (697) and 31.7% (Nalm6) of all loops were interactions between distal cis-regulatory elements and promoters (Supplementary Fig. 7). Notably, 14% and 18% (697) as well as 13% and 28% (Nalm6) of 6hr and 24hr HiChIP promoter loops mapped to DEGs (absolute fold change > 2) at 6hrs and 24hrs respectively, and 65% (697) and 88% (Nalm6) of these loops, on average, involved HGRs (Fig. 1G). GC-responsive H3K27ac alterations further correlated with HiChIP promoter looping at GC-responsive genes. For all loops to DEG promoters (absolute fold change > 2), GC-responsive chromatin sites exhibiting enhanced H3K27ac preferentially looped to up-regulated gene promoters compared to down-regulated gene promoters (697 p < 2.2x10− 16, 6hr DEG odds ratio = 10, 24hr DEG odds ratio = 5.2; Nalm6 p < 8.6x10− 8, 6hr DEG odds ratio = 2.3, 24hr DEG odds ratio = 4.7; Fisher’s Exact), while GC-responsive chromatin sites exhibiting reduced H3K27ac preferentially looped to down-regulated gene promoters compared to up-regulated gene promoters (697 p < 2.9x10− 11, 6hr DEG odds ratio = 6.3, 24hr DEG odds ratio = 7.2; Nalm6 p < 2.2x10− 16, 6hr DEG odds ratio = 13.8, 24hr DEG odds ratio = 7.1; Fisher’s Exact). Taken together, these data indicate that GCs elicit pervasive changes to chromatin state, including the formation of GC-responsive SEs, and these chromatin alterations drive transcriptional programs in ALL cells through long-range looping at sites of GR occupancy.
Identification of transcription factors impacted by GC responses
Transcription factor (TF) footprints42 at accessible chromatin sites harboring GR occupancy were integrated with RNA-seq to identify candidate TFs that cooperate with GR and identified numerous ETS-family TFs (Supplementary Fig. 8,9,10). TF footprints were further used to identify GC-induced changes in TF occupancy (Fig. 2A,B). As expected, the GR motif (represented by NR3C1) was enriched following GC treatment. However, a strong depletion for AP-1 footprints was also observed. Concordant with these findings, accessible chromatin sites with AP-1 footprints were significantly enriched at GC-responsive sites with reduced accessibility following GC treatment compared to accessible chromatin sites devoid of AP-1 footprints (Fisher’s Exact p < 2.2x10− 16, 697 odds ratio = 5.8, Nalm6 odds ratio = 5). Multiple AP-1 TFs were also down-regulated following GC treatment (Fig. 2C,D), and HiChIP loops were uncovered between distal HGRs and the promoters of FOS, FOSB and JDP2 (697) or FOSB and JUN (Nalm6), supporting a direct role for GR in AP-1 transcriptional repression (Fig. 2E).
To explore these data further, PECA statistical analysis43 using chromatin accessibility and gene expression was used to infer changes in TF-gene connections after GC treatment. Validating our findings, GR showed enhanced gene connectivity (Fig. 2F; represented by NR3C1), and numerous AP-1 TFs exhibited reduced connectivity. Moreover, ETS-family TFs exhibited enhanced connectivity, supporting their role as GR cooperating TFs. Using GREAT44, we further determined that GC-responsive AP-1 bound sites were associated with genes involved in cell proliferation and anti-apoptotic processes (Fig. 2G), and concordant pathways were uncovered when limiting associated genes to down-regulated genes (fold change < 2 or < 1.5). Collectively these data suggest that the repression of AP-1 bound cis-regulatory elements and genes precedes GC-induced apoptosis
Identification of transcriptionally active sequences using ATAC-STARR-seq
Self-transcribing active regulatory region sequencing45, 46 at accessible chromatin sites (ATAC-STARR-seq) high-throughput reporter assays were performed in 697 and Nalm6 cells to functionally validate the activity of GC-responsive chromatin sites and HGRs (Fig. 3A). Notably, 98% (697) and 97% (Nalm6) of accessible chromatin sites were cloned into plasmids (Fig. 3B,C). Following transfection of ATAC-STARR-seq plasmid libraries and GC treatment (0hr, 6hr and 24hr), RNA output libraries were sequenced to over one billion reads and STARR-seq active sites were identified using BasicSTARRseq (https://git.bioconductor.org/packages/BasicSTARRseq). On average, we identified 8864 (697) and 3101 (Nalm6) active sites (p < 0.001) at each timepoint. Because active sites at 6hr and 24hr GC treatment showed substantial overlap and read count correlation (Fig. 3D,E), they were combined for all downstream analyses (henceforth named GC; Fig. 3F). Over 66% (697) or 83% (Nalm6) of STARR-seq active sites mapped to HGRs. Notably, 697 and Nalm6 GC-responsive open chromatin sites were significantly enriched at STARR-seq active sites compared to accessible chromatin sites that were not GC-responsive (Fisher’s Exact p < 5.4x1010). However, GC-responsive chromatin sites exhibiting enhanced accessibility where significantly more enriched compared to sites exhibiting reduced accessibility (Chi-square p < 2.2x10− 16; 2.2–3.5 fold; Fig. 3G,H,I).
Shared STARR-seq active sites and sites specific to 0hr or to GC treatment were further examined. Shared active sites between 0hr and GC treatment exhibited substantially greater overlap with promoter locations compared to 0hr- and GC-specific sites, which had a higher proportion of promoter-distal elements (Fig. 3J, Supplementary Fig. 11). GREs were significantly more enriched at GC-specific active sites in 697 cells compared to 0hr-specific active sites (Chi-square p = 1.3x10− 4, 1.5-fold) and exhibited a strong trend for greater enrichment in Nalm6 cells (Chi-square p = 0.077, 1.6-fold). GC-specific active sites were also more closely associated with GC-responsive upregulated genes compared to 0hr-specific sites (24hr DEGs, fold change > 2; K-S test p = 0.05 [697] and p = 0.002 [Nalm6]). Overall, ATAC-STARR-seq validated thousands of GC-responsive accessible chromatin sites, and most harbor GR occupancy.
Genetic disruptions to the GC response impact GC resistance in patient samples
Resistance of primary ALL cells to GCs is predictive of treatment response in patients measured as either persistence of minimal residual disease after remission induction treatment or overall treatment outcome4–6, 10, 47, 48. Thus, ex vivo measurements of primary ALL cell resistance to GCs is concordant with in vivo resistance and predicts treatment outcome in patients. We therefore investigated the impact of inherited genetic variants at HGRs that were associated with ex vivo GC resistance in primary cells.
Using previously published genotyping and ex vivo GC drug sensitivity data in primary ALL cells from patients enrolled in St. Jude Total Therapy XVI17, 18, we intersected variants associated with resistance to prednisolone and/or dexamethasone and variants in high linkage disequilibrium (r2 > 0.8) with HGRs in 697 and Nalm6 cells and fine-mapped 45 variants to HGRs (Supplementary Table 1), including several that were expression quantitative trait loci (eQTLs) for six genes previously implicated in GC resistance (ARHGAP18, ATG10, BFSP2, PPM1E, TLE1 and XRRA1)18. A notable variant was rs7045812 (C/T) which mapped to a Nalm6 HGR and altered a GRE (Fig. 4A). This intragenic variant is located within the TLE1 gene locus, a GC-responsive up-regulated gene that has also been previously correlated with prognostic features in ALL patients49 and functions with Groucho as a canonical Wnt signaling repressor50, 51. In line with the recognized consensus GRE, the alternative T allele, which disrupts the GRE, is associated with greater GC resistance in primary cells (Fig. 4B). In support of these data, luciferase reporter assays testing 300-bp DNA fragments centered on rs7045812 confirmed that the alternative T allele negatively impacted GC-responsiveness compared to the reference C allele (Fig. 4C, Supplementary Fig. 12).
HiChIP further identified a long-range promoter interaction between rs7045812 and the TLE1 promoter (Fig. 4D). To validate these findings, CRISPR/Cas9 genome editing was used to delete the TLE1 HGR in Nalm6 cells (Fig. 4E). Analysis of TLE1 HGR knockout cells determined that TLE1 expression was significantly reduced compared to wild-type cells, thereby validating the regulatory effects of this HGR on TLE1 expression (Fig. 4F). We further tested the impact of TLE1 disruption on GC resistance using CRISPR/Cas9 knockout in Nalm6 cells (Supplementary Fig. 13). GC drug response assays using prednisolone uncovered that TLE1 disruption led to greater GC resistance (Fig. 4G, Supplementary Fig. 14), and this was concordant with lower TLE1 expression being associated with GC resistance in primary cells (Fig. 4H). Collectively, these data highlight a role for inherited genetic variation at sites of GR occupancy in GC resistance.
Integrative analysis identifies epigenetic disruptions to the GC response impacting GC resistance in patient samples
Epigenetic disruptions to HGRs were also explored to better understand their impact on GC resistance. Using published ATAC-seq data in primary ALL cells from our laboratory52, ex vivo GC drug sensitivity assays were performed on 19 of these primary ALL cells. We stratified cells by GC resistance (Fig. 5A) and identified 1929 sites where differential chromatin accessibility associated with GC resistance (FDR < 0.1; Fig. 5B, Supplementary Table 2; henceforth referred to as GC-resistance accessible chromatin sites). To understand factors contributing to chromatin accessibility differences we examined DNA methylation that was available for a subset of patient samples. Of the 908 GC-resistance accessible chromatin sites that overlapped with 2566 CpG probes, 85% of probes exhibited patterns consistent with chromatin accessibility differences between GC-sensitive and GC-resistant samples, highlighting DNA methylation as a contributor to chromatin alterations (Fig. 5C). We also uncovered that GC-resistance accessible chromatin sites were significantly enriched at HGRs compared to accessible chromatin sites not associated with GC-resistance (Fisher’s Exact p < 2.2x10− 16, odds ratio = 1.73). Most of these GC-resistance accessible chromatin sites (78%) exhibited greater occlusion in GC-resistant primary cells, concordant with the role of GCs as a chemotherapeutic drug promoting apoptosis (Fig. 5D). To identify top candidate HGRs impacting GC resistance in patient samples, we performed an integrative analysis that combined GC-resistance accessible chromatin sites with GC response maps, genes implicated in GC resistance18 and CRISPR interference (CRISPRi) screening (Fig. 5E).
CRISPRi using Enhancer-i53 was used to screen all GR occupancy sites at GC-resistance accessible chromatin sites for GC resistance phenotypes using 11038 sgRNAs and 100 non-targeting control sgRNAs (Fig. 5F and Supplementary Table 3). Following a 72-hr drug selection with prednisolone, we identified numerous sgRNAs exhibiting significant enrichment (697 = 1844, Nalm6 = 774; FDR < 0.05; Supplementary Table 4,5), whereas control sgRNAs overall did not show strong or preferential enrichment (Fig. 5G). By further aggregating all sgRNA data at each GR site using MAGeCK33 (~ 6.5 sgRNAs per site), we ranked GR sites by log2-transformed fold change enrichment (Supplementary Table 6,7). As expected, control sgRNAs were situated near the center of the ranking and did not exhibit strong sgRNA enrichment or depletion (Supplementary Fig. 15). We next identified CRISPRi-enriched GR sites that mapped to HGRs (log2 fold change > 0), determined if these HGRs were associated with GC-responsive genes using GREAT44 and examined if these genes had been implicated in GC resistance18. This integrative analysis identified 35 top HGRs associated with 26 GC-responsive genes implicated in GC resistance (Fig. 5H, Supplementary Table 8).
Functional evaluation of epigenetically disrupted HGRs
Closer functional examinations of top HGRs were performed to further associate epigenetic cis-regulatory disruptions in patient samples with GC resistance. We identified a top HGR (Peak1585) in 697 cells that was situated downstream of TLE1 (Fig. 6A). Because we identified TLE1 as a GC-resistance gene harboring an intragenic HGR variant associated with GC resistance (Fig. 4), we functionally investigated this distal HGR. Greater chromatin occlusion was observed in GC-resistant primary cells, and this was supported by greater CpG DNA methylation in GC-resistant samples (Fig. 6B). CRISPR/Cas9-mediated knockout of this HGR significantly reduced TLE1 expression (Fig. 6C,D), supporting the functional role of this distal HGR in TLE1 regulation.
The ROR1 locus was identified as the top hit because it contained (i) three top HGRs (Peak42, Peak43 and Peak44), (ii) an HGR (Peak42) exhibiting HiChIP looping to the ROR1 promoter and (ii) an HGR (Peak42) with STARR-seq activity in 697 cells (Fig. 6E). ROR1 is a receptor tyrosine kinase-like orphan receptor for Wnt5a that can induce activation of non-canonical Wnt signaling54 and has been associated with survival of TCF3-PBX ALL55. All three HGRs exhibited greater chromatin occlusion in GC-resistant samples and these alterations were further supported by greater DNA methylation in GC-resistant samples (Fig. 6F,G,H). CRISPR/Cas9 disruption of ROR1 led to greater GC resistance (Fig. 6I, Supplementary Fig. 16,17). In line with these data, lower baseline ROR1 expression was associated with greater GC resistance in primary cells (Fig. 6J). Although ROR1 is a GC-responsive gene in 697 cells (RNA-seq log2 fold change=-0.61 and Supplementary Fig. 18), CRISPR/Cas9-mediated deletions of the two distal HGRs resulted in ablation of GC-induced ROR1 repression (Fig. 6K,L,M,N). Overall, these data uncovered cis-regulatory epigenetic disruptions to GC responses as a mechanism impacting GC resistance in ALL.