Charting epimutation dynamics in human hematopoietic differentiation

Abstract

DNA methylation plays a critical role in hematopoietic differentiation. Epimutation is a stochastic variation in DNA methylation that induces epigenetic heterogeneity. However, the effects of epimutations on normal hematopoiesis and hematopoietic diseases remain unclear. In this study, we developed a Julia package called EpiMut that enabled rapid and accurate quantification of epimutations. EpiMut was used to evaluate and provide an epimutation landscape in steady-state hematopoietic differentiation involving 13 types of blood cells ranging from hematopoietic stem/progenitor cells to mature cells. We showed that substantial genomic regions exhibited epigenetic variations rather than significant differences in DNA methylation levels between the myeloid and lymphoid lineages. Stepwise dynamics of epimutations were observed during the differentiation of each lineage. Importantly, we found that epimutation significantly enriched signals associated with lineage differentiation. Furthermore, epimutations in hematopoietic stem cells (HSCs) derived from various sources and acute myeloid leukemia were related to the function of HSCs and malignant cell disorders. Taken together, our study comprehensively documented an epimutation map and uncovered its important roles in human hematopoiesis, thereby offering insights into hematopoietic regulation.

Keywords: DNA methylation; Epimutation; Hematopoietic differentiation; Hematopoietic stem cell.

Figure 1.

The development and evaluation of EpiMut. (A) The workflow of EpiMut for epimutation rate measurement. (B) Comparison of reads used in PDR analysis between EpiMut and WSH. (C) Comparison of covered CpG sites in PDR analysis between EpiMut and WSH. (D) Comparison of average PDR values calculated by EpiMut and WSH. (E) The running time of PDR calculation at different sequencing depths in EpiMut and WSH. *P < .05; ***P < .001 (unpaired 2-tailed Student t test). ns = not significant, PDR = proportion of discordant reads.

Figure 2.

Overview of epimutation among human hematopoietic cell types. (A) Boxplots display the average PDR level of 13 hematopoietic cell types, targeting 1000 cells within each library. The gray dashed line represents the mean PDR level of all cell types. (B) Violin plots showing the PDR values across all cell types for different genomic regions. (C–D) Heatmap illustrating absolute PDR level across genomic regions in (C) lymphoid and (D) myeloid cells. (E–F) The average PDR levels of (E) lymphoid cells and (F) myeloid cells around the TSS (±2 kb) of all RefSeq genes. CGI = CpG islands, CLP = common lymphoid progenitors, CMP = common myeloid progenitors, CTCF = CCCTC-binding factor, GMP = granulocyte/monocyte progenitors, HSC = hematopoietic stem cell, LINE = long interspersed nuclear element, LTR = long terminal repeat, MEP = megakaryocyte/erythroid progenitor, MLP = multi-lymphoid progenitor, MPP = multipotent progenitor, NK = natural killer, PDR = proportion of discordant reads, SINE = short interspersed nuclear element, TF = transcription factor, TSS = transcription start site.

Figure 3.

Epimutation distinguishes myeloid and lymphoid progenitors. (A) Heatmap showing the row-scaled PDR of DDRs defined in CMP and CLP. (B) GO enrichment analysis of DDR-related genes between CLP and CMP. The x-axis represents the negative logarithm of the P values. (C) Heatmaps showing the row-scaled PDR values within the promoters of stemness, lymphoid-specific, and myeloid-specific genes across hematopoietic progenitor cells. (D) TF motif enrichment of DDRs between CLP and CMP. CLP = common lymphoid progenitor, CMP = common myeloid progenitor, DDRs = differentially discordant regions, GO = gene ontology, HSC = hematopoietic stem cell, MPP = multipotent progenitor, PDR = proportion of discordant read, TF = transcription factor.

Figure 4.

