DNA Methylation Dynamics of Human Hematopoietic Stem Cell Differentiation

Abstract

Hematopoietic stem cells give rise to all blood cells in a differentiation process that involves widespread epigenome remodeling. Here we present genome-wide reference maps of the associated DNA methylation dynamics. We used a meta-epigenomic approach that combines DNA methylation profiles across many small pools of cells and performed single-cell methylome sequencing to assess cell-to-cell heterogeneity. The resulting dataset identified characteristic differences between HSCs derived from fetal liver, cord blood, bone marrow, and peripheral blood. We also observed lineage-specific DNA methylation between myeloid and lymphoid progenitors, characterized immature multi-lymphoid progenitors, and detected progressive DNA methylation differences in maturing megakaryocytes. We linked these patterns to gene expression, histone modifications, and chromatin accessibility, and we used machine learning to derive a model of human hematopoietic differentiation directly from DNA methylation data. Our results contribute to a better understanding of human hematopoietic stem cell differentiation and provide a framework for studying blood-linked diseases.

Keywords: DNA methylation profiling; bioinformatic lineage reconstruction; cell type prediction; hematopoietic stem cell differentiation; immature lymphoid progenitors; lymphoid-myeloid lineage commitment; megakaryocyte maturation; reference epigenome mapping; single-cell sequencing; whole genome bisulfite sequencing.

Graphical abstract

Figure 1

Charting the DNA Methylation Landscape of Human Hematopoietic Differentiation (A) Conceptual outline of human hematopoietic differentiation, highlighting the 17 hematopoietic cell types whose genome-wide DNA methylation patterns were profiled in this study. Arrows denote established differentiation trajectories, dashed arrows indicate uncertainty about the in vivo differentiation potential of lymphoid progenitors, and the inset illustrates the sorting of four subsets of immature multi-lymphoid progenitors. (B) Fluorescence-activated cell sorting panel used to purify 10 stem and progenitor cell types from peripheral blood. (C) Violin plots and boxplots showing the distribution of DNA methylation levels in 5-kb tiling regions for hematopoietic cell types sorted from peripheral blood. (D) Distribution of DNA methylation levels across cell types for different sets of genomic regions. Gene and promoter annotations are based on GENCODE, CpG islands are from the UCSC Table Browser, enhancer elements are from Ensembl, and tiling regions were calculated with a custom script. (E) Distribution of average DNA methylation levels across cell types for putative regulatory regions annotated by the Ensembl BLUEPRINT Regulatory Build. (F) DNA methylation at putative regulatory regions for illustrative gene loci. Black bars denote the position of regions annotated by the BLUEPRINT Regulatory Build, and dashed horizontal black lines indicate sample medians for the respective regions. Colored vertical bars connect the highest and lowest DNA methylation levels that have been measured in any sample of the indicated cell type. dis, distal element; prox, proximal element; TSS, transcriptional start site. See also Figure S1 and http://blueprint-methylomes.computational-epigenetics.org.

Figure 2

Comparison of DNA Methylation Maps for HSCs and MPPs Isolated from Four Different Sources (A) HSCs and MPPs were sorted from the peripheral blood of three healthy donors (Figures 1A and 1B), and in addition from fetal liver (HSCs only), cord blood (HSCs and MPPs), and bone marrow (HSCs and MPPs). (B) Bar plots showing the numbers of differentially methylated regions in pairwise comparisons between HSCs and MPPs from different sources, based on the BLUEPRINT Regulatory Build regions (FDR-adjusted p ≤ 0.05, absolute difference ≥ 0.167 percentage points), calculated with RnBeads (Assenov et al., 2014). (C) Region set enrichment analysis for genomic regions with lower DNA methylation in peripheral blood-derived HSCs compared with bone-marrow-derived HSCs. Enrichment was determined using LOLA (Sheffield and Bock, 2016). Each dot represents one ChIP-seq dataset, and the horizontal dashed line corresponds to a significance threshold of 0.05 on the adjusted p-value calculated by LOLA using Fisher’s exact test. (D) Source-specific DNA methylation at the IKBKE gene locus. Reduced DNA methylation levels in peripheral blood-derived HSC at two putative regulatory regions (BLUEPRINT Regulatory Build) is associated with detectable expression of the IKBKE gene specifically in this cell population (bar plot). (E) Heatmap showing DNA methylation levels for regions that have lower DNA methylation in peripheral blood-derived HSCs than in bone-marrow-derived HSCs and also overlap with CTCF binding sites in the LOLA Core database (left). The second heatmap (right) shows the overlap of these regions with transcription factor binding sites that a LOLA analysis of this region set identified as enriched. Rows were arranged by hierarchical clustering with complete linkage based on the Euclidean distances between the DNA methylation profiles. See also Figure S2 and http://blueprint-methylomes.computational-epigenetics.org.

