Pre-configuring chromatin architecture with histone modifications guides hematopoietic stem cell formation in mouse embryos

Abstract

The gene activity underlying cell differentiation is regulated by a diverse set of transcription factors (TFs), histone modifications, chromatin structures and more. Although definitive hematopoietic stem cells (HSCs) are known to emerge via endothelial-to-hematopoietic transition (EHT), how the multi-layered epigenome is sequentially unfolded in a small portion of endothelial cells (ECs) transitioning into the hematopoietic fate remains elusive. With optimized low-input itChIP-seq and Hi-C assays, we performed multi-omics dissection of the HSC ontogeny trajectory across early arterial ECs (eAECs), hemogenic endothelial cells (HECs), pre-HSCs and long-term HSCs (LT-HSCs) in mouse embryos. Interestingly, HSC regulatory regions are already pre-configurated with active histone modifications as early as eAECs, preceding chromatin looping dynamics within topologically associating domains. Chromatin looping structures between enhancers and promoters only become gradually strengthened over time. Notably, RUNX1, a master TF for hematopoiesis, enriched at half of these loops is observed early from eAECs through pre-HSCs but its enrichment further increases in HSCs. RUNX1 and co-TFs together constitute a central, progressively intensified enhancer-promoter interactions. Thus, our study provides a framework to decipher how temporal epigenomic configurations fulfill cell lineage specification during development.

Fig. 1. Annotating the large-scale 3D genome re-organization with chromatin states during HSC formation.

a Schematic of experimental design for multi-omics profiling across four populations. b The sampling scheme. Surface markers were provided in the methods. c Pie charts for the ratio of A/B compartment flips. B-to-A flips indicate a resulting compartment B to A change during differentiation, and A-to-B flips likewise indicate A to B; Transient flips indicate the change involving both A to B and B to A along the differentiation path. d Aggregate curves showing features of different histone modification ChIP signals around TAD boundary regions. The values were calculated with signals on ChIP-seq peaks. Each region was ±500 kb around TAD boundary centers, with 10 kb bin size. e Heatmap showing the clustering of top 1000 variable TAD boundaries (40 kb bin) by boundary strengths (row-scaled), and the corresponding H3K4me3 signals at TAD boundaries. Each row was one TAD boundary, identified by insulation scores. The minus value of the insulation score was boundary strength. The numbers of genes within each TAD boundary were shown. The ANOVA P-value (one-side) of H3K4me3 signals at TAD boundaries with increased strength among four stages was 0.922. The ANOVA P-value (one-side) of H3K4me3 signals at TAD boundaries with decreased strength among four stages was 0.353. f Gene ontology analysis for biological processes of TAD boundary clusters. The GO analysis was done by GREAT. P-value was calculated by two-sided binomial test and adjusted by Benjamini–Hochberg (BH) correction (FDR). The X-axis was −log10(FDR) of each term. g Violin plots of gene expression as in (f). Gene expression values were calculated and normalized with log2(TPM/10 + 1) from previous scRNA-seq data (GSE139389 and GSE67120). Each dot represented one single cell. In boundaries with increased strength, n = 166 genes. In boundaries with decreased strength, n = 429 genes. Box-and-whiskers plots represented the maxima, 75th percentile, median, 25th percentile, and minima. Pairwise comparisons between any two populations were performed through two-sample Wilcoxon Rank Sum test (two-side). *P-value < 0.05; **P-value < 0.01; ***P-value < 0.001; ****P-value < 0.0001; “ns” not significant. eAEC early AEC, pre pre-HSC, LT LT-HSC.

Fig. 2. Mapping developmental dynamics of intra-TAD connectivity and chromatin configuration of feature histone modifications.

a Heatmap showing clusters of the top 1000 variable TADs by normalized domain scores. b Heatmap of H3K27ac, H3K27me3, H3K4me1, and H3K4me3 signals in corresponding ChIP-seq peaks within TADs as sorted in (a). c Quantification curves for the asynchronous change of feature histone modifications within TADs and TAD domain scores in three main clusters. Mean values were calculated for TADs in each cluster. d Biological processes enriched in three main clusters of TADs. The GO analysis was done by GREAT. P-value was calculated by two-sided binomial test and adjusted by BH correction (FDR). e Violin plots showing an asynchronous change of TAD domain scores and histone modifications in C3 TADs. The log2(fold change) values between two adjacent stages were calculated. Box-and-whiskers plots represented the maxima, 75th percentile, median, 25th percentile, and minima. The two-sample Wilcoxon Rank Sum test (two-sided) was performed for C3 TADs (n = 599) between two adjacent stages. A early AEC, H HEC, p pre-HSC, L LT-HSC. f Accumulative curves showing the distribution of normalized paired-read counts in Hi-C at 10 kb matrices within different genomic distances. g Exemplification showing the increasing intra-TAD interactions with observed/expected Hi-C matrices (5 kb resolution) at Runx1 gene site, associated with feature histone modifications. h Curves quantifying the changing loop strength between the Runx1 P1 promoter and neighboring regions. The genomic coordinate is chr16:92,614,120–93,400,000, from −212 kb to +574 kb relative to Runx1 P1. The observed/expected Hi-C matrices at 5 kb resolution were used for calculation. The Runx1 P1 region in red was TSS ± 2.5 kb. The region 1 in yellow overlapped with Runx1 P2 promoter, while the region 2 in green represented a distal enhancer. *P-value < 0.05, **P-value < 0.01, ***P-value < 0.001, ****P-value < 0.0001, “ns” not significant. eAEC early AEC, pre pre-HSC, LT LT-HSC.

