High-Resolution Single-Cell DNA Methylation Measurements Reveal Epigenetically Distinct Hematopoietic Stem Cell Subpopulations

Abstract

Increasing evidence of functional and transcriptional heterogeneity in phenotypically similar cells examined individually has prompted interest in obtaining parallel methylome data. We describe the development and application of such a protocol to index-sorted murine and human hematopoietic cells that are highly enriched in their content of functionally defined stem cells. Utilizing an optimized single-cell bisulfite sequencing protocol, we obtained quantitative DNA methylation measurements of up to 5.7 million CpGs in single hematopoietic cells. In parallel, we developed an analytical strategy (PDclust) to define single-cell DNA methylation states through pairwise comparisons of single-CpG methylation measurements. PDclust revealed that a single-cell epigenetic state can be described by a small (<1%) stochastically sampled fraction of CpGs and that these states are reflective of cell identity and state. Using relationships revealed by PDclust, we derive near complete methylomes for epigenetically distinct subpopulations of hematopoietic cells enriched for functional stem cell content.

Keywords: DNA methylation; epigenetics; hematopoietic stem cell; single cell.

Figure 1

Overview of Experimental Methods (A) Schematic of the phenotypes studied. (B) Schematic of the PBAL method. Cells are lysed in 96-well plates to release genomic DNA (gDNA) and then subjected to bisulfite conversion that simultaneously converts and shears the gDNA into fragments. Random hexamers are then used to regenerate double-stranded DNA that is then end-repaired, a-tailed, and ligated with indexed adapters before low-cycle PCR amplification. Libraries are then pooled and sequenced on an Illumina sequencing platform. See also Figures S1 and S2.

Figure 2

Concordance Analysis of Neighboring Methylated CpGs by Their Genomic Separation (A) A schematic for how adjacency is calculated. For a randomly subsampled number of CpGs (CpG1), the concordance of methylation of 100 CpGs before and after CpG1 was recorded along with their separation. (B) Analysis of bulk versus single LSK cells across CpG sites genome wide. Curves represent the weighted average in 100-bp bins. Horizontal lines indicate the probability that two randomly sampled CpG sites have the same methylation state based on their genome-wide methylation (“baseline concordance”). Baseline concordance was calculated as the probability of sampling two methylated or two unmethylated CpGs, which is equal to the square of the average fractional methylation rate plus 1 minus the square of the average fractional methylation rate. (C) Analysis of bulk versus single cells across selected regions. For each panel, we considered only CpGs as CpG1 if they were found within each genomic feature.

Figure 3

Pairwise Analysis of Single Cells Reveals Subsets within the Murine ESLAM Phenotype (A) Schematic showing how the PD values are calculated between all paired comparisons of single cells. (B) ESLAM cells have a lower overall PD compared with LSK cells and all cell types analyzed. Every pairwise comparison between cells denoted on the x axis is summarized as a box plot with the distribution of PD values shown on the y axis. p values were calculated using a two-sided t test. (C) PDclust of CpGs associated with genes implicated in HSC function (Cabezas-Wallscheid et al., 2014) separate ESLAM and LSK cells, with some LSK cells exhibiting an ESLAM epigenetic signature. The rows and columns are symmetrical and represent single cells. The cells are shaded to represent the PD between each pair of cells, with red representing highly dissimilar and yellow representing highly similar. Meth, average genome-wide CpG methylation; cpg_count, number of distinct CpG sites recovered. (D) Same as (C) but instead considering all CpGs regardless of their genomic position. (E) Multidimensional scaling using PD calculated from (C) used directly as input. (F) MDS analysis of (D) reveals group 1 at the epicenter of single ESLAM cells with group 2 surrounding the central cluster. See also Figures S3–S5.

Figure 4

Pairwise Dissimilarity Applied to Existing Datasets Is Able to Distinguish Cells Based on Cell Type and Treatment. (A and B) Clustering (A) and MDS scaling of PD values (B) calculated from HL60 and K562 cells (Farlik et al., 2015) separates cells by cell type, and reveals patterns of treatment-induced differentiation. Cells are shaded to represent the pairwise dissimilarity between each pair of cells, with red representing highly dissimilar and yellow representing highly similar. Meth, average genome-wide CpG methylation; cpg_count, number of distinct CpG sites recovered. (C and D) Clustering (C) and MDS scaling (D) of embryonic stem cells (Smallwood et al., 2014, Angermueller et al., 2016) separates cells by the culture medium in which they were grown, and shows some cells in a transition state.

Figure 5

Pairwise Dissimilarities of Human Hematopoietic Cells with Different Phenotypes (A) PD of CD49f and other CD34⁺ subsets separate cells by phenotype with some overlaps. The rows and columns are symmetrical and represent single cells. Cells are shaded to represent the PD between each pair of cells, with red representing highly dissimilar and yellow representing highly similar. (B) MDS of PD values shows that CD49f cells cluster in the middle, with GMPs, CLPs, and MLPs branching out in separate directions. (C) Single cells have lower PD values compared with cells of the same type versus cells of a different phenotype. The distribution of PD values when a cell is compared with either a cell of the same or different phenotype is plotted as a box plot. See also Figure S6.

Figure 6

Donor Variation in CD49f Cell Methylomes (A) Single CD49f cells still cluster by donor after removing CpGs within 200 bp of non-C and -G SNV locations. (B) MDS projection of PD values onto 2D space remains largely unchanged despite taking into account SNVs. (C) PD values as a function of comparisons versus cells either from the same or the other donor. Boxes are colored if only CpGs near SNVs, outside SNVs, or all CpGs were considered.

Figure 7

Heterogeneity within the Human CD49f Cell Compartment (A) A rare subset of CD49f cells (group 1, 11% and 9% of all CD49f cells for donors 1 and 2, respectively) cluster away from the rest of the cells after projection of PD values with MDS. (B) Distributions of surface markers obtained during index sorting of single cells belonging to cluster 1 or cluster 2, split by donor. (C) Observed over expected enrichment of DMRs. The ratios were calculated by dividing the fraction of DMRs that overlap each region set by the fraction of the genome each region set occupies. (D) Gene enrichment of DMRs that are hypomethylated in each comparison. The comparisons between group 1 and group 2 were separate from the comparisons between CD49f cells and published data for other CD34⁺ phenotypes. Bars are gray if their FDR-corrected binomial q value is <0.1. (E) Example of a DMR near the SERPING1 gene. Methylation values were smoothened for each population of cells with Bsmooth and plotted. Tick marks on the x axis represent the location of CpG dinucleotides. Computationally defined DMRs are highlighted in blue. The genes track shows high confidence protein-coding transcripts obtained from Gencode v75. See also Figure S7.