Advertisement

Analysis of the DNA methylome and transcriptome in granulopoiesis reveals timed changes and dynamic enhancer methylation

Michelle Rönnerblad, Robin Andersson, Tor Olofsson, Iyadh Douagi, Mohsen Karimi, Sören Lehmann, Ilka Hoof, Michiel de Hoon, Masayoshi Itoh, Sayaka Nagao-Sato, Hideya Kawaji, Timo Lassmann, Piero Carninci, Yoshihide Hayashizaki, Alistair R. R. Forrest, Albin Sandelin, Karl Ekwall, Erik Arner and Andreas Lennartsson for the FANTOM consortium

Key Points

  • In granulopoiesis, changes in DNA methylation preferably occur at points of lineage restriction in low CpG areas.

  • DNA methylation is dynamic in enhancer elements and appears to regulate the expression of key transcription factors and neutrophil genes.

Abstract

In development, epigenetic mechanisms such as DNA methylation have been suggested to provide a cellular memory to maintain multipotency but also stabilize cell fate decisions and direct lineage restriction. In this study, we set out to characterize changes in DNA methylation and gene expression during granulopoiesis using 4 distinct cell populations ranging from the oligopotent common myeloid progenitor stage to terminally differentiated neutrophils. We observed that differentially methylated sites (DMSs) generally show decreased methylation during granulopoiesis. Methylation appears to change at specific differentiation stages and overlap with changes in transcription and activity of key hematopoietic transcription factors. DMSs were preferentially located in areas distal to CpG islands and shores. Also, DMSs were overrepresented in enhancer elements and enriched in enhancers that become active during differentiation. Overall, this study depicts in detail the epigenetic and transcriptional changes that occur during granulopoiesis and supports the role of DNA methylation as a regulatory mechanism in blood cell differentiation.

Introduction

Hematopoiesis is the process by which blood cells of all lineages develop from a common hematopoietic stem cell (HSC). This process is controlled by timed expression of lineage-specific combinations of transcription factors. Accumulating evidence demonstrates that epigenetic changes play a central role in differentiation, providing cellular memory and stabilizing lineage choices. Cells undergo dramatic epigenetic changes during hematopoiesis, and epigenetic states define the plasticity of the developing blood cells.1-5 DNA methylation is involved in regulation of transcription and heterochromatin formation and contributes to the control of hematopoiesis.6-8 DNA is predominantly methylated at cytosines in CG dinucleotides (CpG). Historically, the main focus has been on methylation changes in genomic areas with high densities of CpGs, called CpG islands (CGIs). However, it was recently suggested that methylation changes outside of CGIs might be more relevant for cell differentiation.9,10

In both the lymphoid and myeloid lineages, genes that undergo methylation changes in hematopoiesis tend to encode proteins with functions specific for the mature cell type.7,8,11-13 DNA methylation changes are essential in making the choice between the lymphoid and the myeloid lineages in early hematopoiesis. Experiments in mice with mutated DNA methyltransferase 1 (DNMT1), an enzyme involved in maintenance of CpG methylation, demonstrated that it is required for B-cell development but not for myeloid development.6 Moreover, DNA methylation levels have been shown to increase on lymphoid commitment but decrease during myeloid commitment.8,13

Although there have been recent reports describing methylation changes in murine granulopoiesis,7,11,14 no previous studies have been made in the human system, with high cellular resolution ranging from multipotent progenitors to terminally differentiated granulocytes.

Altered DNA methylation is a general hallmark of cancer, including acute myeloid leukemia (AML). AML is caused by malignant transformation of myeloid progenitor cells and is the most common acute leukemia in adults. The DNA methylation machinery is often perturbed in AML including recurrent mutations in DNMT3A (DNA methyltransferse 3A).15-17 Also, DNA methylation changes have prognostic power in AML,18 and the DNMT inhibitor 5-azacitidine has been approved by the US Food and Drug Administration and European Union for treatment of myelodysplastic syndrome and subtypes of AML.19 Consequently, more knowledge of epigenetic mechanisms in myeloid development should provide useful insights into pathology and treatment of AML.

In this study, we investigated the DNA methylome and transcriptome of human cells in 4 separate differentiation stages in granulopoiesis, ranging from the multipotent common myeloid progenitor (CMP) to terminally differentiated bone marrow neutrophils (PMNs). To this end, we used a DNA methylation array with extensive genomic coverage, mRNA expression arrays, and cap analysis of gene expression (CAGE).

We report cell population–specific changes of DNA methylation. The main reduction of CpG methylation coincides with the loss of oligopotency associated with the transition from granulocyte-macrophage progenitor to promyelocyte/myelocyte (GMP-PMC).20 DNA methylation changes were mainly localized outside of CGIs and co-occurred with changes in gene expression and the motif activity,21 ie, transcriptional activity of target genes, of key hematopoietic transcription factors. Enhancer elements specifically activated during granulocytic maturation displayed enrichment of methylation changes, indicating a regulatory function of DNA methylation in these regions possibly relevant to differentiation control. Thus, by combining a set of genome-wide molecular techniques with a high cell stage resolution, we are able to present a detailed characterization of epigenetic and transcriptional changes during granulopoiesis. This work is part of the Functional Annotation Of The Mammalian Genome 5 (FANTOM5) project. Data downloads, genomic tools, and copublished manuscripts are summarized at http://fantom.gsc.riken.jp/5/.

