Super-Enhancer-Mediated RNA Processing Revealed by Integrative MicroRNA Network Analysis

Abstract

Super-enhancers are an emerging subclass of regulatory regions controlling cell identity and disease genes. However, their biological function and impact on miRNA networks are unclear. Here, we report that super-enhancers drive the biogenesis of master miRNAs crucial for cell identity by enhancing both transcription and Drosha/DGCR8-mediated primary miRNA (pri-miRNA) processing. Super-enhancers, together with broad H3K4me3 domains, shape a tissue-specific and evolutionarily conserved atlas of miRNA expression and function. CRISPR/Cas9 genomics revealed that super-enhancer constituents act cooperatively and facilitate Drosha/DGCR8 recruitment and pri-miRNA processing to boost cell-specific miRNA production. The BET-bromodomain inhibitor JQ1 preferentially inhibits super-enhancer-directed cotranscriptional pri-miRNA processing. Furthermore, super-enhancers are characterized by pervasive interaction with DGCR8/Drosha and DGCR8/Drosha-regulated mRNA stability control, suggesting unique RNA regulation at super-enhancers. Finally, super-enhancers mark multiple miRNAs associated with cancer hallmarks. This study presents principles underlying miRNA biology in health and disease and an unrecognized higher-order property of super-enhancers in RNA processing beyond transcription.

Keywords: Brd4; DGCR8; Drosha; broad H3K4me3 domain; cancer; microRNA; super-enhancer.

Figure 1. Super-enhancers mark cell-type-specific and abundant miRNAs

(A) Assignment of super-enhancers (SEs) and typical enhancers (TEs) to the proximal genes and miRNA genes. The distances from enhancers to the closest genes (X axis) and the closest miRNA genes (Y axis) are compared. Super-enhancers (SEs) and typical enhancers (TEs) spatially closely associated with miRNA genes (i) and with non-miRNA genes (ii) represent dots near and far from the red diagonal, respectively. (B) Contribution of enhancer-miRNA assignment to miRNA expression status in mESCs. The linkage score (b/a) along X axis, determined by two parameters (a, b) in Figure 1A, represents a relative relationship between distances from enhancers to the closest miRNA genes and the closest genes, as an indicator of relative closeness between SE/TEs and the nearest miRNA genes (See STAR methods). For each miRNA gene, only the closest enhancer is considered, and expression levels of the closest miRNAs from each enhancer in TT-FHAgo2 ESC cells (Y axis, right) are plotted versus the linkage score. Red lines indicate cumulative distribution curve of total miRNA expression (Y axis, left). Dots in the pink and yellow shades indicate SE- and TE-associated miRNAs, respectively. (C) Box plots of miRNA expression and Ago2-bound miRNA reads from SE-associated, TE-associated, and other miRNAs. (D) Plots of expression levels and tissue-type specificity scores (See STAR methods) for SE-miRNAs (red), TE-miRNAs (yellow), and other miRNAs (black). (E) ChIP-seq profiles for master transcription factors, OSN in ESCs, PU.1 in Pro-B cells, and MyoD in myotubes, around miR-290-295, miR-148a, and miR-1/133a2. (F) Venn diagrams of SE-miRNAs and TE-miRNAs in ESCs (blue border), Pro-B cells (green border), and myotubes (orange border). See also Table S1 and Figure S1.

Figure 2. SE-miRNAs in many human cell types and target avoidance phenomenon

