Abstract
Glucocorticoids (GCs) are a mainstay of contemporary, multidrug chemotherapy in the treatment of childhood acute lymphoblastic leukemia (ALL), and resistance to GCs remains a major clinical concern. Resistance to GCs is predictive of ALL relapse and poor clinical outcome, and therefore represents a major hurdle limiting further improvements in survival rates. While advances have been made in identifying genes implicated in GC resistance, there remains an insufficient understanding of the impact of cis-regulatory disruptions in resistance. To address this, we mapped the gene regulatory response to GCs in two ALL cell lines using functional genomics and high-throughput reporter assays and identified thousands of GC-responsive changes to chromatin state, including the formation of over 250 GC-responsive super-enhancers and a depletion of AP-1 bound cis-regulatory elements implicated in cell proliferation and anti-apoptotic processes. By integrating our GC response maps with genetic and epigenetic datasets in primary ALL cells from patients, we further uncovered cis-regulatory disruptions at GC-responsive genes that impact GC resistance in childhood ALL. Overall, these data indicate that GCs initiate pervasive effects on the leukemia epigenome, and that alterations to the GC gene regulatory network contribute to GC resistance.
Similar content being viewed by others
Introduction
Acute lymphoblastic leukemia (ALL) is the most common cancer in children, and represents 25% of all childhood cancers [1]. Despite treatment advances, relapsed ALL remains the fifth most common pediatric cancer [1]. Unlike newly diagnosed ALL, relapsed ALL has a low overall survival of only 40% [2], and is commonly characterized by chemotherapeutic drug resistance [3]. Consequently, understanding mechanisms of antileukemic drug resistance is critical for improving overall survival in childhood ALL.
Glucocorticoids (GCs) are a mainstay of contemporary multidrug chemotherapy in ALL and resistance to GCs is predictive of ALL relapse and poor clinical outcome [4,5,6,7,8,9,10,11]. GCs are unique chemotherapeutic agents because they function through activation of the glucocorticoid receptor (GR) transcription factor (TF), which is encoded by the NR3C1 gene [12]. Once activated by GCs, GR translocates into the nucleus where it binds to glucocorticoid response element (GRE) sequences or forms protein–protein interactions at cis-regulatory elements to drive a transcriptional program that leads to leukemic cell apoptosis [12, 13]. Although the exact molecular mechanism is not fully understood, reduced glucose metabolism [14] and metabolic reprogramming [15], or the suppression of B-cell development genes [16] has been suggested.
Previous studies have identified transcriptional signatures and genes implicated in GC resistance in ALL [10, 16,17,18,19]. Up-regulation of the NALP3 inflammasome was associated with CASP1-mediated cleavage of GR, leading to GC resistance [17], and loss of GR expression through mutation was also uncovered as a mechanism of GC resistance [20]. Another mechanism involves the BTG1 gene, which is commonly deleted in pediatric ALL and impacts GR autoinduction [21]. Transcriptional co-factors are also implicated in GC resistance. CREBBP and NCOR1 are genes frequently mutated at relapse in ALL, and CREBBP mutations impair the activation of GC-responsive genes [22]. TBL1XR1 gene repression is also observed at relapse and impacts GR chromatin recruitment [23], and decreased expression of SWI/SNF subunits was associated with GC-resistance [24]. Phosphorylation and inactivation of the EHMT1/2 coregulator complex via up-regulation of AURKB at relapse were also shown to negatively impact GC-induced activation of genes that drive cell apoptosis [19].
Because the mechanism of action for GCs involves activation of the GR TF and the impact of transcriptional co-factors on GC resistance, disruption of GR binding sites would also be predicted to impact GC resistance. In support of this, a recent study mapped lymphocyte-specific chromatin accessibility in ALL cells and identified chromatin alterations that correlate with GC resistance, including a GR site within the BIM gene locus that harbored enhanced DNA methylation in GC-resistant cells [25]. Apart from this study, there has otherwise been limited investigation on how alterations to GR binding impact GC sensitivity in ALL.
To better understand the gene regulatory responses to GCs in ALL, and the impact of cis-regulatory disruptions on GC resistance, we mapped epigenomic responses to GCs in two ALL cell lines (697 and Nalm6) and integrated these maps with genomic data in primary ALL cells from patients. These two ALL cell lines were chosen for technical and biological reasons. They grow well in culture and cellular transfections are more efficient in these cell lines compared to other established ALL cell lines. They are also widely used in ALL research studies and are representative of ALL molecular subtypes that are common in children (697 = TCF3-PBX1 and Nalm6 = DUX4/ERG). Notably, both cell lines are sensitive to GCs, thereby allowing the mapping of gene regulatory effects leading to cellular apoptosis. Because of epigenetic and gene regulatory heterogeneity among diverse ALL subtypes [26], the use of two cell lines that are representative of distinct subtypes additionally provides for a more comprehensive evaluation of the ALL cell GC response. Collectively, our study delineated extensive GC-mediated effects on the chromatin landscape and identified genetic and epigenetic cis-regulatory disruptions of GC-responsive genes as a mechanism impacting GC resistance in childhood ALL.
Materials/subjects and methods
Patient samples
Written informed consent was obtained from all patients or their legal guardians. Gene expression, DNA methylation, SNV genotyping and ex vivo drug sensitivity data [17, 18] were collected as part of St. Jude Total Therapy XVI (NCT00549848) [27]. The use of these samples was approved by the institutional review board at St. Jude Children’s Research Hospital. Leukemia blasts were isolated from bone marrow obtained prior to treatment and subject to Ficoll gradient centrifugation. Samples underwent further enrichment by magnetic-activated cell sorting if blast percent was <85%.
Functional genomics assays
Fast-ATAC was performed [28] on 10,000 fresh cells from ALL cell lines. For RNA-seq, total RNA was purified from patient samples using the Norgen Total RNA Purification Kit (Norgen, 35300). ChIP-seq using anti-GR-alpha (BD, 611227) and anti-H3K27ac (Active Motif, 39133) antibodies was performed as previously described [29, 30]. Fast-ATAC data from primary ALL cells were downloaded from NCBI Gene Expression Omnibus (GSE161501). Additional details are provided in Supplementary Methods.
ATAC-STARR-seq
Fast-ATAC transposed DNA from 697 and Nalm6 cells was cloned into the hSTARR-seq_ORI vector (Addgene plasmid #99296). Cells were transfected and treated with vehicle control or prednisolone. Additional details are provided in Supplementary Methods.
Western blot
Cells were plated 1 × 106 cells/mL and treated with prednisolone (697 = 10 µM; Nalm6 = 5 µM) for 24 h or vehicle control. Cells were collected, washed with PBS and lysed in RIPA Buffer with protease and phosphatase inhibitors. Protein (40µg) was loaded and run on 4–12% Bis-Tris gradient gels, transferred onto PVDF membrane, blocked with 5% Milk in TBST and incubated with primary antibodies overnight at 4 °C. The following antibodies were used: TLE1 (Abcam, ab183742), ROR1 (Cell Signaling Technologies, 4102S), Actin (Sigma, A5441). Membranes were washed with TBST and incubated with HRP-conjugated secondary antibody.
CRISPR/Cas9 genome editing
Gene and HGR knockouts were generated using CRISPR/Cas9 genome editing after transient transfection with precomplexed ribonuclear proteins. Additional details are provided in Supplementary Methods.
CRISPR interference screens
All sgRNAs and control non-targeting sgRNAs were designed and synthesized by Custom Array, amplified, purified (Qiagen PCR Purification Kit, #28104), cloned into the lentiviral pXPR_003-puro-IRES-CFP vector using Gibson Assembly Cloning Kit (NEB, #E5510S) and packaged into viral particles. Additional details are provided in Supplementary Methods.
Drug viability assays
Drug viability assays were performed in 96-well plates and treated for 72-h (ALL cell lines) or 96-h (primary ALL cells) with prednisolone as previously outlined [17, 18]. Additional details are provided in Supplementary Methods.
Luciferase reporter assays
A 300-bp of sequence centered on reference or the alternative allele of rs7045812 was cloned upstream of the minimal promoter in pGL4.23 (Promega, E841A). Nalm6 cells were transfected with constructs using the Neon transfection system (Thermo Fisher, MPK5000). Following an overnight incubation after transfection, cells were treated with 5uM prednisolone or vehicle control for 6 h before luciferase was measured using the Dual Luciferase Reporter Assay System (Promega, E1960).
Data analysis
All NGS reads were mapped to the hg19 reference genome using bowtie2 [31] and ATAC-seq and ChIP-seq peaks were identified using MACS2 using default parameters [32]. Differentially enriched sites and expressed genes were identified using DESeq2 [33]. sgRNA enrichments were determined using DESeq2 [33] and aggregate sgRNA log2 fold changes were calculated using MAGeCK [34]. Additional details are provided in Supplementary Methods.
Results
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 h window (0, 2, 6, 12, and 24 h timepoints; schematic in Supplementary Fig. 1) following treatment with prednisolone [13]. 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 elements [35,36,37]. 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 cell lines, respectively across all timepoints (FDR < 0.01; Fig. 1A). These two chromatin features also demonstrated opposing temporal patterns, with most changes to acetylation occurring early within 6 h of GC treatment, whereas most accessibility changes occurred after 6 h (Supplementary Fig. 2). Notably, the extent of GC-responsive chromatin alterations was in line with previous investigations in the non-leukemic A549 cell line [38] (Supplementary Fig. 3). However, we identified substantial cell type-specificity, with only 5.8% (ATAC-seq) and 8.3% (H3K27ac) of GC-responsive sites shared, on average, between leukemic and non-leukemic cell lines, which is consistent with our previous findings demonstrating low nuclear receptor binding concordance between distinct cell types [39].
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) [40, 41] 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 showing >2-fold H3K27ac enrichment following GC treatment (Supplementary Fig. 4). Notably, the top GC-responsive SE in each cell line was located within the ZBTB16 gene locus, a GC response gene identified in vivo [42] (Fig. 1B).
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 studies [39] (Supplementary Fig. 5). Notably, most GC-responsive accessible chromatin (61%) and H3K27ac (52%) sites harbored GR occupancy (Supplementary Fig. 6), 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 timepoint (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.2 × 10−16; Supplementary Fig. 7), consistent with a role in transcriptional activation.
RNA-seq (0, 6, and 24 h) identified 3414 and 2767 differentially expressed genes (DEGs, FDR < 0.01, absolute fold change >2; combined DEGs from 6 and 24 h) 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 upregulated genes (K-S test p < 2.2 × 10−16; Supplementary Figs. 8, 9). By contrast, HGRs with reduced H3K27ac following GC treatment were enriched near downregulated genes (K-S test p < 7.9 × 10−10). To more directly link HGRs and GC-responsive chromatin alterations with transcriptional effects we performed H3K27ac HiChIP [43] to map long-range, three-dimensional interactions between distal cis-regulatory elements and promoters (0, 6, and 24 h). Most loops (FitHiChIP q < 0.01) were identified after 24 h 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. 10). Notably, 14% and 18% (697) as well as 13% and 28% (Nalm6) of 6 and 24 h HiChIP promoter loops mapped to DEGs (absolute fold change >2) at 6 and 24 h respectively, and 65% (697) and 88% (Nalm6) of these loops, on average, involved HGRs (Fig. 1C). 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 upregulated gene promoters compared to downregulated gene promoters (697 p < 2.2 × 10−16, 6 hr DEG odds ratio = 10, 24 h DEG odds ratio = 5.2; Nalm6 p < 8.6 × 10−8, 6 h DEG odds ratio = 2.3, 24 h DEG odds ratio = 4.7; Fisher’s Exact), while GC-responsive chromatin sites exhibiting reduced H3K27ac preferentially looped to downregulated gene promoters compared to upregulated gene promoters (697 p < 2.9 × 10−11, 6 h DEG odds ratio = 6.3, 24 h DEG odds ratio = 7.2; Nalm6 p < 2.2 × 10−16, 6 h DEG odds ratio = 13.8, 24 h 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 cell lines through long-range looping at sites of GR occupancy.
Identification of transcription factors impacted by GC responses
Transcription factor (TF) footprints [44] at accessible chromatin sites harboring GR occupancy were integrated with RNA-seq to identify candidate TFs that cooperated with GR and identified numerous ETS-family TFs (Supplementary Figs. 11–13). TF footprints were further used to identify GC-induced changes in TF occupancy (Fig. 2A, Supplementary Figs. 14, 15). As expected, the GR motif (represented by NR3C1) was enriched following GC treatment. However, a strong depletion of 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.2 × 10−16, 697 odds ratio = 5.8, Nalm6 odds ratio = 5). Multiple AP-1 TFs were also downregulated following GC treatment (Fig. 2B), 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. 2C).
To explore these data further, PECA statistical analysis [45] 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 (Supplementary Fig. 16), and numerous AP-1 TFs exhibited reduced connectivity. Moreover, ETS-family TFs exhibited enhanced connectivity, supporting their role as GR cooperating TFs. Using GREAT [46], we further determined that GC-responsive AP-1 bound sites were associated with genes involved in cell proliferation and anti-apoptotic processes (Supplementary Fig. 17), and concordant pathways were uncovered when limiting associated genes to downregulated 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 sequencing [47, 48] at accessible chromatin sites (ATAC-STARR-seq) were performed in 697 and Nalm6 cell lines to functionally validate the activity of GC-responsive chromatin sites and HGRs (Fig. 3A). However, this assay can identify activity for all cloned ATAC-seq accessible chromatin sites, irrespective if they are occupied by GR or GC-responsive. To control for effects from variable DNA sequence length, ATAC-seq transposed DNA was size selected prior to cloning into the STARR-seq plasmid. As validation of our plasmid library, we determined that 98% (697) and 97% (Nalm6) of accessible chromatin sites were cloned into plasmids (Supplementary Fig. 18). Cellular transfections in each ALL cell line were performed using 4 bulk transfections that were subsequently split to create the cell populations for the 3 different treatment conditions (GCs for 0, 6, or 24 h; n = 4 replicates per treatment condition in each cell line) to control for variation in transfection efficiencies between treatment conditions. Following transfection of ATAC-STARR-seq plasmid libraries and GC treatment (for 0, 6, or 24 h), RNA output libraries from replicate treatment conditions were sequenced and merged to increase total sequence depth and coverage. RNA output libraries were also downsampled to control for read count differences between treatment conditions and active STARR-seq sites at 0, 6, and 24 h were identified by significant increases in RNA output versus DNA input using BasicSTARRseq (https://git.bioconductor.org/packages/BasicSTARRseq).
On average, we identified 8864 (697) and 3101 (Nalm6) active STARR-seq sites (p < 0.001) at each timepoint (examples in Fig. 2B). As additional verification of these findings, we implemented a second approach to identify active STARR-seq sites using DESeq2 [33] differential analysis that compared enrichment at accessible chromatin sites of ATAC-STARR-seq RNA output replicates (n = 4) over ATAC-seq pooled input replicates (n = 3) that consisted of ATAC-seq replicates (#1, 2 or 3) pooled across the time course (e.g., pooled replicate #1 for 0, 2, 6, 12, and 24 h, etc.; Supplementary Methods). Supporting this approach, ATAC-STARR-seq plasmid DNA input was highly correlated with pooled ATAC-seq signal (697 r2 = 0.95, Nalm6 r2 = 0.94; Supplementary Fig. 19). Consistently, nearly all genomic regions identified as active STARR-seq sites using the DESeq2-based approach (FDR < 0.05) were also identified by the BasicSTARRseq approach (697 = 95.9%, Nalm6 = 88.4%), validating our BasicSTARRseq findings and supporting the overall reproducibility and robustness of our ATAC-STARR-seq results.
Because active STARR-seq sites after 6 and 24 h of GC treatment showed substantial overlap and read count correlation (Supplementary Fig. 20), they were combined for all downstream analyses (henceforth named GC-active STARR-seq sites; Fig. 3C). Over 66% (697) or 83% (Nalm6) of GC-active STARR-seq sites mapped to HGRs. Moreover, 697 and Nalm6 GC-responsive open chromatin sites were significantly enriched at GC-active STARR-seq sites compared to accessible chromatin sites that were not GC-responsive (Fisher’s Exact p < 5.4 × 1010; Fig. 3C, Supplementary Fig. 21). Notably, GC-responsive chromatin sites exhibiting enhanced accessibility were significantly more enriched at GC-active STARR-seq sites compared to sites showing reduced accessibility (Chi-square p < 2.2 × 10−16; 2.2–3.5-fold).
We further compared active STARR-seq sites before (0 h) and after (GC) GC treatment to identify sites that are active only after GC exposure (i.e., GC-specific active STARR-seq site). Shared active STARR-seq sites between 0 h and GC conditions exhibited substantially greater overlap with promoter locations compared to 0 h-specific and GC-specific active STARR-seq sites, which had a higher proportion of promoter-distal elements (Supplementary Fig. 22). Importantly, GREs were significantly more enriched at GC-specific active STARR-seq sites in 697 cells compared to 0 h-specific active STARR-seq sites (Chi-square p = 1.3 × 10−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 STARR-seq sites were also more closely associated with GC-responsive upregulated genes compared to 0 h-specific active STARR-seq sites (24 h DEGs, fold change >2; K-S test p = 0.05 [697] and p = 0.002 [Nalm6]). Overall, ATAC-STARR-seq validated gene regulatory activity for 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 outcome [4,5,6, 10, 49, 50]. 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 XVI [17, 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 cell lines 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 upregulated gene that has also been previously correlated with prognostic features in ALL patients [51] and functions with Groucho as a canonical Wnt signaling repressor [52, 53]. In line with the recognized consensus GRE, the alternative T allele, which disrupts the GRE, is associated with greater GC resistance in primary cells (Supplementary Fig. 23). 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. 4B, Supplementary Fig. 24). CRISPR/Cas9 genome editing was further used to delete the TLE1 HGR in Nalm6 cell lines and resulted in greater GC resistance (Fig. 4C, Supplementary Fig. 25).
Notably, HiChIP identified long-range promoter interactions between rs7045812 and the TLE1 promoter (Supplementary Fig. 26). In line with this observation, TLE1 expression was significantly reduced in TLE1 HGR deleted cells compared to wild-type cells, thereby validating the gene regulatory effects of this HGR on TLE1 expression (Supplementary Fig. 27). We further tested the impact of TLE1 disruption on GC resistance using CRISPR/Cas9 knockout in Nalm6 cell lines (Supplementary Fig. 28). GC drug response assays using prednisolone uncovered that TLE1 disruption led to greater GC resistance (Supplementary Fig. 29), and this was concordant with both lower TLE1 expression being associated with GC resistance in primary cells (Supplementary Fig. 30) and with greater GC resistance in TLE1 HGR deleted cells (Fig. 4C). 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 laboratory [26], ex vivo GC drug sensitivity assays were performed on 19 of these primary ALL cells. We stratified cells by GC resistance (Supplementary Fig. 31) and identified 1929 sites where differential chromatin accessibility was associated with GC resistance (FDR < 0.1) through inter-subtype (all samples) and/or intra-subtype (only ETV6-RUNX1 or hyperdiploid samples) analysis (Fig. 5A, 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 (Supplementary Fig. 32). 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.2 × 10−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 (Supplementary Fig. 33). 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 resistance [18] and CRISPR interference (CRISPRi) screening (Supplementary Fig. 34).
CRISPRi using Enhancer-i [54] 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. 5B and Supplementary Table 3). Following a 72 h drug selection with prednisolone, we identified numerous sgRNAs exhibiting significant enrichment (697 = 1844, Nalm6 = 774; FDR < 0.05; Supplementary Tables 4, 5), whereas control sgRNAs overall did not show strong or preferential enrichment (Fig. 5C). By further aggregating all sgRNA data at each GR site using MAGeCK [34] (~6.5 sgRNAs per site), we ranked GR sites by log2-transformed fold change enrichment (Supplementary Tables 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. 35). 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 GREAT [46] and examined if these genes had been implicated in GC resistance [18]. This integrative analysis identified 35 top HGRs associated with 26 GC-responsive genes implicated in GC resistance (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 cell lines 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 (Supplementary Fig. 36). CRISPR/Cas9-mediated deletion of this HGR significantly reduced TLE1 expression (Fig. 6B, Supplementary Fig. 37) and led to greater GC resistance (Fig. 6C), supporting the functional role of this distal HGR in TLE1 gene regulation and GC resistance.
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 a neighboring HGR (Peak43) and (ii) an HGR (Peak42) with STARR-seq activity in 697 cell lines (Fig. 7A). ROR1 is a receptor tyrosine kinase-like orphan receptor for Wnt5a that can induce activation of noncanonical Wnt signaling [55] and has been associated with survival of TCF3-PBX ALL [56]. However, ROR1 is also expressed in other ALL molecular subtypes (Supplementary Fig. 38). All three HGRs exhibited greater chromatin occlusion in GC-resistant samples and these alterations were further supported by greater DNA methylation and lower ROR1 expression in GC-resistant samples (Supplementary Figs. 39, 40). In line with these data, lower ROR1 expression in primary cells was associated with greater GC resistance in an independent patient cohort (Supplementary Fig. 41) and CRISPR/Cas9 disruption of ROR1 in 697 cell lines led to greater GC resistance (Supplementary Figs. 42, 43). Because ROR1 is repressed by GCs, GC activity at this gene locus is anti-apoptotic (RNA-seq log2 fold change = −0.61, Supplementary Fig. 44). GC-induced ROR1 repression is also consistent with decreased STARR-seq activity at Peak42 at 24 h and a significant reduction in H3K27ac enrichment at Peak43 after GC treatment (FDR < 0.01). Concordant with these observations, individual CRISPR/Cas9-mediated deletions of the two distal HGRs ablated GC-induced ROR1 repression and decreased GC resistance (Fig. 7B, C, Supplementary Fig. 45), but baseline ROR1 expression was not significantly impacted (Supplementary Fig. 46). Overall, these data uncovered that cis-regulatory epigenetic disruptions to GC responses are a mechanism impacting GC resistance in ALL.
Discussion
In this study we map the genome-wide response to GCs in ALL cell lines and identified pervasive effects on the chromatin landscape, including the identification of GC-responsive chromatin sites and SEs, with most harboring GR occupancy. Transcriptomic and three-dimensional looping information further established a role for GC-responsive chromatin sites and HGRs in gene regulation, whereas investigations of TFs involved in the GC response uncovered a repression of AP-1 genes and AP-1 bound cis-regulatory elements implicated in cell proliferation and anti-apoptotic processes. We identified a direct role for GR in the repression of AP-1 TFs through HiChIP looping between HGRs and AP-1 promoters, and in the repression of AP-1 cis-regulatory elements through direct GR occupancy. To functionally evaluate cis-regulatory activities, STARR-seq was employed and validated thousands of GC-responsive chromatin sites, most of which mapped to HGRs.
To determine the functional impact of cis-regulatory disruptions to GC responses in drug resistance we integrated our GC response maps with genomic data from patient samples. We identified genetic and/or epigenetic disruptions to HGRs at the TLE1 and ROR1 gene loci that impact GC resistance. Importantly, because TLE1 and ROR1 are involved in canonical and noncanonical Wnt signaling [52, 53, 55], respectively, these analyses suggest that cis-regulatory disruptions to Wnt signaling is a mechanism impacting GC resistance in childhood ALL. However, the effects of GCs on TLE1 are pro-apoptotic, with greater GC-mediated TLE1 expression associated with lower GC resistance, whereas the effects of GCs on ROR1 are anti-apoptotic. ROR1 is repressed by GCs, and ROR1 disruption leads to greater GC resistance. Consistently, deletion of the distal ROR1 HGRs that are utilized for GC-induced repression promotes increased GC sensitivity. Our data further suggest that both distal HGRs appear to be necessary for robust ROR1 repression, but not for maintaining baseline ROR1 expression, which may be regulated by additional redundant cis-regulatory elements. Collectively, these results suggest that the gene regulatory activities of GCs are not exclusively pro-apoptotic, as has been previously reported [19], and alterations to GC-mediated anti-apoptotic processes also appear to play a role in resistance. Future directions should center on defining molecular mechanisms that link Wnt signaling to GC resistance, as well as follow-up investigations of the top epigenetically altered HGRs and variants we identified in this study (Supplementary Tables 1, 8). In addition, the identification of both genetic and epigenetic alterations at HGRs at the TLE1 gene locus as well as a correlation between TLE1 expression and prognostic features in ALL patients [51] all highlight TLE1 as an important gene for further follow-up studies.
Collectively, our study mapped gene regulatory responses to GCs in ALL cells using orthogonal functional genomic assays and high-throughput reporter assays. Our data further suggest that genetic and epigenetic disruptions to this gene regulatory response impact GC resistance in primary cells from patients.
Data availability
All functional genomic data from cell lines have been deposited into the Gene Expression Omnibus (GSE175484).
References
Moriyama T, Relling MV, Yang JJ. Inherited genetic variation in childhood acute lymphoblastic leukemia. Blood. 2015;125:3988–95.
Cooper SL, Brown PA. Treatment of pediatric acute lymphoblastic leukemia. Pediatr Clin North Am. 2015;62:61–73.
Tasian SK, Loh ML, Hunger SP. Childhood acute lymphoblastic leukemia: Integrating genomics into therapy. Cancer. 2015;121:3577–90.
Den Boer ML, Harms DO, Pieters R, Kazemier KM, Gobel U, Korholz D, et al. Patient stratification based on prednisolone-vincristine-asparaginase resistance profiles in children with acute lymphoblastic leukemia. J Clin Oncol. 2003;21:3262–8.
Kaspers GJ, Veerman AJ, Pieters R, Van Zantwijk CH, Smets LA, Van Wering ER, et al. In vitro cellular drug resistance and prognosis in newly diagnosed childhood acute lymphoblastic leukemia. Blood. 1997;90:2723–9.
Pieters R, Huismans DR, Loonen AH, Hahlen K, van der Does-van den Berg A, van Wering ER, et al. Relation of cellular drug resistance to long-term clinical outcome in childhood acute lymphoblastic leukaemia. Lancet. 1991;338:399–403.
Dordelmann M, Reiter A, Borkhardt A, Ludwig WD, Gotz N, Viehmann S, et al. Prednisone response is the strongest predictor of treatment outcome in infant acute lymphoblastic leukemia. Blood. 1999;94:1209–17.
Hunger SP, Mullighan CG. Acute lymphoblastic leukemia in children. N Engl J Med. 2015;373:1541–52.
Pui CH. Genomic and pharmacogenetic studies of childhood acute lymphoblastic leukemia. Front Med. 2015;9:1–9.
Holleman A, Cheok MH, den Boer ML, Yang W, Veerman AJ, Kazemier KM, et al. Gene-expression patterns in drug-resistant acute lymphoblastic leukemia cells and response to treatment. N Engl J Med. 2004;351:533–42.
Bhojwani D, Pui CH. Relapsed childhood acute lymphoblastic leukaemia. Lancet Oncol. 2013;14:e205–17.
Weikum ER, Knuesel MT, Ortlund EA, Yamamoto KR. Glucocorticoid receptor control of transcription: precision and plasticity via allostery. Nat Rev Mol Cell Biol. 2017;18:159–74.
Inaba H, Pui CH. Glucocorticoid use in acute lymphoblastic leukaemia. Lancet Oncol. 2010;11:1096–106.
Buentke E, Nordstrom A, Lin H, Bjorklund AC, Laane E, Harada M, et al. Glucocorticoid-induced cell death is mediated through reduced glucose metabolism in lymphoid leukemia cells. Blood Cancer J. 2011;1:e31.
Dyczynski M, Vesterlund M, Bjorklund AC, Zachariadis V, Janssen J, Gallart-Ayala H, et al. Metabolic reprogramming of acute lymphoblastic leukemia cells in response to glucocorticoid treatment. Cell Death Dis. 2018;9:846.
Kruth KA, Fang M, Shelton DN, Abu-Halawa O, Mahling R, Yang H, et al. Suppression of B-cell development genes is key to glucocorticoid efficacy in treatment of acute lymphoblastic leukemia. Blood. 2017;129:3000–8.
Paugh SW, Bonten EJ, Savic D, Ramsey LB, Thierfelder WE, Gurung P, et al. NALP3 inflammasome upregulation and CASP1 cleavage of the glucocorticoid receptor cause glucocorticoid resistance in leukemia cells. Nat Genet. 2015;47:607–14.
Autry RJ, Paugh SW, Carter R, Shi L, Liu J, Ferguson DC, et al. Integrative genomic analyses reveal mechanisms of glucocorticoid resistance in acute lymphoblastic leukemia. Nat Cancer. 2020;1:329–44.
Poulard C, Kim HN, Fang M, Kruth K, Gagnieux C, Gerke DS, et al. Relapse-associated AURKB blunts the glucocorticoid sensitivity of B cell acute lymphoblastic leukemia. Proc Natl Acad Sci USA. 2019;116:3052–61.
Wandler AM, Huang BJ, Craig JW, Hayes K, Yan H, Meyer LK, et al. Loss of glucocorticoid receptor expression mediates in vivo dexamethasone resistance in T-cell acute lymphoblastic leukemia. Leukemia. 2020;34:2025–37.
van Galen JC, Kuiper RP, van Emst L, Levers M, Tijchon E, Scheijen B, et al. BTG1 regulates glucocorticoid receptor autoinduction in acute lymphoblastic leukemia. Blood. 2010;115:4810–9.
Mullighan CG, Zhang J, Kasper LH, Lerach S, Payne-Turner D, Phillips LA, et al. CREBBP mutations in relapsed acute lymphoblastic leukaemia. Nature. 2011;471:235–9.
Jones CL, Bhatla T, Blum R, Wang J, Paugh SW, Wen X, et al. Loss of TBL1XR1 disrupts glucocorticoid receptor recruitment to chromatin and results in glucocorticoid resistance in a B-lymphoblastic leukemia model. J Biol Chem. 2014;289:20502–15.
Pottier N, Yang W, Assem M, Panetta JC, Pei D, Paugh SW, et al. The SWI/SNF chromatin-remodeling complex and glucocorticoid resistance in acute lymphoblastic leukemia. J Natl Cancer Inst. 2008;100:1792–803.
Jing D, Huang Y, Liu X, Sia KCS, Zhang JC, Tai X, et al. Lymphocyte-specific chromatin accessibility pre-determines glucocorticoid resistance in acute lymphoblastic leukemia. Cancer Cell. 2018;34:906–21.e8.
Diedrich JD, Dong Q, Ferguson DC, Bergeron BP, Autry RJ, Qian M, et al. Profiling chromatin accessibility in pediatric acute lymphoblastic leukemia identifies subtype-specific chromatin landscapes and gene regulatory networks. Leukemia. 2021;35:3078–91.
Jeha S, Pei D, Choi J, Cheng C, Sandlund JT, Coustan-Smith E, et al. Improved CNS control of childhood acute lymphoblastic leukemia without cranial irradiation: St Jude total therapy study 16. J Clin Oncol. 2019;37:3377–91.
Corces MR, Buenrostro JD, Wu B, Greenside PG, Chan SM, Koenig JL, et al. Lineage-specific and single-cell chromatin accessibility charts human hematopoiesis and leukemia evolution. Nat Genet. 2016;48:1193–203.
Savic D, Ramaker RC, Roberts BS, Dean EC, Burwell TC, Meadows SK, et al. Distinct gene regulatory programs define the inhibitory effects of liver X receptors and PPARG on cancer cell proliferation. Genome Med. 2016;8:74.
Ramaker RC, Savic D, Hardigan AA, Newberry K, Cooper GM, Myers RM, et al. A genome-wide interactome of DNA-associated proteins in the human liver. Genome Res. 2017;27:1950–60.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.
Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9:R137.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Li W, Xu H, Xiao T, Cong L, Love MI, Zhang F, et al. MAGeCK enables robust identification of essential genes from genome-scale CRISPR/Cas9 knockout screens. Genome Biol. 2014;15:554.
Shlyueva D, Stampfel G, Stark A. Transcriptional enhancers: from properties to genome-wide predictions. Nat Rev Genet. 2014;15:272–86.
Creyghton MP, Cheng AW, Welstead GG, Kooistra T, Carey BW, Steine EJ, et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc Natl Acad Sci USA. 2010;107:21931–6.
Heintzman ND, Hon GC, Hawkins RD, Kheradpour P, Stark A, Harp LF, et al. Histone modifications at human enhancers reflect global cell-type-specific gene expression. Nature. 2009;459:108–12.
McDowell IC, Barrera A, D’Ippolito AM, Vockley CM, Hong LK, Leichter SM, et al. Glucocorticoid receptor recruits to enhancers and drives activation by motif-directed binding. Genome Res. 2018;28:1272–84.
Gertz J, Savic D, Varley KE, Partridge EC, Safi A, Jain P, et al. Distinct properties of cell-type-specific and shared transcription factor binding sites. Mol Cell. 2013;52:25–36.
Whyte WA, Orlando DA, Hnisz D, Abraham BJ, Lin CY, Kagey MH, et al. Master transcription factors and mediator establish super-enhancers at key cell identity genes. Cell. 2013;153:307–19.
Loven J, Hoke HA, Lin CY, Lau A, Orlando DA, Vakoc CR, et al. Selective inhibition of tumor oncogenes by disruption of super-enhancers. Cell. 2013;153:320–34.
Schmidt S, Rainer J, Riml S, Ploner C, Jesacher S, Achmuller C, et al. Identification of glucocorticoid-response genes in children with acute lymphoblastic leukemia. Blood. 2006;107:2061–9.
Mumbach MR, Rubin AJ, Flynn RA, Dai C, Khavari PA, Greenleaf WJ, et al. HiChIP: efficient and sensitive analysis of protein-directed genome architecture. Nat Methods. 2016;13:919–22.
Bentsen M, Goymann P, Schultheis H, Klee K, Petrova A, Wiegandt R, et al. ATAC-seq footprinting unravels kinetics of transcription factor binding during zygotic genome activation. Nat Commun. 2020;11:4267.
Duren Z, Chen X, Jiang R, Wang Y, Wong WH. Modeling gene regulation from paired expression and chromatin accessibility data. Proc Natl Acad Sci USA. 2017;114:E4914–23.
McLean CY, Bristor D, Hiller M, Clarke SL, Schaar BT, Lowe CB, et al. GREAT improves functional interpretation of cis-regulatory regions. Nat Biotechnol. 2010;28:495–501.
Wang X, He L, Goggin SM, Saadat A, Wang L, Sinnott-Armstrong N, et al. High-resolution genome-wide functional dissection of transcriptional regulatory regions and nucleotides in human. Nat Commun. 2018;9:5380.
Arnold CD, Gerlach D, Stelzer C, Boryn LM, Rath M, Stark A. Genome-wide quantitative enhancer activity maps identified by STARR-seq. Science. 2013;339:1074–7.
Bosanquet AG. Correlations between therapeutic response of leukaemias and in-vitro drug-sensitivity assay. Lancet. 1991;337:711–4.
Hongo T, Yajima S, Sakurai M, Horikoshi Y, Hanada R. In vitro drug sensitivity testing can predict induction failure and early relapse of childhood acute lymphoblastic leukemia. Blood. 1997;89:2959–65.
Brassesco MS, Pezuk JA, Cortez MA, Bezerra Salomao K, Scrideli CA, Tone LG. TLE1 as an indicator of adverse prognosis in pediatric acute lymphoblastic leukemia. Leuk Res. 2018;74:42–6.
Chodaparambil JV, Pate KT, Hepler MR, Tsai BP, Muthurajan UM, Luger K, et al. Molecular functions of the TLE tetramerization domain in Wnt target gene repression. EMBO J. 2014;33:719–31.
Daniels DL, Weis WI. Beta-catenin directly displaces Groucho/TLE repressors from Tcf/Lef in Wnt-mediated transcription activation. Nat Struct Mol Biol. 2005;12:364–71.
Carleton JB, Berrett KC, Gertz J. Multiplex enhancer interference reveals collaborative control of gene regulation by estrogen receptor alpha-bound enhancers. Cell Syst. 2017;5:333–44.e5.
Karvonen H, Perttila R, Niininen W, Hautanen V, Barker H, Murumagi A, et al. Wnt5a and ROR1 activate non-canonical Wnt signaling via RhoA in TCF3-PBX1 acute lymphoblastic leukemia and highlight new treatment strategies via Bcl-2 co-targeting. Oncogene. 2019;38:3288–300.
Bicocca VT, Chang BH, Masouleh BK, Muschen M, Loriaux MM, Druker BJ, et al. Crosstalk between ROR1 and the Pre-B cell receptor promotes survival of t(1;19) acute lymphoblastic leukemia. Cancer Cell. 2012;22:656–67.
Acknowledgements
We would like to thank the St. Jude Hartwell Center for next-generation sequencing and the St. Jude Center for Advanced Genome Engineering for CRISPR/Cas9 genome editing. We also thank Jeremy Hunt, Dylan Maxwell, Brandon Smart, and Monique Payton for technical support. This work was supported by the National Institutes of Health (R01CA234490, P30CA021765, P50GM115279) and the American Lebanese Syrian Associated Charities (ALSAC). The content of this paper is solely the responsibility of the authors and does not necessarily represent the official views of the U.S. National Institutes of Health.
Author information
Authors and Affiliations
Contributions
Conceptualization, DS; Methodology, BPB, JDD, SMP, CL, DS; investigation, BPB, JDD, BSH; formal analysis, BPB, JDD, YZ, KRB, QD, DCF, RJA, WJ, CS, YF, DS; patient sample collection and sample acquisition, KRC, CHP, MVR, JJY, WEE; writing—original draft, BPB, DS; writing—review and editing, BPB, JDD, YZ, KRB, QD, DCF, RJA, WJ, BSH, CS, KRC, YF, CHP, SMP, MVR, JJY, CL, WEE, DS; funding acquisition, WEE, MVR, DS.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Bergeron, B.P., Diedrich, J.D., Zhang, Y. et al. Epigenomic profiling of glucocorticoid responses identifies cis-regulatory disruptions impacting steroid resistance in childhood acute lymphoblastic leukemia. Leukemia 36, 2374–2383 (2022). https://doi.org/10.1038/s41375-022-01685-z
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1038/s41375-022-01685-z