Whole-genome fingerprint of the DNA methylome during chemically induced differentiation of the human AML cell line HL-60/S4

Abstract

Epigenomic regulation plays a vital role in cell differentiation. The leukemic HL-60/S4 [human myeloid leukemic cell line HL-60/S4 (ATCC CRL-3306)] promyelocytic cell can be easily differentiated from its undifferentiated promyelocyte state into neutrophil- and macrophage-like cell states. In this study, we present the underlying genome and epigenome architecture of HL-60/S4 through its differentiation. We performed whole-genome bisulphite sequencing of HL-60/S4 cells and their differentiated counterparts. With the support of karyotyping, we show that HL-60/S4 maintains a stable genome throughout differentiation. Analysis of differential Cytosine-phosphate-Guanine dinucleotide methylation reveals that most methylation changes occur in the macrophage-like state. Differential methylation of promoters was associated with immune-related terms. Key immune genes, CEBPA, GFI1, MAFB and GATA1 showed differential expression and methylation. However, we observed the strongest enrichment of methylation changes in enhancers and CTCF binding sites, implying that methylation plays a major role in large-scale transcriptional reprogramming and chromatin reorganisation during differentiation. Correlation of differential expression and distal methylation with support from chromatin capture experiments allowed us to identify putative proximal and long-range enhancers for a number of immune cell differentiation genes, including CEBPA and CCNF Integrating expression data, we present a model of HL-60/S4 differentiation in relation to the wider scope of myeloid differentiation.

Keywords: DNA methylation; Differentiation; Epigenomic regulation; HL60; Long range interactions; Promyelocyte.

Fig. 1.

Analysis of DNA methylome upon chemical induction of differentiation of HL-60/S4 cells. (A) Schematic diagram of the experimental design of the study. (B) Whole-genome CpG methylation rate density plot. The upper left density plot shows that all three cell states (UN, RA and TPA) have very similar genome-wide CpG methylation rates. The subsequent density plots show the CpG methylation rates for each cell state separately. (C) Box plots summarising the distribution of CpG methylation rates per sample replicates for the ∼22 million CpGs with coverage ≥10× in all samples. The upper and lower limits of the boxes represent the first and third quartiles, respectively, and the black horizontal line is the median. The whiskers indicate the variability outside the upper and lower quartiles. (D) Principal component analysis of the WGBS data for the three cell states. Principal component 1 and 2 separate TPA from UN and RA cells. (E) Circular representation of DNA methylation rates for the different treatments. CpG methylation rates (colour scale beige–blue) were averaged over 10-Mb windows and are presented as heatmap tracks. The heatmaps show the DNA methylation change (heatmap black–white-red) with respect to the samples in the adjacent tracks.

Fig. 2.

DMP analysis. (A) Number of DMPs identified with Fisher exact test for each comparison. In the bar plot, the x-axis labels indicate the comparisons, and the y-axis indicates the number of identified DMPs. The number of hypermethylated DMPs are shown by the red bars, and the hypomethylated by blue bars. (B) Enrichment of genomic features in the hypermethylated (left) and hypomethylated (right) DMPs in RA and TPA compared UN cells. Genes, exon and TSS features are of all genes in the gencode v19 gene models. PTSS denotes the TSS of protein coding genes. CpGI denotes CpG islands. Enhancers were identified from ENCODE chromHMM chromatin segmentation. CTCF binding sites were from ENCODE ChIP-seq experiments. Epichromatin and LADs have been identified to associate with the nuclear envelope. Integrations denote chromatin interacting regions as defined to chromatin conformation capture experiments in ENCODE. DNA, LINE, LTR, SINE, satellite and simple repeats are all classes of repeats. (C) Density plot of the methylation rates of RA DMPs in UN and RA. Hypermethylated DMPs increased from methylation rates of 0 to 0.2, and hypomethylated posited decreased from 0.2 to 0. Hyper and hypomethylated DMPs found in RA compared to UN are denoted by (RA hyper) in pink and (RA hypo) in grey respectively; the methylation rate of the DMPs in UN and RA are denoted as solid and dashed lines, respectively. (D) As panel C, but for TPA. TPA hypermethylated DMPs (TPA hyper, in red) and hypomethylated DMPs (TPA hypo, in blue). (E) Unsupervised cluster analysis reveals 12 DMP ‘modules’. 12 modules were identified, and a heatmap representation showed distinct methylation rates in each of the three cell states. Modules are demarcated by alternating grey and beige bars on the right. (F) Genomic feature enrichment in the 12 modules. Module M6 shows a strong enrichment for enhancers, but not TSS regions. Simple repeat elements are enriched in states M7, M9 and M12. The genomic features are explained in B.