Figure 3

DNA Methylation Differences Associated with Myeloid-Lymphoid Lineage Commitment (A) Scatterplot showing average DNA methylation levels for myeloid progenitors (CMP, MEP, GMP) and lymphoid progenitors (CLP, MLP0, MLP1, MLP2, MLP3) across BLUEPRINT Regulatory Build regions (Pearson’s r = 0.97). Differentially methylated regions were identified with RnBeads (FDR-adjusted p ≤ 0.05, absolute difference ≥ 0.167 percentage points). (B) Average DNA methylation in each cell type relative to the mean over all samples aggregated over all regions with lower DNA methylation in lymphoid progenitors (top) and myeloid progenitors (bottom). Error bars correspond to the standard error. (C) Region set enrichment analysis for regions with lower DNA methylation in lymphoid progenitors (left) or myeloid progenitors (right). Enrichment was determined using LOLA. Colored dots represent ChIP-seq experiments for transcription factors in the indicated lineage. Dot size denotes the log-odds ratio, and the numbers in the legend (“X/Y”) refer to significantly enriched region sets (X) versus all analyzed region sets (Y). The horizontal dashed line represents the significance threshold (adjusted p ≤ 0.05). (D) Mean-adjusted DNA methylation relative to the average CpG methylation levels for each individual 10-, 50-, and 1,000-cell sample averaged across ChIP-seq peaks for GATA1 (left) and TAL1 (right). ^∗∗∗p ≤ 0.001 (two-tailed Wilcoxon test). (E) Two-dimensional projection of all 10-, 50-, and 1,000-cell samples from peripheral blood using principal component analysis based on the mean-adjusted DNA methylation across all transcription factor binding datasets identified by LOLA. The first two principal components are shown, and the numbers in parentheses indicate the percentage of variance explained. (F) Heatmap displaying mean-adjusted DNA methylation for all single-cell DNA methylation profiles across the same region sets as in (E). Rows and columns are arranged by hierarchical clustering with Euclidean distance and complete linkage. The labels on the right summarize the transcription factors and cell types for the major branches of the row dendrogram. (G) Distribution of HSC (top) and MPP (bottom) samples derived from fetal liver (FL), cord blood (CB), bone marrow (BM), and peripheral blood (PB) when projected onto the first principal component from (E). p.p., percentage points. See also Figure S3 and http://blueprint-methylomes.computational-epigenetics.org.

Figure 4

Characterization of MLP Populations by DNA Methylation and In Vitro Differentiation Assays (A) Sorting panel for purifying four MLP populations from peripheral blood. (B) Two-dimensional projection of all 10-, 50-, and 1,000-cell MLP samples using principal component analysis based on the mean-adjusted DNA methylation relative to the average CpG methylation levels across all region sets in the LOLA Core database. The first two principal components are shown, the numbers in parentheses indicate the percentage of variance explained, and the density plots (top and right) summarize the distribution of cell types along the two principal components. (C) Heatmap displaying the mean-adjusted DNA methylation for all MLP samples across the 2 × 25 genomic region sets that contributed most strongly to the first principal component (PC1, top) and the second principal component (PC2, bottom). Rows and columns are arranged by hierarchical clustering with Euclidean distance and complete linkage. The row labels indicate the cell type and ChIP-seq target of the corresponding LOLA region sets. (D) Differentiation potential of CLPs and four MLP populations measured by in vitro culture of single cells on MS-5 stroma with cytokines promoting lymphoid and myeloid differentiation. The percentage of colonies that show lymphoid (CD19⁺ or CD56⁺) as well as myeloid (CD14⁺ or CD15⁺) markers was determined by flow cytometry. (E) Differentiation potential of the MLP populations measured as the percentage of colonies containing B cells (CD19⁺), granulocytes (CD15⁺), or NK cells (CD56⁺) in flow cytometry. ^∗p ≤ 0.05 (Fisher’s exact test). See also Figure S4 and http://blueprint-methylomes.computational-epigenetics.org.

Figure 5

DNA Methylation Analysis of Megakaryocyte Maturation (A) Conceptual outline of megakaryocyte development from MEPs via a maturation phase involving endomitotic genome replication and a concomitant increase in ploidy. (B) Scatterplot comparing mean-adjusted DNA methylation (relative to the average CpG methylation level in each sample) for all region sets in the LOLA Core database between megakaryocyte at the 2N and at the 32N stage of ploidy. Region sets that were significantly less methylated in 32N (n = 14, bottom right) or in 2N megakaryocyte (n = 6, top left) are highlighted with filled circles (p ≤ 0.05, Wilcoxon test, absolute difference ≥ 10 p.p.). Point colors indicate the magnitude of difference, and the size is proportional to statistical significance [−log₁₀(p)]. (C) Mean-adjusted DNA methylation in region sets that gain (top) or lose DNA methylation (bottom) as identified in (B), plotted across ploidy stages of megakaryocyte maturation. Error bars correspond to the standard error. p.p., percentage points. See also Figure S5 and http://blueprint-methylomes.computational-epigenetics.org.

