Histone hyperacetylation disrupts core gene regulatory architecture in rhabdomyosarcoma

Abstract

Core regulatory transcription factors (CR TFs) orchestrate the placement of super-enhancers (SEs) to activate transcription of cell-identity specifying gene networks, and are critical in promoting cancer. Here, we define the core regulatory circuitry of rhabdomyosarcoma and identify critical CR TF dependencies. These CR TFs build SEs that have the highest levels of histone acetylation, yet paradoxically the same SEs also harbor the greatest amounts of histone deacetylases. We find that hyperacetylation selectively halts CR TF transcription. To investigate the architectural determinants of this phenotype, we used absolute quantification of architecture (AQuA) HiChIP, which revealed erosion of native SE contacts, and aberrant spreading of contacts that involved histone acetylation. Hyperacetylation removes RNA polymerase II (RNA Pol II) from core regulatory genetic elements, and eliminates RNA Pol II but not BRD4 phase condensates. This study identifies an SE-specific requirement for balancing histone modification states to maintain SE architecture and CR TF transcription.

Figure 1.. Core regulatory circuitry includes SOX8 and is critical for FP-RMS

a, Core regulatory circuitry identified by analysis of motif-networks in super-enhancer (SE) associated transcription factors (TFs). Heatmap shows predicted CR TFs found in FP-RMS cell lines and tumors (n = 9 independent samples), clustered and colored by degree of connectivity (scaled to 1 = maximum connectivity in the sample). Expression of CR TFs is plotted on the right (box plots show quartiles, with whiskers showing the 1.5 × inter-quartile range, with data distribution as violin plots). b, Core regulatory TF validation ChIP-seq in RH4 ChIP-seq (a FP-RMS cell line). Metagene plots at ATAC-seq peaks in SEs were divided into SOX8-bound (solid lines, n = 839 peaks) and SOX-unbound (dotted lines, n = 1,190 peaks). c, Prevalence and co-occurrence of CR TFs in SEs in RH4. d, Functional genetic screening to evaluate TF dependencies in FP-RMS cells, using pooled CRISPR-cas9 methodology. Each sgRNA (n = 8,837) is designed to recruit Cas9 to DNA sequence coding for the DNA-binding domain(s) of each human TF (n = 1433). e, Depletion of sgRNAs targeting CR TFs, compared to all non-CR TFs, over 6 passages of negative selection. Top panel shows a histogram of the number of CR TFs (red bar). Lower panel shows a CRISPR score for each TF as an average log₂ fold change of sgRNAs (n = 6 distinct sequences per TF) targeting DNA-binding domains. Right panel shows summary statistics in box (quartiles) and whisker (1.5 × inter-quartile range) plot, with P value comparing extent of depletion between CR TFs and all other TFs calculated with an unpaired, two-sided student’s t test with Welch’s correction.

Figure 2.. Genetic dissection of core regulatory network reveals a SOX8-mediated myogenic blockade in RMS

(a) Connectivity of CR TFs shown by ChIP-seq of PAX3-FOXO1, MYOD, MYOG and SOX8 at super-enhancers (SEs) that interact with the promoters of PAX3, MYOD1, MYOG and SOX8. Heatmaps depict interaction frequency between genomic elements, assayed by H3K27ac HiChIP (one representative of two experiments in RH4 cells). (b) Co-regulation of CR TFs evaulated in RH4 cells with CRISPR-cas9 targeted to the DNA binding domains of PAX3, MYOD, MYOG or SOX8. Experiments were repeated orthogonally with shRNA with similar results (Supplementary Fig. 2e). (c) Gene set enrichment analysis shows divergent transcriptional impact of individual CR TF disruption in RH4 cells (gene sets available in Supplemental Table 2). Similar results were obtained using shRNA (Supplementary Fig. 2f). P value is generated using the GSEA algorithm of the enrichment score relative to the null distribution calculated with 1,000 permutations. Each bubble represents one gene set analyzed against an individual RNA-seq experiment (sgRNA target v. non-targeting sgRNA) at the indicated time points. (d) Activation and enrichment of myoblast differentiation genes upon sgRNA mediated disruption of PAX3-FOXO1 or SOX8, in constrast to depletion of the myogenic program upon MYOD or MYOG knockout, and was also seen upon shRNA knockdown (Supplementary Fig. 2g). (e) Model of auto-regulatory feed-forward (PAX3-FOXO1, MYOD1, MYOG, MYCN and other CR TFs) and negative feedback (from SOX8) circuitry in FP-RMS.

