Nested epistasis enhancer networks for robust genome regulation

Abstract

Mammalian genomes have multiple enhancers spanning an ultralong distance (>megabases) to modulate important genes, but it is unclear how these enhancers coordinate to achieve this task. We combine multiplexed CRISPRi screening with machine learning to define quantitative enhancer-enhancer interactions. We find that the ultralong distance enhancer network has a nested multilayer architecture that confers functional robustness of gene expression. Experimental characterization reveals that enhancer epistasis is maintained by three-dimensional chromosomal interactions and BRD4 condensation. Machine learning prediction of synergistic enhancers provides an effective strategy to identify noncoding variant pairs associated with pathogenic genes in diseases beyond genome-wide association studies analysis. Our work unveils nested epistasis enhancer networks, which can better explain enhancer functions within cells and in diseases.

Fig. 1.. High-resolution multiplexed CRISPRi perturbation of ultralong-distance enhancers at the MYC locus reveals a nested two-layer epistasis network.

A, Top: The MYC locus regulated by multiple enhancers distributed over an ultralong distance (~1.9Mb). Bottom: Diagram showing the multiplexed CRISPRi screening for high-resolution dissection of enhancer interactions. K562 cells expressing the doxycycline (Dox)-inducible dCas9-KRAB are transduced by a pooled sgRNA library targeting single or double MYC enhancers. Cells are harvested to sequence the pairwise sgRNA enrichment before and after 30 doublings. sgT, targeting sgRNA; sgC, control sgRNA. B, A quantitative epistasis map of sgRNA pairs targeting all enhancer combinations in the MYC locus. Each dot represents the epistasis interaction score of a pair of sgRNAs smoothed by adjacent sgRNAs. C, A quantitative enhancer epistasis map at the MYC locus. D, A nested two-layer model for the enhancer epistasis network. E,F, qRT-PCR of MYC mRNA expression for perturbing SREs e3&e7 or e4&e7 (E), or non-SREs e1&e4 or e5&e7 (F). P=0.02, 1.13E-05, 0.13, 0.61, for e3&e7, e4&e7, e1&e4, e5&e7, respectively. Data are represented as individual biological replicates (dots) and the mean value (black bar). The purple area indicates the expected additive effect by plotting mean ± one standard derivation. P values are calculated by t-test.

Fig. 2.. A machine learning model for analyzing determinants of the SRE synergy.

A, An elastic net regularized linear regression model for predicting epistasis interaction scores. We selected features including the chromatin spatial interaction (SI_ij) and co-occupancy (CO_ij,k) of 38 TFs and 8 HM profiles. B, The relative importance of each feature group for predicting epistasis interaction scores. The representative feature has the highest correlation in that group (fig. S8A). m.s.e., mean squared error. C-F, Correlation between epistasis interaction scores and Z-scores normalized spatial contact (C) and BRD4 co-occupancy (E). D, Heatmap of normalized HiChIP interaction intensity between enhancers. F, Correlation between predicted SRE scores and observed epistasis interaction scores. In C, D, F: red, SREs; blue, non-SREs. The Pearson correlation coefficient (R) and P value are shown.

Fig. 3.. Experimental validation of predicted SREs at other genomic loci in different cell types.

A-D, Prediction and validation of SREs at BCL9 (A) and KTN1 loci (B) in K562 cells, and COX6C (C) and FOXP1 loci (D) in Jurkat cells. Top: Diagram showing multiple enhancers spanning an ultralong distance at each genomic locus. Bottom left: Rank of predicted SREs using the model. Dashed line represents the empirical threshold from the MYC locus. Orange dots indicate the validated SREs. Bottom right: qRT-PCR of mRNA expression for each gene when perturbing the predicted SREs. Data are represented as individual biological replicates (dots) and the mean value (black bar). The purple area indicates the expected additive effect by plotting mean ± one standard derivation. P values are calculated by t-test.