Fig. 3. Pre-configuration with feature histone modifications precedes chromatin loop dynamics.

a Heatmap showing clusters of top 1000 variable loops among enhancers and promoters by loop strengths (normalized CPM). Loops contained E–E, E–P, or P–P loops. b Heatmaps showing cognate H3K4me3, H3K27ac, and gene expression on anchors of each loop. Only reads on ChIP-seq peaks were used to calculate normalized signals of H3K4me3 and H3K27ac. The eAEC and HEC scRNA-seq data were from GSE139389. Pre-HSC and LT-HSC scRNA-seq data were from GSE67120 and GSE66954. Gene expression values were calculated and normalized with log2(TPM/10 + 1). c Percentages of genes in three loop clusters with down-regulated, transient, and up-regulated expression from eAECs to LT-HSCs. Transient genes meant those whose expression was up-regulated in some cases and down-regulated in others. d Biological process analysis of looping regions in three clusters. The inputs were these looping anchor regions for GO term analysis with GREAT. The −log10(FDR) values were shown after BH correction. e Track views showing the loop change trend among enhancers and promoters at the Ets1, Cebpg, and Meis1 loci. f Violin plots showing expression of representative genes in three clusters with log2(TPM/10 + 1). The scRNA-seq data contain 42 eAECs, 33 HECs, 50 pre-HSCs and 16 LT-HSCs. Box-and-whiskers plots represented the maxima, 75th percentile, median, 25th percentile, and minima. Wilcoxon Rank Sum test (two-sided) was performed between two adjacent stages. *P-value < 0.05, **P-value < 0.01, ***P-value < 0.001, ****P-value < 0.0001, “ns” not significant. eAEC early AEC, pre pre-HSC, LT LT-HSC.

Fig. 4. RUNX1 is engaged in looping between promoters and enhancers for priming genome interactions during HSC formation.

a Bubble plots showing gene expression and TF motif enrichment, identified at anchor regions of C3 loop clusters as in Fig. 3a, with most variable gene expression. TF motif analysis was performed by HOMER. P-value was calculated by two-sided binomial test. b Pie chart showing the portion of E–E, P–P, and E–P interactions with or without RUNX1 occupancy. c Heatmap showing hierarchical clustering based on the expression of genes with RUNX1-engaged E–P interactions. d Heatmap showing the corresponding RUNX1 binding and histone modifications on looping sites as in (b). Representative genes were listed on the right. The ANOVA P-values (one-sided) of TF or histone modification signals in C2 hematopoietic cluster were 1.27e−05 **** for RUNX1, 3.92e–12 **** for H3K27ac, 0.000216 *** for H3K27me3, and 0.0606 ns for H3K4me3. e Gene Ontology analysis of biological processes supporting the clustering. f Lines reflecting strengths of three sub-clusters of RUNX1-engaged E–P interactions among C2 genes. Those RUNX1-engaged E–P interactions existing from eAEC to LT-HSC were used for clustering. Smoothly fitted interactions dynamics into a function (method = ‘loess’) in each cluster, the black line indicates the predicted value. The thick white bars represented a 95% confidence interval. g Violin plots quantifying the strengths of three sub-clusters of RUNX1-engaged E–P interactions among C2 genes, corresponding to (f). For decreased interactions, n = 91. For ephemeral interactions, n = 60. For increased interactions, n = 113. Box-and-whiskers plots represented the maxima, 75th percentile, median, 25th percentile, and minima. The two-sample Wilcoxon Rank Sum test (two-sided) was performed between two adjacent stages. *P-value < 0.05; **P-value < 0.01; ***P-value < 0.001; ****P-value < 0.0001; “ns” not significant. eAEC early AEC, pre pre-HSC, LT LT-HSC.

Fig. 5. Identifying potential RUNX1 co-TFs in mediating enhancer–promoter interactions during HSC formation.

a Bar plot showing the portion of distal enhancers (with genomic distance to TSS > 5 kb) and promoters (with genomic distance to TSS ≤ 5 kb) within RUNX1-engaged E–P interactions. b Pie chart showing the portion of E–P interactions with RUNX1 occupancy on promoters, enhancers, or both. Schematic on the left shows the three types of genomic locations of RUNX1 occupancy. c Hierarchical clustering based on gene expression of potential RUNX1-interacting TFs identified through de novo TF motif discovery on the E–P looping regions as above. TF genes used for clustering were genes in families whose protein binding motifs were top 11 significant in all candidates identified in Supplementary Fig. 7b. Gene expressions were calculated by log2(TPM/10 + 1) from scRNA-seq data. d Bubble plots showing significances of TF motif enrichments as in (c), together with gene expression. TF motif analysis was performed by HOMER. P-value was calculated by two-sided binomial test. e Working model of multi-omics dynamics during HSC formation. RUNX1 occupies target promoters or/and enhancers primed with H3K27ac in eAECs. Repressive H3K27me3 gradually diminishes accompanying later H3K27ac increase. Other hematopoiesis related TFs are recruited to form a TF module with RUNX1 to facilitate enhancer–promoter interactions driving the core HSC gene expression. eAEC early AEC, pre pre-HSC, LT LT-HSC.