Materials and methods

Isolation of cells

Bone marrow samples were aspirated from voluntary healthy donors, and cell populations representing different maturation stages of granulopoiesis were isolated using fluorescence-activated cell sorter (FACS) as previously described.22 The cell sorting and 4,6 diamidino-2-phenylindole staining of nuclei are described in supplemental Figure 1 and supplemental Materials on the Blood Web site. The study was approved by the ethics committee (permission numbers LU-281-00 and 1012/408-31/3). Informed consent was obtained in accordance with the Declaration of Helsinki.

RNA isolation and microarray analysis

RNA was extracted with Trizol from CMP, GMP, PMC, and PMN cells in biological replicates from 5 donors. Microarray analysis was performed using Illumina HumanHT-12 v3 and v4 arrays according to the manufacturer’s recommendations. Differentially expressed genes were selected based on P < .05 (repeated-measures analysis of variance with Benjamini-Hochberg adjustment for multiple hypothesis testing) and paired fold change ≥2 in any pairwise comparison between 2 cell populations. At least 1 of the two cell stages compared was required to have ≥3 replicates with detection P < .05. Gene ontology (GO) analysis was performed using DAVID Bioinformatics Resources.23 Array data are available at Gene Expression Omnibus (GEO) (accession number GSE50797). For more information, see supplemental Materials.

CAGE

RNA from 2 biological replicates of each cell stage were analyzed with CAGE,24 modified to the single molecule sequencing system.25 CAGE tags were normalized to library size as tags per million (TPM). The following FANTOM5 libraries passed quality control and were used in analysis of CAGE data: CMP (library ID CNhs12523); GMP (CNhs12524, CNhs12528); PMC (CNhs12525, CNhs12529); and PMN (CNhs12526, CNhs12530). CAGE tags were clustered using a single linkage method, where CAGE tags within 20 bp of each other on the same strand were merged into clusters. CAGE library preparation, sequencing, and genome mapping is described in detail in the supplement of Forrest et al.26

Identification, expression quantification, and additional analysis of enhancers from CAGE-derived bidirectionally transcribed loci were performed as in Andersson et al.27 Identification of putative enhancer-promoter associations were derived as described elsewhere,27 using a set of 88 myeloid FANTOM5 CAGE libraries. For details, see supplemental Materials. Motif activity response analysis (MARA) was performed as described in ref. 21. Briefly, assuming that expression of a promoter (CAGE cluster) is proportional to the number of binding sites for transcription factors in the proximal region and the activity of these transcription factors, for a given DNA binding motif, an activity measure (motif activity) is fitted using the expression of all promoters in the dataset. DNA binding motifs from SwissRegulon28 were used in this analysis.

DNA isolation and DNA methylation analysis