Fig. 4.. Perturbation of SREs leads to synergistic reduction of spatial contacts and BRD4 condensation at the genomic locus.

A-B, Spatial contacts between the promoter and enhancers measured by Trac-looping for the MYC locus upon perturbation of e3, e7, and e3&e7. Colors represent the log₂ fold change of spatial contacts normalized to the wildtype cells. Black boxes in (B) indicate synergistically decreased (more than additive) spatial contacts of e3&e7 pair perturbation. C-F, DNA-FISH colocalization between BRD4 and the MYC locus of representative K562 cells for 2D (C-D) and 3D image analysis (E-F) upon perturbation of e3, e7, and e3&e7. In (C&E), red, BRD4 immunofluorescence (IF) staining; green, DNA-FISH at the MYC locus; blue dashed line, nuclear periphery determined by DAPI staining (not shown); scale bars, 5μm. The rightmost column in (C) shows insets in the yellow boxes. Scale bars, 500 nm. Quantification of BRD4 and the MYC locus colocalization are shown for 2D (D) and 3D image analysis (F). In (D), percentage of loci with colocalization is shown on the top and percentage of cells (≥ 2 colocalization loci) is shown on the bottom; data is represented as mean ± standard error of the mean. In (F), each dot represents an individual locus. n = total loci, N = total cells. **** P < 0.0001 in Fisher’s exact test (D) or t-test (F) versus the expected additive effect (dashed line). G, A model to explain the synergy between SREs.

Fig. 5.. Synergistic interactions between predicted SRE variants influence gene expression and disease risk in an epistatic manner.

A-B, Analysis of predicted SRE variants at the MYC locus in K562 cells for influence on gene expression. A, quantile-quantile (QQ) plot showing the distribution of P values for the epistasis influence on MYC expression between e4&e6 variants (red) in LAML patients, compared to random permutations (grey); P value in Kolmogorov–Smirnov (KS) test. B, MYC expression in LAML patients stratified by e4&e6 SRE variants. * P < 0.05 in Wilcoxon test. C-G, Analysis of predicted SRE variants at the MYC locus in GM12878 cells for influence on gene expression and associated disease risk. C, Diagram showing the rank of predicted SREs; orange dots show top SREs. D, QQ plot showing the distribution of P values for the epistasis influence of Be1&Be7 variants (red) on MYC expression in the B lymphoblasts of 373 European individuals, compared to random permutations (grey). P value in KS test. E, MYC expression in the B lymphoblasts from individuals stratified by Be1&Be7 variant. ** P < 0.01 in Wilcoxon test. F-G, Calculated odds ratio on the relapse risk in acute lymphoblastic leukemia (ALL) (F) and Crohn's disease (CD) (G). Odds ratios are calculated by considering the genotypes of individual variants or both SRE variants. Colors represent the odds ratios. H-I, Analysis of predicted SRE variants at the CHD7 locus in GM12878 cells for influence on ALL. H, Diagram showing the rank of predicted SREs; orange dots show top SREs. I, Calculated odds ratio on the relapse risk in ALL. Odds ratios are calculated by considering the genotypes of individual variants or both SRE variants. Colors represent the odds ratios.

Fig. 6.. Genome-wide analysis of epistatic influence of SRE variants on disease risk.

A, Percentage of enhancer pairs with observed epistatic effects on ALL relapse risk for predicted SREs and non-SREs. B-C, Percentage of enhancer pairs (B) and genes (C) exhibiting interactive effects on ALL relapse risk. SRE pairs: enhancer pairs with top 40% SRE predicted score; non-SRE pairs: enhancer pairs with bottom 10% SRE predicted score. D, Comparison of identified ALL pathogenic genes between the SRE model and the traditional locus-by-locus model at different odds ratio levels. In all figures, *: P < 0.05; **: P < 0.01; ***: P < 0.001; ****: P < 0.0001 in Fisher’s exact test.