(A) Heatmap showing the classification of SE-miRNAs across 26 human cell and tissue types. Each row is a miRNA gene and red color indicates the miRNA being associated with SEs in the respective cell type. (B) Specificity of SE-miRNAs and TE-miRNAs. Y and X axis indicate number of miRNA genes and number of observed cell types when a miRNA is associated with a super-enhancer (SE-miRNA) and a typical enhancer (TE-miRNA). (C) Examples of representative SE-miRNAs observed in diverse tissue types. miR-15a/16-1 are associated with SEs in multiple hematopoietic cell types and other cell types (*). (D) Gene Ontology terms for target genes of SE-miRNAs with corresponding p values. (E) Cross-correlation analysis between SE-miRNA association and miRNA target density map. Pearson correlation coefficients between the SE-miRNA association matrix and target density map including negative log10 p value for depletion of miRNA targets assessed by Kolmogorov–Smirnov test (Farh et al., 2005) are shown (See STAR methods). See also Table S2 and Figures S2 and S3.

Figure 3. Network-based visualization of super-enhancer composition

(A) Quantitative contribution of expression of miRNAs associated with OSN-Med1 SEs and H3K27Ac SEs to the total miRNA pool in V6.5 mESCs. (B) Analysis strategy for classification of SE constituents. (C) Correlation networks of ChIP-seq profiles of 64 factors for 2401 enhancer constituents included in OSN-Med1 and H3K27Ac SEs, determined by coregulation analysis. (D) Heatmap showing background-subtracted and normalized ChIP-seq density for enhancer constituent clusters. Each row represents a enhancer constituent. (E) Pie charts showing distribution of numbers of enhancer constituents (C1–6) included in each OSN-Med1 or H3K27Ac SE. (F) Gene track of Nanog locus. (G) Distribution of H3K4me3 breadth in ESCs, with a subset of exceptionally broad H3K4me3 peaks. (H) Genomic interval overlap between SE constituents (C1, C2) and broad H3K4me3 peaks. (I, J) Overlap (I) and expression levels (J) of miRNAs associated with OSN-Med1 SEs, H3K27Ac SEs, and broad H3K4me3 peaks. (K) Relationship between Mediator/H3K27Ac-defined SEs and broad H3K4me3 domains. See also Tables S3 and S4 and Figure S4.

Figure 4. CRISPR/Cas9 functional genomics of miRNA super-enhancers

(A) Top: ChIP-seq binding profiles of Oct4, Sox2, Nanog (OSN), and Mediator (Med1) at the miR-290-295 locus in mESCs. Bottom: Scheme of miR-290-295 pathway. (B) RT-qPCR showing expression levels of mature miR-290-295 and Lin28a in TT-FHAgo2 mESCs deleted of SE constituent with/without Ago2 induction by doxycycline (Dox). (C) Top: ChIP-seq binding profiles of MyoD at the miR-1a-1/133a-2 locus in myotubes. Bottom: Scheme of miR-1/133 pathway. (D) RT-qPCR showing expression levels of mature miR-1a and myogenic differentiation markers in C2C12 cells deleted of SE constituent with/without differentiation medium (DM). (E) Top: ChIP-seq binding profiles of PU.1 and Mediator (Med1) at the miR-148a locus in Pro-B cells. Bottom: Scheme of miR-148a pathway. (F) RT-qPCR showing expression levels of mature miR-148a, Bach2, and Blimp1 in Pro-B cells deleted of SE constituent with/without stimulation by imatinib (IM). (G–I) Comparison of pri-miRNA, pre-miRNA, and mature miRNA expression levels in Ago2-induced TT-FHAgo2 cells (G), differentiated C2C12 myotubes (H), and imatinib-stimulated Pro-B cells (I) with SE constituent deletion. Data are represented as mean + SD from two biological replicate experiments. See also Table S5 and Figure S5.

Figure 5. Association between super-enhancers and DGCR8/Drosha