Different epimutation between progenitors and terminally differentiated cells. (A) Illustration depicting the numbers of DDRs in pairwise comparisons between GMP and neutrophil/monocyte. (B–C) GO enrichment analysis of genes related to DDRs with lower PDR values in (B) neutrophils and (C) monocytes compared to GMPs. The x-axis represents the negative logarithm of the P values. (D) Boxplots showing PDR values of cell type-specific genes across GMP, neutrophil, and monocyte. *P < .05, **P < .01 (unpaired 2-tailed Student t test). (E) Illustration depicting the numbers of DDRs in pairwise comparisons between MLP and mature cells (NK, B, CD4⁺, and CD8⁺ T cells). (F) GO enrichment analysis of genes with lower PDR values in lymphoid mature cells (NK, B, CD4⁺, and CD8⁺ T cells) than MLP, respectively. The size of the dots represents the enrichment score, while the color indicates the negative logarithm of the P values. (G) TF motif enrichment of DDRs with higher PDR values in mature cells than progenitors. Both the colors and sizes of the dots indicate the motif enrichment scores. DDRs = differentially discordant regions, GMP = granulocyte/monocyte progenitor, GO = gene ontology, MLP = multi-lymphoid progenitor, NK = natural killer, PDR = proportion of discordant read, TF = transcription factor.

Figure 5.

HSCs exhibit differential epimutation across 4 tissues. (A) Boxplots display the average PDR levels of HSCs from different tissues. The gray dashed line indicates the mean PDR level of all libraries. *P < .05; **P < .01; ***P < .001 (unpaired 2-tailed Student t test). (B) Boxplots display the PDR values across different tissues for different genomic regions. (C) Region set enrichment analysis for tissue-specific-lower DDRs (top) and tissue-specific-higher DDRs (bottom) of HSCs. The y-axis represents the negative logarithm of the P values, which is calculated by LOLA. The horizontal dashed line corresponds to a significance threshold of .05 for the P values. (D) Heatmap showing the row-scaled PDR values of tissue-specific lowest DDRs. The right column shows the GO enrichment of DDR-related genes. (E) Boxplots compare the PDR scores from each selected HALLMARK gene set across different tissues of HSCs. PDR scores are defined as the mean PDR levels of the promoters within the corresponding gene set. *P < .05, **P < .01, ***P < .001 (unpaired 2-tailed Student t test). (F) Heatmaps showing the scaled PDR values within the promoters of HSC maintenance and stemness markers across different tissues of HSCs. (G) TF motif enrichment of tissue-specific DDRs with highest PDR. The color indicates the motif enrichment scores. BM = bone marrow, CB = cord blood, CGI = CpG islands, CTCF = CCCTC-binding factor, DDRs = differentially discordant regions, FL = fetal liver, GO = gene ontology, HSC = hematopoietic stem cell, LINE = long interspersed nuclear element, LOLA = locus overlap analysis, LTR = long terminal repeat, ns = not significant, PB = peripheral blood, PDR = proportion of discordant read, SINE = short interspersed nuclear element, TF = transcription factor.

Figure 6.

Epimutation associated with AML. (A) Boxplots display the average PDR levels of normal HSC (n = 3) and AML (n = 3). P values are shown above the box (unpaired 2-tailed Student t test). (B) Boxplots showing the PDR values for different genomic regions across normal HSC and AML. *P < .05, **P < .01, ***P < .001 (unpaired 2-tailed Student t test). (C) GO enrichment analysis of DDR-related genes between normal HSC and AML. The x-axis represents the negative logarithm of the P values. (D) Boxplots display the PDR levels within the promoters of ARHGAP22 and PARD6A across normal HSC and AML (left). Survival analysis was performed based on the expressions of ARHGAP22 and PARD6A in the TCGA-LAML patient cohort (right). *P < .05; **P < .01 (unpaired 2-tailed Student t test). AML = acute myeloid leukemia, CGI = CpG islands, CTCF = CCCTC-binding factor, DDR = differentially discordant region, GO = gene ontology, HSC = hematopoietic stem cell, LINE = long interspersed nuclear element, LTR = long terminal repeat, PDR = proportion of discordant read, SINE = short interspersed nuclear element, TF = transcription factor.