Fig. 3.

Association between promoter methylation and gene expression. (A) Scatter plot of genes with more than 1.5× log2 fold change in expression against the DNA methylation different (minimum of 0.2) of upstream transcription factor binding region between RA and UN cells. The number of genes is identified in the title of the plot. In each quadrant the number of genes in reported, with the upper right quadrant and lower left indicating genes with possible positive correlation with DNA methylation, and the upper left and lower right quadrants indicating genes that have negative correlation with DNA methylation. (B) As panel A, but for TPA. The scatter plots shows a comparable number of genes with increased TSS methylation with up and downregulation, however genes with reduced TSS methylation also exhibit reduced expression. (C) The distribution of Pearson correlation coefficient values for each gene between the gene expression value and the average methylation of DMPs overlapping with TSS in UN, RA and TPA cells. Increase in DNA methylation is associated with nearly as many up and downregulated genes as with decrease in methylation. While we expect to see primarily negative correlation, the histogram shows comparable numbers of highly positive and highly negative correlations. (D) As panel C, but for TPA. This histogram shows an elevated number of genes with positive correlation between 0.5 and 0.8.

Fig. 4.

Chromatin interacting regions explain gene expression correlation with distal DMPs. (A) CEBPE shows a strong inverse correlation between its expression and methylation of DMPs in its downstream region. Tracks (from top to bottom): genomic coordinates; gene model of CEBPE; simplified chromatin segmentation based on E029 (primary monocyte cells from peripheral blood) from the ROADMAP epigenome project; chromatin interacting regions (K562 cells from the ENCODE project using IM-PET); bar plot depicts correlations of DMPs with CEBPE differential gene expression, indicating a number of downstream DMPS negatively correlating with CEBPE expression (where a bar below 0 indicates negative correlation and a bar above 0 indicates positive correlation between DNA methylation and gene expression); line plot of DNA methylation rate of CpGs in UN, RA and TPA cells; coverage plot of RNAseq expression in UN, RA and TPA cells. (B) The cyclin-F-box protein coding gene CCNF interacts with a distant upstream region, which regulated its expression through methylation. Panel description is similar to A, except that the regions between PGP and CCNF are cut, the CHIA-PET (from K562 cells from the ENCODE project) interaction between the upstream of PGP and CCNF are connected, the top DMP-gene expression correlation panel is of proximal CpGs with the distal gene (‘PGP DNAme∼CCNF Exp’ and ‘CCNF DNAme∼PGP Exp’) and the lower correlation panel is of proximal CpGs with the proximal gene (‘PGP DNAme∼PGP Exp’ and ‘CCNF DNAme∼CCNF Exp’). For CCNF, no proximal DMP correlated highly with CCNF expression (right track ‘CCNF DNAme∼CCNF Exp’), however a number of DMPs upstream of PGP (in a region interacting with the CCNF promoter in K562 cells) exhibited negative correlation with CCNF differential expression (right track ‘CCNF DNAme∼PGP Exp’). Together, this could imply that these DMPs upstream of PGP play a role in CCNF expression. Please note that the chromatin interaction data are derived from different cells than the chromatin state data (K562 and primary monocytes respectively). While both are derived from myeloid cells, they are different from our myeloid HL-60/S4 cells.

Fig. 5.

Chemical differentiation model of HL-60/S4 showing the transcription factors that may play an essential role in determining cell fate. Downregulation or upregulation of gene expression are denoted by ‘−’ or ‘+’ respectively. Genes with no sign attached show that their levels are maintained at similar levels as in UN (promyelocytic) state. Genes with associated DMPs either being increased or decreased in DNA methylation are coloured as red and blue, respectively.