(A) ChIP-qPCR analysis showing recruitment of DGCR8 and Drosha to pri-miRNA loci in Ago2-induced TT-FHAgo2 mESCs, differentiated C2C12 myotubes, and imatinib-stimulated Pro-B cells with deletion of SE constituents. % input values were normalized with pri-miRNA expression levels. Data are represented as mean + SD. (B) Distribution of p values of DGCR8 and Drosha ChIP peaks mapping to various genomic regions. Hairpin indicates genomic loci corresponding to pre-miRNAs. Insets show fraction of SEs and TEs containing at least one DGCR8 or Drosha peak. (C) ChIP-seq profiles of DGCR8 and Drosha at the miR-290-295 locus in mESCs. (D) Metagene profiles of ChIP density of DGCR8 and Drosha across TEs, SEs, and gene bodies in mESCs. (E) Hierarchical clustering of ChIP-seq binding profiles of Drosha/DGCR8 and other factors in SE constituents defined in Figures 3B–D. Pearson correlation coefficients between each factor are shown. (F) Identification of Drosha-enhanced endonucleolytic mRNA cleavage events in mESCs. Cumulative frequency of Deg-seq read abundance at Deg-seq peaks in wild type (WT) and Drosha KO mESCs is shown. (G) Impact of Drosha-enhanced mRNA cleavage events on the response of abundance of chromatin-associated RNA (ChrRNA) for DGCR8 depletion (left) and mRNA half-life (right) of SE-associated genes in mESCs. See also Table S6 and Figure S6.

Figure 6. Inhibition of super-enhancer-associated miRNA processing by JQ1

(A) Boxplot showing changes in ChIP-seq density of DGCR8 and Drosha in TSS and hairpin (pre-miRNA) regions of SE-miRNAs, enhancers, and promoters after JQ1 treatment (500 nM, 12 hr). (B) Boxplot showing changes in ChIP-seq density of DGCR8/Drosha in the promoter-proximal region and elongating gene body region and DGCR8/Drosha traveling ratio for transcriptionally active and DGCR8/Drosha-bound genes after JQ1 treatment (500 nM, 12 hr). Left panel shows schema of the calculation of traveling ratio. (C) Gene track of rRNA-depleted ChrRNA sequencing and pre-miRNA sequencing at miR-290-295 locus (JQ1: 500 nM, 12 hr). (D) Effects of JQ1 (500 nM) on ChrRNA abundance (12 hr) at gene levels in V6.5 mESCs. (E) Effects of JQ1 (500 nM) on chromatin-associated pri-miRNA expression (12 hr) in V6.5 mESCs and pre-miRNA expression (6 and 12 hr) in Dicer (Dcr) knockout mESCs. (F) Interdependency of SE constituents for DGCR8 and Drosha recruitment. ChIP-qPCR results for SE constituents are shown. Data are represented as mean + SD. (G) Model of SE-regulated pri-miRNA processing.

Figure 7. Super-enhancer-associated miRNAs in cancer

(A) Comparison of distribution of SEs around miRNA genes in normal tissues and in cancer cells. Each row is a miRNA gene. In the left part, colored box indicates miRNAs associated with SEs in the respective samples, and color indicates distinct tissue types (blue: neuronal, red: hematopoietic, green: digestive tract, orange: pancreas, purple: breast, light purple: lung, gray: others). In the right part, red and blue boxes indicate gain and loss of SEs in the respective tissue types, respectively. (B) Pie charts showing functional assignment of miRNAs with SE loss or gain into tumor-suppressive or oncogenic groups. (C) Cancer hallmark-related Gene Ontology terms enriched in target genes of miRNAs with SE loss or gain, with corresponding p values. (D) Association between SE-miRNAs and cancer hallmark traits adapted from (Hanahan and Weinberg, 2011). (E, F) Survival impacts of SE-miRNAs. The distribution of −log10 (p value) (E) and Cox coefficient (F) in a single miRNA-based univariate Cox proportional hazard model are shown for individual miRNAs. * P < 0.05 with Wilcoxon signed rank sum test. (G) Compound SE-miRNA signatures. Subsets of miRNAs with an individual prognostic value of 0.05 or 0.1 are selected in select cancer cells and used as a stratifier for Kaplan–Meier survival analysis. P values were calculated with two-tailed log-rank test. See also Tables S2 and S7 and Figure S7.