DNA was isolated from interphase after Trizol extraction of RNA according to the manufacturer’s recommendations. DNA methylation assays were performed on biological triplicate samples from 3 healthy donors with the Infinium Human Methylation450 BeadChip (Illumina) according to the manufacturer’s instructions. For quality control results, see supplemental Figure 2. The methylation status of each site is represented by an average β-value. An average β of 1 corresponds to complete methylation and an average β of 0 corresponds to no methylation. Probes with a detection P > .01 in any sample were discarded. Sites were considered to be differentially methylated if the average change in β between 2 cell stages was ≥0.17,29,30 with a Benjamini-Hochberg (BH)–adjusted P ≤ .05 (repeated-measures analysis of variance). The principal component analysis (PCA), heatmap, and clustering were performed using R (http://www.rproject.org), version 2.15.1. The 450K array annotation was used to identify CpG islands, shores, shelves, open sea, and enhancer regions. Sixteen differentially methylated sites and 1 additional CpG were validated in biological triplicate samples of CMP, GMP, and PMC cells using bisulfite pyrosequencing. 450K array data are available at GEO (accession number GSE50797). For details, see supplemental Materials.

Results

Isolation of granulocytic progenitors and mature cells

Bone marrow was obtained from healthy donors and FACS sorted based on surface marker expression to isolate the CMP, GMP, PMC, and metamyelocyte/band-myelocyte/PMN cell populations (Figure 1A; supplemental Figure 1). 4,6 Diamidino-2-phenylindole staining of sorted cells demonstrated the characteristic change of nuclear morphology from rounded nuclei in the immature CMP and GMP cells to increasingly lobular in the PMC and PMN cells (supplemental Figure 1A).

Figure 1

Genomic distribution of DNA methylation in granulopoietic cells. (A) A schematic picture of hematopoiesis showing the position of CMP, GMP, PMC, and PMN cell populations with regard to other hematopoietic cells. Dashed lines represent existence of intermediate populations. (B) Frequency plots showing the distribution of β-values in CGIs, shores (0-2 kb flanking islands), shelves (2-4 kb flanking islands), and open sea (unrelated to islands). DNA methylation was measured with the IlluminaHumanMetylome 450 Beadchip using biological triplicate samples for CMP, GMP, PMC, and PMN cells. Completely unmethylated CpG: β = 0; completely methylated CpG: β = 1.

Distribution of DNA methylation

DNA methylation was assessed using the HumanMethylation450 BeadChip (450K array) from Illumina (for array quality control results, see supplemental Figure 2). This platform probes >480 000 cytosines at single cytosine resolution including sites in 96% of known CGIs, but also CpGs in the flanking shores (0-2 kb from island), shelves (2-4 kb from island), and regions called open sea (>4 kb from nearest island).31 The accuracy of the 450K array has been validated by 2 independent laboratories.31,32 The cell populations were analyzed in biological triplicates, and the methylation level for each cytosine is represented by β-values defined as the ratio of the methylation. Correlations between replicates were >0.98 (supplemental Figure 3), and PCA of the 1000 most variable sites separated samples primarily on cell type and donor, especially CD34+ (CMP and GMP) and CD34 (PMC and PMN) samples (supplemental Figure 4).

Our analysis shows that cytosines located in CGIs are generally hypomethylated, whereas shores show bimodal methylation (Figure 1B). In contrast, cytosines in shelves and open sea are generally hypermethylated. Hence, the methylation profile of the myeloid-granulocytic lineage shows an overall bimodal pattern that has also been observed in other cell types.33,34

Temporally distinct methylation changes during granulopoiesis

Comparing methylation data for the different cell populations, we identified 10 156 sites that change in methylation level during granulopoiesis. Of these differentially methylated sites (DMSs), the majority showed decreased methylation (8973 CpGs) and a smaller set gained methylation (930 CpGs) (Figure 2A-B). Interestingly, 253 CpGs both gained and lost methylation at different stages of cell maturation (henceforth referred to as oscillating DMSs). When all DMSs were used for hierarchical clustering, the samples clustered based on cell stage rather than donor (Figure 2A).

Figure 2

Changes in DNA methylation occur at distinct time points. (A) A total number of 10 156 DMSs were identified and subjected to hierarchical clustering for all cell populations from all 3 donors. C, CMP; G, GMP; P, PMC; N, PMN. DMSs were defined as CpGs with mean Δβ ≥ 0.17 between any 2 cell populations (eg, CMP vs GMP and CMP vs PMN). (B) All 10 156 DMSs were divided into (left) DMS up, (center) DMS down, or (right) DMS oscillating based on direction of methylation change. β-values were plotted between consecutive cell stages, ie, GMP vs CMP, PMC vs GMP, and PMN vs PMC. β-values of DMSs are shown as dots, and β-values for unchanged CpGs are shown as smooth scatter (blue and black dots).

Strikingly, increases in methylation predominantly occur between CMP and GMP, the 2 more immature cell types (Figure 2A-B, left). Some CpGs also show increased methylation in the GMP-PMC transition, whereas very few CpG sites increase at the final stage of differentiation from PMC to PMN. Although reduction of methylation occurs at all stages of granulopoiesis, the far greatest change happens between GMP and PMC (Figure 2A-B, center). The oscillating DMSs generally change methylation level between CMP and GMP and then reverse in the transition to PMC (Figure 2B, right). In general, the changes in methylation are least dramatic between the 2 final stages PMC and PMN after commitment to neutrophil differentiation, suggesting that methylation is involved in lineage choice.

Validation of DMSs

We compared the β-values for our DMSs in PMNs with 450K data from a recent study on peripheral blood cells.35 Our dataset showed exceptionally good correlation with both neutrophils and granulocytes but poor correlation with T cells (Figure 3A; supplemental Figure 5A). This observation validates our DMS values in PMN cells and shows that DNA methylation changes during granulopoiesis persist when the cells enter the bloodstream. Also, the poor correlation with T cells indicates that the DNA methylation changes are neutrophil specific.

Figure 3

DMS validation and distribution of DMSs. (A) Correlation plot comparing the β-values for the 10 156 DMSs in PMNs with published 450K data from peripheral neutrophils.35 The plot shows the average of the 3 PMN replicates from this study vs the average of the 6 neutrophil replicates. Lower density, red; higher density, blue. The Pearson correlation coefficient is given above the plot. (B) Correlation plot comparing the average β-values of 16 DMSs in CMP, GMP, and PMC, with pyrosequencing results (average of biological triplicates). The DMSs were located in the MPO, PRTN3, AZU1, ELANE, GFI1, TBCD, and HDAC4 genes. One additional 450K CpG (probe cg07805029) in GFI1 was also included in pyro analysis but not in the list of DMSs (P = .07). The Pearson correlation coefficient is given above the plot. (C) Bar chart showing the methylation levels of the DMSs plotted in B (including the additional probe cg07805029) in GMP (white) and PMC (gray) measured by bisulfite pyrosequencing. The bars show the average of biological triplicates and error bars display standard deviation. Red dots mark the corresponding average β-value from 450K analysis. The 450K probe ID is stated on the x-axis. (D) The bar chart displays the relative distribution of CpGs over CpG islands, shores, shelves, and open sea. Blue bars represent CpGs present on the 450K methylation array and green bars represent DMSs. The numbers of CpGs are given in the adjacent table. Total number of cytosines on array after filtering for detection P < .01 was 477 112. Total number of differentially methylated CpGs was 10 156. P values were calculated with the χ2 test.

For the CMP, GMP, and PMC datasets, we performed extensive validation using pyrosequencing of 16 DMSs and 1 additional CpG (probe cg07805029 in pyro assay “GFI1 reg5” was excluded from DMS list, P = .07). We observed >95% correlation between pyrosequencing and 450K results (Figure 3B). For all sites, the methylation change between GMP and PMC was highly reproducible (Figure 3C) and extended to adjacent CpGs included in the pyrosequencing but not on the 450K array (supplemental Figure 5B). Thus, in silico and pyrosequencing validation support the accuracy of our methylation results. For detailed information about pyro-assay design, locations, and results, see supplemental Tables S1 to S3 and supplemental Figure 7.

Genomic distribution of differentially methylated cytosines

To examine the distribution of methylation changes over genomic context, we compared the frequencies of sites in CGIs, shores, shelves, and open sea on the platform with the frequencies in the list of DMSs (Figure 3D). DMSs within CGIs were greatly underrepresented: 4% of the probes in the DMS list vs 31.2% present on the 450K array (P < .001 with χ2 test), whereas DMSs were overrepresented in shelves (15.9% vs 9.6%, P < .001) and open sea (57.3% vs 36.1%, P < .001). The frequency of DMSs in shores was approximately the same on the platform and in the DMS set (22.8% vs 23.1%). Thus, methylation appears to be more dynamic outside of CGIs during granulocytic development as has previously been reported for other systems.34

GO analysis of differentially expressed genes reflects terminal differentiation

To investigate changes in gene expression during granulocyte development, we used 2 complementary genome-wide techniques to determine mRNA abundance: microarray analysis and CAGE, in which a short sequence adjacent to the 5′cap is sequenced, giving information about the transcription level and the precise transcription start site (TSS). Transcriptome microarray analysis was performed on 5 biological replicates, 2 of which were also analyzed with CAGE. For validation, we examined the expression for some of the surface markers used in the FACS sorting of the cell populations including CD34, KIT (CD117), and CD11B (supplemental Figures 1A and 6A-B). For both methods, the results match the expected transcription patterns of these genes associated with neutrophil differentiation.

Using the microarray data, we created lists of differentially expressed genes between any pairwise comparison of cell stages (eg, CMP vs GMP or CMP vs PMC; fold change ≥2, P ≤ .05 BH adjusted, repeated-measures analysis of variance). GO analysis showed that down-regulated genes (827 genes) were enriched in GO terms related to biosynthesis, translation, and cell division reflecting terminal differentiation (Figure 4A). The up-regulated genes (704 genes) were significantly enriched in genes involved in neutrophil functions, including the GO terms defense response, response to bacterium, and neutrophil. Seventy-two genes were present in the lists of both up-regulated and down-regulated genes, and this list of oscillating genes was enriched for the terms lysosome and immune response genes (supplemental Figure 6C), likely because many granule proteins, such as elastase and myeloperoxidase, are expressed during a brief period in granulopoeisis20 (Figure 5A-B).

Figure 4

Differentially expressed genes are enriched in neutrophil function-related factors and overlap with DNA methylation changes. (A) Lists of differentially expressed genes were generated from microarray analysis based on a BH adjusted P < .05 and a fold change >2 between any 2 cell populations (eg, CMP vs GMP and CMP vs PMC). The lists of up-regulated genes from all comparisons and the lists of down-regulated genes from all the comparisons were combined. Genes appearing on both lists (oscillating genes) were removed and analyzed separately (supplemental Figure 6C). Gene lists were used for gene ontology analysis using the DAVID Bioinformatics Resources. The tables include selected GO terms. Only terms not present in both lists are shown. BH, Benjmini-Hochberg adjusted P value. (B) Venn diagrams showing the overlap between differentially methylated and differentially expressed genes. Lists of differentially methylated genes include all genes with ≥1 DMS from 1500 bp upstream of TSS to the end of the 3′UTR. Separate lists for genes with increased and decreased methylation were made, excluding genes with oscillating DMSs or >1 DMS changing in separate directions. Lists of differentially methylated genes and differentially expressed genes were filtered for genes present on both platforms and compared for overlap. P values were calculated using the hypergeometric distribution test. METdwn, genes with decreased methylation; METup, genes with increased methylation; EXPdwn, genes with decreased expression; EXPup, genes with increased expression; METall, all differentially methylated genes; EXPosci, genes with oscillating expression.

Figure 5

DNA methylation is reduced when MPO, ELA2, CST7, and PRTN3 expression increases. Transcription and DNA methylation levels of MPO, ELA2, CST7, and PRTN3. (A) Transcription levels measured with microarray analysis. The boxplots show the signals of all 5 replicates (arbitrary units). (B) Transcription levels measured with CAGE analysis (TPM). Mean expression of biological duplicates including tags over the entire gene normalized to the population with highest expression. (C) DNA methylation changes in MPO, ELA2, CST7, and PRTN3. The plot displays the change in β-value of all DMSs identified in each gene (from 1500 bp upstream of TSS to the end of the 3′UTR). Each dot represents 1 DMS.

Changes in gene expression overlap with DNA methylation changes

Lists of differentially methylated genes, ie, genes with ≥1 DMS located between 1500 bp upstream of the TSS and the end of the 3′untranslated region (UTR), were constructed yielding 355 genes with increased methylation (METup) and 3827 genes with decreased methylation (METdwn). There was a small but significant overlap between genes showing decreased methylation and genes with increased expression, as well as between genes with increased methylation and decreased expression (hypergeometric distribution test; Figure 4B). In contrast, there was no statistically significant overlap between increased expression and methylation or decreased expression and methylation. Therefore, our results support an inverse relationship between DNA methylation and gene expression. Genes with oscillating expression also showed a significant overlap with changes in methylation.

Changes in methylation occur concomitantly with changes in gene expression

Next, we investigated the timing of DNA methylation and expression changes with regard to each other. The genes encoding the azurophilic granule proteins myeloperoxidase (MPO), elastase (ELA2), and proteinase 3 (PRTN3) are expected to be expressed in the PMC stage.20 Expression of these genes peaked in PMC (Figure 5A-B), and CpG methylation levels concomitantly decreased (Figure 5C). The same was true for CST7 encoding a protease inhibitor with immunoregulatory functions (Figure 5). Collectively, this demonstrates an overlapping timing of expression and DNA methylation changes on single gene level.

DNA methylation of key hematopoietic transcription factors correlate with expression and motif activity

In addition to measuring transcription levels, CAGE analysis identifies the precise TSS and consequently gene promoters (method reviewed in ref. 36). MARA assesses the transcriptional activity of genes with the DNA sequence motif recognized by a specific transcription factor in proximity of the TSS.21 Using the CAGE data, we screened for transcription factors that showed changes in motif activity during granulopoiesis and cross-referenced these with transcription and DNA methylation data. We found several transcription factors with changes in motif activity that also showed changes in expression and reciprocal changes in DNA methylation. Notably, these transcription factors include PU.1 (SPI1), GATA2, GFI1, and ETS1, all of which are important in regulation of hematopoiesis,37-40 and MEF2D, which is implicated in macrophage differentiation.41 The expression of PU.1, a regulator of myeloid development, increases gradually with differentiation and shows reciprocal changes in DNA methylation while the motif activity, ie, activity of PU.1 target genes, increases (Figure 6). GATA2 is expressed in the CMP population. On differentiation, the GATA2 gene shows increased methylation, decreased expression, and a decrease in GATA motif activity (Figure 6). Similarly, increased expression of GFI1, encoding a transcriptional repressor important for neutrophil development and function,38,40 is accompanied by decreased methylation and decreased motif activity (Figure 6). Similar trends were observed for ETS1 and MEF2D (supplemental Figure 8). These results suggest that DNA methylation may regulate the expression, and thereby the activity, of several important hematopoietic transcription factors.

Figure 6

DNA methylation of hematopoietic transcription factors is reciprocal to changes in expression and motif activity. DNA methylation changes, transcription, and motif activity of PU.1, GATA2, and GFI1. (A) Graphs showing the β-values of DMSs annotated to respective genes. Each dot represents 1 DMS. (B) Box plots displaying the signals from expression microarray analysis for all 5 replicates (arbitrary units). (C) Bar chart demonstrating the relative motif activity for all CAGE samples (arbitrary units), ie, transcriptional activity of genes with promoters containing potential binding sites for PU.1, GATA2, or GFI1. Error bars indicate standard deviation. The GATA2 binding motif is shared with GATA1 and GATA3.

Methylation is more dynamic in enhancer elements

On the 450K array, 100 377 sites are annotated to regions predicted to be enhancers based on their chromatin state.31,42,43 We compared the variation of probe signals within and outside of enhancers across all samples divided into the categories CGIs, shores, shelves, and open sea (Figure 7A). For all regions, the variation in methylation level was greater within than outside of enhancers (P < 1 × 10−15, Student t test). In addition, CpGs in enhancer regions are significantly enriched in the list of DMSs (3774 DMSs, P < 2.2 × 10−16, Fisher’s exact test), further supporting the observation that enhancer regions display dynamic DNA methylation during granulopoiesis.

Figure 7

Methylation is more variable in enhancer elements and DMSs are enriched in differentiation-induced neutrophil enhancers. (A) Box plot demonstrating the β-value variation of CpGs investigated with the 450K methylation array. All sites present on the 450K platform (filtered for detection P < .01) were divided in groups based on location in CpG island, shore, shelf, or open sea. Each group was further divided based on location in enhancer regions (white) or outside enhancer regions (gray). The log10 of the standard deviations of probe β-values for all replicates were plotted. P values (Student t test) for comparison between enhancer and not enhancer for each group are indicated between boxes. E, enhancer; noE, not enhancer; SD, standard deviation. (B-D) Bar graphs displaying the relative distribution of 450K methylation array probes/DMSs in enhancers identified by bidirectional transcription in CAGE analysis. Overlapping 450K array probes, ie, percentages of 450K probes located in enhancers, are shown in dark gray (total 477 112 probes) and percentages of DMSs located in enhancers are shown in light gray (total 10 156 probes). Absolute numbers are given at the top of each bar. P values were calculated by Fisher’s exact test. (B) Overlap in all enhancers identified as being active (expressed) in neutrophils. (C) Comparing enhancer data for CD34+ cells vs neutrophils. The graph displays the relative number of overlapping probes in enhancers specifically expressed/up-regulated in (left) CD34+ cells or in enhancers specifically expressed/up-regulated in neutrophils. (D) Same as in C, but for the CD133+ vs neutrophils comparison. (E) Expression (ie, bidirectional transcription) changes of enhancers with decreased methylation in neutrophils compared with progenitor cells. A total of 294 CAGE-derived enhancers active in CD133+ cells, CD34+ cells or neutrophils overlapped with DMSs showing decreases in β-value during differentiation. The expression of these enhancers was measured by CAGE. Box plot displays log2 fold change between neutrophils and (left) CD133+ cells or (right) neutrophils and CD34+ cells. P values were 1.16 × 10−35 and 4.24 × 10−41, respectively (1-sided Wilcoxon signed-rank test). (F) Expression changes of genes associated with differentially methylated enhancers. A total of 105 of the enhancers in E were associated with genes. Box plot shows expression (log2 fold change) of genes associated with enhancers with DMSs that decrease during differentiation in PMN cells compared with CMP, GMP, and PMC cells. P values were 1.031 × 10−10, 1.694 × 10−10, and 2.17 × 10−9 using the 1-sided Wilcoxon signed-rank test based on the average of replicates. Expression was measured by microarray analysis of 5 biological replicate samples. Genes were required to have 3 of 5 replicates with a detection P < .01 in ≥1 cell type.

We also investigated whether DMSs were preferentially located in enhancers defined by an alternative method based on the observation that active enhancers are often bidirectionally transcribed.27 The enhancer RNA abundance is a strong indicator of enhancer activity and can be used to determine cell-specific enhancer use.27 The strength of this CAGE-based method compared with the annotated chromatin-defined enhancers used above is that we can assess enhancer activity in the actual cell lineage used in this study. First, we considered enhancers expressed in neutrophils and found that there was a highly significant overlap between these regulatory elements and our DMS set (P < 1 × 10−50, Fisher's exact test; Figure 7B). We then looked closer at DMS-overlapping enhancers that change activity during granulopoiesis. Enhancers were identified by bidirectional expression in CAGE libraries (“Materials and methods”; supplemental Materials) from 2 immature cell populations, CD34+ and CD133+, and neutrophils. Hierarchical clustering based on enhancer expression readily separated progenitors from neutrophils, confirming cell type–specific enhancer activity (supplemental Figure 9). We found a significant enrichment of DMSs in enhancers that are specifically or significantly higher expressed in neutrophils compared with enhancers active in progenitor populations (Figure 7C-D). Most of the progenitor and neutrophil CAGE-derived enhancers that overlapped with DMSs (294 of 311) showed a decrease in methylation during differentiation. When comparing expression levels of these 294 enhancers in progenitor cells and in neutrophils, we observed increased transcription in neutrophils indicating that enhancers that decrease in methylation during granulopoiesis concomitantly increase their activity (Figure 7E). Of these enhancers, 105 had a significant CAGE expression correlation (BH false discover rate ≤1 × 10−5, Pearson correlation test across 88 myeloid FANTOM5 CAGE libraries) to ≥1 RefSeq TSS within 500 kb, suggesting regulatory interaction.27 Corroborating this, these enhancer-associated genes were significantly higher expressed (P < 2.2 × 10−9, 1-sided Wilcoxon signed-rank test) in PMN cells compared with the earlier populations (Figure 7F). Interestingly, several of these genes encoded factors important for neutrophil development and function, such as CSF2Rb, IL1b, IL8Ra/b, TLR1, and TLR6.

We conclude that DMSs are overrepresented in enhancers that are active during granulopoiesis and that enhancers with decreased methylation show increased expression. Furthermore, genes associated with demethylated enhancers show increased expression, suggesting a gene regulatory role for enhancer methylation.

Discussion

CpGs distant from islands show most changes in DNA methylation during granulopoiesis

Most previous studies of DNA methylation in myeloid development have been performed in the murine system and/or analyzed few or no intermediate cell types between multipotent cells and terminally differentiated neutrophils. Although there are many similarities between human and murine hematopoiesis, there are also known differences (reviewed in ref. 44). We present data from human granulocytic development with high cellular resolution yielding detailed information of DNA methylation changes throughout differentiation.

Promoter proximal CGIs have been the primary focus in DNA methylation research. However, in agreement with the massive under-representation of DMSs in CGIs in our study, it has recently been proposed that methylation changes in development and malignancies mainly occur outside of CGIs, especially in the shores.9,10 Although more frequent in shores than in CGIs, we see no over-representation of DMSs in shore regions. In fact, the greatest number of changes occurs in the shelf and open sea regions, >2 kb from the nearest CGI.

DNA methylation changes coincide with cell lineage restriction

Previous studies have demonstrated a general loss of methylation on CMP entry, as opposed to an increase on CLP entry.6,8,13 This trend continues throughout granulopoiesis, as we observe predominantly decreased DNA methylation in the more mature cell populations. Moreover, the methylation changes in granulopoiesis appear to be distinct in time. Our analysis shows that gain of DNA methylation mainly occurs between CMP and GMP, whereas the greatest reduction occurs between GMP and PMC. Between PMC and PMN, changes are comparatively minute. In agreement with our data, recent articles report that neutrophils do indeed have more hypomethylated regions than HSCs and lymphoid cells.11,13,14 However, the intermediate cell stages analyzed were fewer than in this study, and the exact time of methylation loss was not established.

It is striking that the DNA methylation preferentially changes at points of lineage restriction with the greatest change occurring on loss of oligopotency between GMP and PMC. This suggests a role of DNA methylation in regulating cell plasticity. Mice with hypomorphic DNMT1 have decreased HSC self-renewal capacity, implicating DNA methylation in maintenance of multipotency.6 Furthermore, these mice were protected against induced AML, indicating a requirement of DNA methylation for plasticity and self-renewal of leukemic cells. AML is caused by malignant transformation of cells mainly corresponding to the GMP population. It is tempting to speculate that the loss of DNMT1 mimics the loss of methylation that we observed in this study between GMP-PMC. If so, the protection against AML in the DNMT1-impaired mice could be caused by an artificial reduction of methylation at this stage, forcing the cells to differentiate and thereby limiting self-renewal potency and preventing leukemic transformation. Recently 5-hydroxymethylation has gained interest both as a demethylation intermediate and as a potential regulatory mark.45 As the 450K methodology is based on bisulfite treatment, we cannot distinguish between DNA methylation and hydroxylmethylation, but such analysis would be highly interesting.

DMSs are enriched in neutrophil activated enhancer elements

CpGs located in putative enhancer regions were significantly enriched in our DMS list, both based on the 450K array annotation of chromatin-derived enhancers and on enhancers active in neutrophils identified by bidirectional transcription (Figure 7). Importantly, DMSs were especially enriched in enhancers that become increasingly active when cells differentiate from progenitor stage to neutrophils. Enhancers with decreased methylation showed increased expression. This suggests that DNA methylation may be involved in repression of neutrophil specific enhancers in progenitor cells and regulation of differentiation and lineage-specific enhancer activity. This is supported by a finding in T cells, where investigation of DMRs at cell type–specific genes revealed methylation-sensitive enhancer activity.46 It was suggested that methylation contributes to lineage restriction by regulating lineage-associated enhancer elements. In addition, another group reported that several putative enhancer elements were specifically unmethylated in T cells.11

DNA methylation changes are coordinated and overlap with expression changes

We observed a small but significant concordance between changes in DNA methylation and in gene expression, supporting the general view of DNA methylation as a repressive modification. Importantly, we observe altered DNA methylation levels in genes coding for several key hematopoietic transcription factors, including ETS1, PU.1, and GATA2, with reciprocal changes in transcription. Moreover, these changes were reflected by the transcriptional activity of these factors’ target genes, suggesting a role of DNA methylation in regulating the important transcription factor circuitry that governs hematopoiesis. Indeed, methylation of the genes and binding sites of myeloid transcription factors in lymphoid cells has been proposed to function as a safeguard against lineage inappropriate activation.11 For several of the genes, such as MPO, ELA, PRTN3, and CST7 (Figure 5), reduction of expression at the late stages of granulopoiesis is not followed by increased DNA methylation. This could be due to alternative modes of epigenetic silencing (eg, histone deacetylation or methylation) or down-regulation of controlling transcription factors.38

We found a small but significant overlap between genes with decreased in DNA methylation and increased gene expression and vice versa (Figure 4B). However, it has been suggested that while anticorrelation is true for promoter DNA methylation, intragenic DNA methylation relates positively to transcriptional activity.33,47,48 Although roughly half of the gene-associated DMSs identified in our study are located after the first exon in the gene body or in the 3′UTR (data not shown), we fail to detect significant overlaps between increases/decreases of both transcription and DNA methylation. This is in agreement with 2 recent studies showing anticorrelation either irrespective of CpG position46 or when located in intragenic CGIs.7 An additional study showed that the best correlation with expression was with hypomethylation of DNA around 1000 bp downstream of the TSS.13

In conclusion, the current study highlights the importance of using a high cell stage resolution concerning cell maturation when studying differentiation. By examining several closely related cell populations in myeloid development, we were able to pinpoint the timing of DNA methylation changes and show that it coincides with loss of oligopotency accompanying the GMP-PMC transition. DNA methylation is generally more dynamic in enhancer regions and in areas outside of CGIs during granulopoiesis. There is a strong inverse relationship between changes in DNA methylation and gene expression, and these changes are coordinated in time. We also find that the expression and activity of several key hematopoietic transcription factors and differentiation activated neutrophil enhancers may be regulated by DNA methylation. These data suggest an important regulatory role for DNA methylation in granulopoiesis.

Authorship

Contribution: M.R. did the experiments, the analysis of expression, and 450K arrays and wrote the manuscript; T.O., I.D., and S.L. isolated the cells; R.A., I.H., and A.S. analyzed and characterized the CAGE enhancers; M.K. supervised the pyrosequencing validation; M.I. and S.N-S. were responsible for CAGE data production; T.L. was responsible for CAGE tag mapping; H.K. managed the data handling; P.C., Y.H., and A.R.R.F. were responsible for FANTOM5 management and concept; E.A. and M.d.H. performed additional analysis of CAGE and methylation array data; S.L., R.A., I.H., A.S., E.A., K.E., and A.L. contributed to the manuscript writing; and E.A., K.E., and A.L. initiated, planned, and supervised the study.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Andreas Lennartsson, Karolinska Institute, Center of Biosciences, Dept. of Biosciences and Nutrition, Halsovagen 7-9, 141 57 Huddinge, Sweden; e-mail: andreas.lennartsson{at}ki.se; Karl Ekwall, Karolinska Institute, Center of Biosciences, Dept. of Biosciences and Nutrition, Halsovagen 7-9, 141 57 Huddinge, Sweden; e-mail: karl.ekwall{at}ki.se; and Erik Arner, RIKEN Center for Life Science Technologies, Division of Genomic Technologies, RIKEN Yokohama Institute, 1-7-22 Suehiro-cho, Tsurumi-ku, Yokohama, Kanagawa 230-0045, Japan; e-mail: arner{at}gsc.riken.jp.

Acknowledgments

The authors thank Bodil Rosberg, Stefan Deneberg, Christer Nilsson, and Sofia Bengtzen for technical assistance with isolation of primary cell populations; Per Lönnbro for nuclear staining; and the Bioinformatics and Expression Analysis (BEA) core facility for help with Illumina arrays. The authors thank all members of the Functional Annotation Of The Mammalian Genome 5 (FANTOM5) consortium for contributing to generation of samples and analysis of the dataset and thank Gene Network Analysis Service (GeNAS) for data production.

M.R., A.L., and K.E. were supported by grants from the Swedish Cancer Foundation, the Swedish Children Cancer Foundation, the Swedish Research Council (VR), the Göran Gustafsson Foundation for Research in Natural Sciences, the Sigurd and Elsa Goljes Foundation, the Åke Olsson Foundation for Hematology, and the Knut and Alice Wallenberg Foundation. R.A., I.H., and A.S. were supported by funds from FP7/2007-2013/ERC grant agreement 204135, the Novo Nordisk foundation, the Lundbeck Foundation, and the Danish Cancer Society. FANTOM5 was supported by a Research Grant from the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) to the RIKEN Center for Life Science Technologies, RIKEN Preventive Medicine and Diagnosis Innovation Program (to Y.H.), RIKEN Omics Science Center from MEXT (to Y.H.), and a Grant of the Innovative Technology (Cell Innovation Program) from MEXT to (Y.H.).

Footnotes

  • K.E., E.A., and A.L. contributed equally to this work.

  • *RIKEN Omics Science Center ceased to exist as of April 1, 2013 due to RIKEN reorganization.

  • This article contains a data supplement.

  • There is an Inside Blood Commentary on this article in this issue.

  • The publication costs of this article were defrayed in part by page charge payment. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

  • Submitted February 4, 2013.
  • Accepted October 15, 2013.

References

View Abstract