Figure 3.. Core regulatory TFs require HDAC1/2/3 for their transcription

a, Genome browser view at the SOX8 locus, including 3 super-enhancers (red bars) shown by ChIP-seq to be heavily decorated by CR TFs, HDACs (nuclear isoforms 1, 2 and 3), histone acetylation, other co-activators and RNA Pol2. Each track is one representative of two independent ChIP-seq experiments, except for HDAC3, MYCN, YY1 and MED1 which were performed once, and PAX3-FOXO1 and H3K27ac ChIP-seq data which were performed more than 4 times across different RH4 cell passages. b, Enrichment of CR TFs in SEs associated with genes that have high amounts of HDAC bound to their promoters (top, dotted line indicate genome-wide average), genes divided into 5 categories of HDAC ChIP-seq density (middle, violin plots), which correlated positively with expression (bottom, boxplots with median and quartiles, whiskers showing the 1.5 × inter-quartile range). Bins were defined using a single representative HDAC2 ChIP-seq, and similar binding was seen with HDAC1 and HDAC3 (see Supplementary Fig. 3a–b). c, Changes in gene expression upon HDAC1/2/3 inhibition with Entinostat for 6 hours in RH4 cells. Violin plots are overlaid with box plots of median and quartiles, whiskers showing the 1.5 × inter-quartile range. P values were calculated with a two-tailed unpaired t test with Welch’s correction. d, Time-course RNA-seq followed by geneset enrichment reveals that Entinostat causes an initial increase in core regulatory transcription (1 hr), followed by a selective downregulation at 6 and 24 hours in RH4 cells. FDR was calculated by GSEA with 1,000 permutations for each timepoint, and were for 1 hour, q = 0.006; 6 hours, q = 0; 24 hours, q = 0. e, Chromatin run-on and sequencing (ChRO-seq) to measure nascent RNA shows transcriptional increase at CR TFs MYOD1 and MYOG in 10 minutes, sustained for the first hour, and decreased at 6 hours of Entinostat treatment (1 μM in RH4 cells). f, Nascent RNA changes over time course treatment with Entinostat, as meansured by ChRO-seq and summarized by box plots (median with quartiles, whiskers representing 1.5 × interquartile range in the data). A two-sided, unpaired t test with Welch’s correction was used to calculate P values. g, Single-cell RNA-seq analysis of cells treated with DMSO (n = 2,925 cells) or Entinostat for 1 hr (n = 3,805 cells) or 6 hrs (n = 3,240 cells). Top panel, shows increased number of cells expressing the CR TFs SOX8 and MYOD1 at 1 hour of Entinostat and a decreased number at 6 hours. Bottom left panel, single-cell RNA-seq data of CR TFs abundance per cell in RH4 cells. These experiments were performed once. Bottom right shows heuristic boxplots and describes some relevant advantages of single cell RNA-seq over bulk RNA-seq.

Figure 4.. AQuA-HiChIP shows disruption of super-enhancer architecture by hyperacetylation