Figure 6

Integrative Analysis of Gene Expression, Histone Modifications, and Chromatin Accessibility (A) Heatmap showing row-normalized expression levels for 656 differentially expressed genes (FDR-adjusted p ≤ 0.05, |log₂FC| ≥ 1) between myeloid progenitors (CMP, GMP) and lymphoid progenitors (CLP, MLP0, MLP1, MLP2, MLP3) determined using DEseq2 (Love et al., 2014). Rows and columns were arranged by hierarchical clustering with Euclidean distance and complete linkage. (B) Top 10 most highly enriched Gene Ontology terms associated with genes overexpressed in lymphoid (top) and myeloid (bottom) progenitors based on the Enrichr software (Kuleshov et al., 2016). (C) Density scatterplot contrasting myeloid-lymphoid differences in DNA methylation at gene promoters with expression differences of the corresponding genes. Selected genes with strong differences in DNA methylation (absolute difference ≥ 10 p.p.) and gene expression (|log₂FC| ≥ 1) are highlighted. (D) Boxplots showing histone modification levels for open chromatin-associated H3K4me1, active enhancer-linked H3K27ac, and polycomb-associated H3K27me3 in regions that were differentially methylated between myeloid and lymphoid progenitors (Figure 3A). Histone modification levels were calculated from multiple ChIP-seq datasets for myeloid cells (neutrophils, monocytes, and macrophages, in red) and lymphoid cells (NK cells, B cells, and CD4⁺/CD8⁺ T cells, in blue). Brackets identify two-tailed Mann-Whitney U tests. ^∗p ≤ 0.05, ^∗∗∗p ≤ 0.001. (E) Heatmap showing DNA methylation levels (columns) in regions with cell-type-specific chromatin accessibility based on published ATAC-seq data for hematopoietic cell types (rows). Numbers in parentheses denote the number of chromatin accessible regions specific to each cell type. (F) Distribution of DNA methylation levels across regions with cell-type-specific chromatin accessibility. (G) Composite plots showing DNA methylation averages across regions with cell-type-specific chromatin accessibility. CpGs in the neighborhood of these regions were annotated with coordinates relative to their start and end (x axis). CpGs with a relative coordinate of 0 and 1 are located at the start and end of a region, respectively, and the coordinates −1 and 2 correspond to one region length upstream and downstream of the region. The curves show cubic spline smoothing of DNA methylation levels per cell type across accessible regions. p.p., percentage points; n.s., not significant. See also Figure S6 and http://blueprint-methylomes.computational-epigenetics.org.

Figure 7

Data-Driven Reconstruction of the Human Hematopoietic Lineage using Machine Learning (A) Conceptual outline of the machine learning approach used to predict cell type, to identify signature regions, and to infer cellular differentiation landscapes. (B) Confusion matrix showing the frequency of misclassification based on 10-fold cross-validation of cell-type classifiers trained and evaluated on 319 stem and progenitor samples (all 10-, 50-, and 1,000-cell pools) from peripheral blood. (C) Heatmap showing average DNA methylation levels of merged replicates (one column for each cell type in each donor) for the 1,234 signature regions extracted from a classifier trained on all peripheral blood-derived stem and progenitor samples. Regions (rows) were arranged using hierarchical clustering with Euclidean distance and complete linkage. (D) Region set enrichment analysis for the signature regions using LOLA. Colored dots represent ChIP-seq experiments for transcription factors in the indicated lineage. Dot size denotes the log-odds ratio, and the numbers in the legend (“X/Y”) refer to significantly enriched region sets (X) versus all analyzed region sets (Y). The vertical dashed line represents the significance threshold (adjusted p ≤ 0.05). (E) Two-dimensional projection of merged replicates (one point for each cell type in each donor) using principal component analysis based on average DNA methylation levels in the signature regions. The first two principal components are shown, and the numbers in parentheses indicate the percentage of variance explained. (F) Hematopoietic lineage reconstruction using the prediction propensities of DNA methylation-based classifiers as a measure of similarity between cell types. Nodes in the graph represent cell types, and edges are weighted by class probabilities of cross-prediction. An automated edge-weighted graph layout algorithm was used to define the positions of the nodes. See also Figure S7 and http://blueprint-methylomes.computational-epigenetics.org.