a, Super-enhancer dynamics upon HDAC inhibition, revealing an acute increase in H3K27ac after 1 hour, followed by a spread beyond the endogenous SE boundary at 6 hours. ChIP-Rx with exogenous spike in is reported as Reference-normalized Reads Per Kilobase per Million mapped reads (RRPKM). Shading shows the SEM. b, Acetylation changes quantified on diverse lysine residues on histone subunits H3 (K36, K27), H2B (K5) and H4 (K16) by ChIP-Rx. Reference-normalized Reads Per Million mapped reads (RRPM) differences for diverse acetylation sites on histones are shown for Core Regulatory Domains after treatment with either DMSO or Entinostat for 6 hours (left), or for H3K27ac upon treatment with HDAC1/2 inhibitor Merck60, HDAC3 inhibitor LW3 or HDAC1/2/3 inhibitor Entinostat (right). All experiments were in RH4 cells treated for 6 hours with 1 μM of the indicated inhibitors. Box plots show the median and quartiles, with whiskers showing the 1.5 × inter-quartile range. c, Absolute Quantification of Architecture HiChIP (AQuA-HiChIP) identifies spreading of SE-mediated contacts within the insulated neighborhood of CR TF MYOD1. Contact map is shown at 5-kb resolution (5-kb by 5-kb contact squares), and is scaled to AQuA normalized contacts per million. d, Gained aberrant contacts and lost SE-to-SE contacts are visualized by AQuA-Virtual 4C at MYOD1 from viewpoint anchor SE₁ (top) and SE₂ (bottom). Virtual 4C is representative of two replicate biotin captures and library preparations, both with similar results, and agrees with non-virtual 4C experiments at MYOD1 SEs under the same treatment conditions. e, SE contact spreading as seen by AQuA-HiChIP Aggregate Peak Analysis (APA) plots of all SE-to-SE contacts within insulated neighborhoods (SE_n to SE_intra, top) or SEs nearby but outside insulated neighborhood CTCF boundaries (SE_n to SE_inter, bottom). Resolution is shown at 10-kb by 10-kb squares. f, Endogenous p300 recruitment to MYOD1 SE elements. Chemical Epigenetic Modifier-114 (CEM-114, bi-functional FKBP-binder and p300 bromodomain binder) enables dCas9-guided recruitment of p300 to MYOD1 SE epicenters (sgEpicenters) or SE boundaries (sgBoundaries) in RH4 cells. Triplicate values are one representative of two independent cell treatments each with two RT-qPCR replicates that gave similar results.

Figure 5.. SE clusters and phase condensates are disrupted by hyperacetylation

a, ChIP-Rx binding sites for CR TFs, ranked by change in binding upon 6 hours of Entinostat treatment in RH4 cells. b, Entinostat-induced changes in binding of CR TFs shown at cis-regulatory elements for SOX8, PAX3-FOXO1, MYOD1, and MYOG in RH4 cells. c, RNA Pol2 unloading along all genic positions of MYOD1 at 1 and 6 hours of Entinostat treatment, and associated changes in H3K27ac, as measured by ChIP with spike-in normalization (ChIP-Rx). d, Clusters of RNA Pol2 tagged with GFP (in live FP-RMS cells, RH4), imaged in a time course with or without addition of HDAC inhibitor Entinostat. Scale bar equals 5 μM. Representative single puncta of RNA Pol2 clusters over time are shown (2.8 μM squares) from control or HDAC inhibitor treated RH4 cells. Cross-section profiles of pixel brightness for corresponding time points are shown below image frames. Images are from wide-field (see Online Methods) and are representative of 50 images across 2 independent experiments with similar results. Alternative imaging modality iSIM gave similar results to widefield and is shown in Supplementary Figure 7f. e, Live cell imaging of BRD4 (endogenously tagged with mEGFP in RH4 cells) treated for 6 hours with DMSO (left) or Entinostat (right). Scale bar represents 5 μM. Images are representative of 20 images across 2 independent experiments with similar results. f, Quantification of RNA Pol2 (top) and BRD4 (bottom) binding changes upon HDAC1/2/3 inhibition with Entinostat (1 μM, RH4 cells) by ChIP-Rx signal. g, HDAC inhibition induced changes in binding of key factors within Core Regulatory Domains, measured by ChIP-Rx in RH4 cells, summarized by box plots representing the median and quartiles, whiskers showing 1.5 × interquartile range. h, A model whereby histone hyperacetylation disrupts normal looping interactions, reduces binding of certain CR TFs, reduces binding of RNA Pol2 and causes loss of transcription at super-enhancer-dependent CR TF genes.