Methylome Profiling Reveals Distinct Alterations in Phenotypic and Mutational Subgroups of Myeloproliferative Neoplasms

The myeloproliferative neoplasms (MPN), essential thrombocytosis, polycythemia vera, and primary myelofibrosis (PMF) share the same acquired genetic mutation, JAK2V617F, but vary with respect to epidemiology, disease phenotype, and prognosis. PMF carries the worst prognosis within the MPN due to complications of bone marrow failure or leukemic transformation. Within the JAK2V617F-positive MPNS, variation in gene dosage of JAK2V617F due to varying rates of chromosome 9p uniparental disomy (UPD) offers some rationale as to these discrepancies in phenotype, yet significant overlap still occurs between polycythemia vera and myelofibrosis. Acquired genomic lesions in other signal transduction pathway genes, JAK2 Exon 12, MPL, and LNK, recapitulate myeloproliferative phenotypes but still do not segregate phenotypes within the MPNs. More recently, lesions in genes central to epigenetic regulation, TET2 and ASXL1, have been found to be prevalent within the MPNs, and in the case of ASXL1, segregate with myelofibrosis phenotypes, implicating the importance of epigenetic regulation as a novel pathway relevant to myelofibrosis biology (1, 2)

Recent evidence suggests that MPNs are characterized by aberrant transcriptional profiles, and some of these changes are driven by the mutant JAK2V617F kinase (3). One of the ways that gene expression may be dysregulated is through aberrant DNA methylation. Methylation of cytosine has been implicated as a way to silence genes epigenetically and indicates an attractive target for potential therapeutics. Single locus studies have shown aberrant methylation of promoters of genes such as SOC1, SOCS3 and CXCR4 can be seen in MPNs, although these changes have not been studied at a genome-wide level (4, 5). Genome-wide methylome studies in other disease have revealed surprising patterns and have the potential to change existing paradigms (6, 7). The recent discovery of direct epigenetic effects of the JAK2 mutation on histone H3 phosphorylation further point to a role of epigenetic disturbances in the pathogenesis of these diseases. (8). Thus, the purpose of this study was to assess genome-wide patterns of DNA methylation across the MPNs stratified by disease class and to determine the effect of JAK2V617F on the methylome.

The study population consisted of 29 patients with MPNs and controls evaluated at either Johns Hopkins University (Baltimore, MD) or Albert Einstein School of Medicine (Bronx, NY). The controls were aged 40, 79, and 83 years, respectively. The study was approved by an Institutional Review Board, and all patients gave written consent. Genomic DNA from neutrophils isolated from peripheral blood via Ficoll gradient density separation was prepared using Qiagen reagents.

The HELP assay was conducted as previously published (9). Intact DNA of high molecular weight was corroborated by electrophoresis on 1% agarose gel in all cases. One microgram of gDNA was digested overnight with either HpaII or MspI (NEB). On the following day, the reactions were extracted once with phenol–chloroform and resuspended in 11 μL of 10 mmol/L Tris-HCl, pH 8.0, and the digested DNA was used to set up an overnight ligation of the JHpaII adapter using T4 DNA ligase. The adapter-ligated DNA was used to carry out the PCR amplification of the HpaII- and MspI-digested DNA as previously described (9). Both amplified fractions were submitted to Roche-NimbleGen, Inc. for labeling and hybridization onto a human hg17 custom-designed oligonucleotide array (50-mers) covering 25,626 HpaII-amplifiable fragments (HAF) located at gene promoters. HAFs are defined as genomic sequences contained between 2 flanking HpaII sites found within 200 to 2,000 bp from each other. Each fragment on the array is represented by 15 individual probes distributed randomly spatially across the microarray slide. Thus, the microarray covers 50,000 CpGs corresponding to 14,000 gene promoters. HELP microarray data have been submitted to the Gene Expression Omnibus (GEO) database for public access (GSE42721).

Validation of HELP microarray findings was conducted by matrix-assisted laser desorption/ionization–time-of-flight (MALDI-TOF) mass spectrometry using EpiTyper by massARRAY (Sequenom) on bisulfite-converted DNA as previously described (10, 11). massARRAY primers were designed to cover the flanking HpaII sites for a given HAF, as well as any other HpaII sites found up to 2,000 bp upstream of the downstream site and up to 2,000 bp downstream of the upstream site, to cover all possible alternative sites of digestion. Primers are available on request.

All microarray hybridizations were subjected to extensive quality control using the following strategies. First, uniformity of hybridization was evaluated using a modified version of a previously published algorithm (12) adapted for the NimbleGen platform, and any hybridization with strong regional artifacts was discarded and repeated. Second, normalized signal intensities from each array were compared against a 20% trimmed mean of signal intensities across all arrays in that experiment, and any arrays displaying a significant intensity bias that could not be explained by the biology of the sample were excluded.

Signal intensities at each HAF were calculated as a robust (25% trimmed) mean of their component probe-level signal intensities. Any fragments found within the level of background MspI signal intensity, measured as 2.5 mean absolute differences (MAD) above the median of random probe signals, were categorized as “failed.” These “failed” loci therefore represent the population of fragments that did not amplify by PCR, whatever the biologic (e.g., genomic deletions and other sequence errors) or experimental cause. On the other hand, “methylated” loci were so designated when the level of HpaII signal intensity was similarly indistinguishable from background. PCR-amplifying fragments (those not flagged as either “methylated” or “failed”) were normalized using an intra-array quantile approach wherein HpaII/MspI ratios are aligned across density-dependent sliding windows of fragment size–sorted data. The log₂(HpaII/MspI) was used as a representative for methylation and analyzed as a continuous variable. For most loci, each fragment was categorized as either methylated, if the centered log HpaII/MspI ratio was less than zero or hypomethylated if, on the other hand, the log ratio was greater than zero.

Unsupervised clustering of HELP data by hierarchical clustering was conducted using the statistical software R version 2.6.2. A 2-sample t test was used for each gene to summarize methylation differences between groups. Genes were ranked on the basis of this test statistic, and a set of top differentially methylated genes with an observed log fold change of >1 between group means was identified. Genes were further grouped according to the direction of the methylation change (hypo- vs. hypermethylated), and the relative frequencies of these changes were computed among the top candidates to explore global methylation patterns. Validations with massARRAY showed good correlation with the data generated by the HELP assay. massARRAY analysis validated significant quantitative differences in methylation for differentially methylated genes selected by our approach.

Using the Ingenuity Pathway Analysis (IPA) software, we conducted an analysis of the biologic information retrieved by each of the individual platforms alone and compared it with the information obtained by the integrated analysis of all 3 platforms. Enrichment of genes associated with specific canonical pathways was determined relative to the Ingenuity knowledge database for each of the individual platforms and the integrated analysis at a significance level of P < 0.01. Biologic networks captured by the different microarray platforms were generated using IPA and scored on the basis of the relationship between the total number of genes in the specific network and the total number of genes identified by the microarray analysis. The list of hypermethylated genes was examined for enrichment of conserved gene–associated transcription factor–binding sites using the Molecular Signatures Database (MSigDB; ref. 13). Their functional gene sets were obtained from Gene Ontology (GO; ref. 14).

This analysis was conducted by Gene Set Enrichment Analysis (GSEA; ref. 13), a computational method that determines whether an a priori defined set of genes (commonly hypermethylated genes in myelodysplastic syndrome) shows statistically significant, concordant differences between 2 biologic states. Same method was applied to determine whether transcription-binding sites are randomly distributed in the differentially methylated genes. The a priori defined gene sets used in this analysis is transcription factor target (TFT), which contains genes that share a transcription factor–binding site defined in the TRANSFAC database (15). Using GSEA ‘Pre-ranked’ algorithm, 1,000 permutations were applied to sample labels to test whether genes from each TFT gene sets were randomly distributed along the differentially methylated gene list.

Neutrophil isolation and DNA preparation were conducted as previously described (16). The JAK2V617F neutrophil allele burdens were measured using an allele-specific, quantitative real-time PCR assay sensitive to a lower limit of detection of 5% of either the wild-type or mutant JAK2 allele as previously described (16). Single-nucleotide array karyotyping assay and analysis were conducted as described previously (17) Briefly, Gene Chip Mapping Affymetrix 250K arrays (Affymetrix) were used for single-nucleotide polymorphism karyotyping (SNP-K) analysis and used per the manufacturer's instructions. Signal intensity was analyzed and SNP calls determined using Gene Chip Genotyping Analysis Software Version 4.0 (GTYPE). Copy number and areas of UPD were investigated using a hidden Markov model and copy number Analyzer for Affymetrix GeneChip Mapping 250K arrays (CNAG v3.0). Exon 12 of ASXL1 was amplified from neutrophil genomic DNA and resequenced as previously described (18) For TET2 mutations, the all the exons of the gene were PCR-amplified and sequenced using similar methodology. TET2 amplification and sequencing primers have been described previously (19). Identification of known SNPs was established via searching the NCBI dbSNP database and the 1000 Genomes Project.

Levels of methylated and hydroxymethylated DNA were assessed using the Methyl Flash Methylated/Hydroxymethylated DNA Quantification Kit (Colorimetric, catalog No. P-1034/P-1036; Epigentek), respectively. Briefly, 100 ng (methylation assay) and 200 ng (hydroxymethylation assay) of genomic DNA were assayed according to the manufacturer's instructions as previously described (20). Experiments were carried out in triplicates. Adsorption at 450 nm was read on a Versamax Tunable Microplate Reader (Molecular Devices). Groups were compared using the Mann–Whitney U test. P < 0.05 was considered statistically significant. Statistics were conducted using PASW/SPSS Version 18 (IBM).

Cell lines HL60, K562, and HEL were purchased from American Type Culture Collection (ATCC) and SET-2 was obtained from DSMZ. These had been authenticated at ATCC and DSMZ. All cell lines were cultured in RPMI medium (Invitrogen) supplemented with 10% heat-inactivated FBS and penicillin (100 U/mL), streptomycin (100 mg/mL), and 4 mmol/L glutamine. Cells were maintained at 37°C and humidified with 95% air and 5% CO₂ for cell culture. For the viability assays, the cells were cultured in 0.5, 1, and 5 μmol/L decitabine (Sigma) for 2 days. Decitabine was added to the culture daily, dimethyl sulfoxide (DMSO) served as control. Viability was measured on day 3 using MTT assay (Promega) according to the manufacturer's instruction.

The gDNA (200–500 ng) was cleaved with HpaII + EcoRI or MspI + EcoRI in 2 separate 20-mL reactions containing 33 mmol/L Tris-acetate, 10 mmol/L Mg acetate, 66 mmol/L K acetate, pH 7.9, 0.1 mg/mL bovine serum albumin, and 5 units of each restriction enzymes. The reactions were set up in a 96-well format and incubated at 37°C for 4 hours. Then, 20-mL annealing buffer (20 mmol/L Tris-acetate, 2 mmol/L Mg acetate, pH 7.6) was added to the cleavage reactions, and samples were placed in a PSQ96TMMA system (Biotage AB). The instrument was programmed to add dNTPs in 4 consecutive steps including step 1: dATP (the derivative dATPaS is used as it will not react directly with luciferase and prevents nonspecific signals); step 2: mixture of dGTP + dCTP; step 3: dTTP; and step 4: mixture of dGTP + dCTP. Peak heights were calculated using the PSQ96TMMA software. The HpaII/EcoRI and MspI/EcoRI ratios were calculated as (dGTP + dCTP)/dATP for the respective reactions. The HpaII/MspI ratio was defined as (HpaII/EcoRI)/(MspI/EcoRI) (6)

We conducted the HELP assay to study genome-wide methylation patterns in primary polycythemia vera, essential thrombocytosis, and PMF samples. The HELP assay uses differential methylation specific digestion by HpaII and MspI followed by amplification, 2-color labeling and hybridization to quantitatively determine individual promoter methylation of 50,000 CpGs loci covering 14,000 promoters (9, 21). We selected neutrophils to study for many reasons. First, neutrophils are obtained from all patients with a more than 95% purity, and abundant DNA is available. Second, in the MPNs, neutrophils are highly enriched for the malignant clone, more so than purified progenitor cells (16). Neutrophils are a fully differentiated tissue, representing a single hematopoietic lineage, such that differences in lineage representation that may exist in a marrow sample would not confound this analysis. Analysis of 26 MPN neutrophil samples comprising 9 cases of essential thrombocytosis, 5 cases of polycythemia vera, and 12 cases of PMF was conducted and compared with healthy controls (Table 1). Unsupervised clustering based on global methylation profiles showed that whereas polycythemia vera and essential thrombocytosis cases were more similar to the normal controls, the PMF cases were epigenetically distinct from these groups. Five samples of PMF formed a cluster with similar methylation profiles, whereas the rest of the samples exhibited greater heterogeneity. Interestingly, the PMF samples that clustered tightly all had mutated or deleted ASXL1, suggesting that disruption of this gene was driving the epigenetic similarity between these samples (Fig. 1). Overall, methylation profiling of MPNs did not correlate with demographics of these patients. These samples were also examined by high-resolution SNP-array karyotyping, and the methylation patterns did not correlate with the presence and absence of cytogenetic alterations found in these patients.

To further analyze the epigenetic differences between these subtypes of MPNs, we conducted supervised clustering of polycythemia vera and essential thrombocytosis cases and compared them with controls. We saw that these samples had 141 genes that were uniformly hypermethylated compared with controls (Fig. 2A). Bioinformatic analysis revealed that hypermethylated genes displayed specific genomic characteristics and were enriched for binding sites for GATA1 and other transcription factors (Table 2) showing potential regulatory disturbances in these pathways. Genes that were significantly hypermethylated included several novel candidates such as transcription factor HNF4-α, histone acetyltransferase MYST2, interleukin-1, and others (Supplementary Table S1). Functional pathway analysis revealed that these genes are involved in pathways regulated by the NF-κB and HNF4-α transcription factors (Supplementary Fig. S1).

PMF was characterized by both aberrantly hypermethylated (n = 162) and hypomethylated (n = 106) loci when compared with controls (Fig. 2B). These changes were further evaluated by the luminometric methylation assay (LUMA), a quantitative assay that interrogates all HpaII sites in the genome (22). LUMA also revealed significantly more hypomethylation in PMF than in polycythemia vera and essential thrombocytosis cases (49% mean hypomethylation in PMF vs. 38% in polycythemia vera/essential thrombocytosis, P < 0.05, t test; Fig. 2C). The findings of the HELP assay were validated by bisulfite massARRAY analysis and showed good correlation (Supplementary Fig. S2A–S2C). Important functional pathways were found to be affected by aberrant methylation in PMF. The pathways affected by hypomethylated genes included cell signaling, hematopoiesis, and immunologic pathways (Supplementary Fig. S3A), whereas the pathways affected by hypermethylated genes included those governing inflammatory responses and others (Supplementary Fig. S3B and Supplementary Tables S2 and S3 listing all significantly hyper- and hypomethylated genes).

These changes also affected specific genomic regions, and the hypermethylated genes were enriched for binding sites for LEF1 and POU6F1 transcription factors, whereas the hypomethylated genes were enriched for the TFAP2A and NF1 transcription factors (Tables 3 and 4). Interestingly, the significantly, aberrantly hypomethylated regions in PMF were found to be preferentially located outside of CpG islands (Fig. 2D) further showing the unique nature of these changes in PMF.

Because we observed that the ASXL1-mutated/deleted samples were epigenetically similar (Fig. 1), we wanted to determine how the methylation profiles were impacted by disruptions of this gene. Supervised analysis revealed that cases of PMF with mutated or deleted ASXL1 were relatively more hypermethylated than the cases with no disruptions in ASXL1 (Fig. 3A and B). The genes that were uniformly hypermethylated in ASXL1-mutated/deleted cases in comparison to controls included various important candidates such as NPM2, RUNX1, HOXB3, SMAD3, and others (listed in Supplementary Tables S4 and S5) These results reveal that disruption in functioning of chromatin-binding protein, ASXL1, is associated with specific methylation patterns in PMF.

TET2 is an important regulator of hydroxymethylation and subsequent demethylation and has been shown to be mutated in a notable proportion of MPNs (23). Because the effects of TET2 mutation on the methylome of MPNs have not been studied, we sequenced this gene in all MPN samples and correlated it with the methylation patterns in these samples. Supervised analysis revealed that cases of MPNs with mutated TET2 were aberrantly hypermethylated when compared with controls (Fig. 4A and B) and had significantly decreased hydroxymethylation levels (Fig. 4C). The TET2-mutant cases were characterized by a distinct set of commonly methylated genes that included HOXB3, HOXD3, TRAF6, MYST, and others (Supplementary Tables S6 and S7).

The JAK2V617F mutation is seen in a large proportion of patients with MPNs, and recent evidence has shown that the mutant JAK2 kinase can translocate to the nucleus and can directly affect the histone epigenetic machinery (8). Thus, we wanted to determine whether this mutation led to any effects on the methylome of these patients. Because disease class (PMF vs. polycythemia vera or essential thrombocytosis) was the strongest determinant of epigenetic clustering and unsupervised clustering did not show any epigenetic clustering due to JAK2V617F mutation (Fig. 1), we next examined the polycythemia vera and essential thrombocytosis cohort individually by unsupervised clustering. We did not observe any epigenetic clustering based on the presence and magnitude of JAK2 mutation even in this cohort (Supplementary Fig. S4). Supervised clustering could also not show any significant gene-specific methylation signature in the mutant cases. Because patients with MPNs have multiple genetic abnormalities, these data suggest that the JAK2V617F mutation does not exert dominant effect on the methylome of primary samples.

Having shown that MPNs are characterized by aberrant hypermethylation, we next wanted to determine whether these cells were sensitive to growth inhibition by DNMT inhibitors. We compared leukemic cell lines that were derived from patients with MPN (HEL, SET-2, both have JAK2V617F mutation and HEL has ASXL1 deletion; ref. 24) to non-MPN cell lines K562 and HL60. Genome-wide methylation was found to be significantly more in HEL and SET-2 cell lines when assessed with the LUMA assay (Fig. 5). Treatment with decitabine also resulted in significant growth inhibition in both HEL and SET-2 cells at all doses, thus showing functional importance of hypermethylation observed in these MPN-derived cell lines.

The discovery of the most prevalent genomic lesion specific to the chronic MPN, JAK2V617F, has provided important insights into pathogenesis and treatment of these diseases. However, even within the JAK2V617F-positive MPNs, variation in disease phenotype and natural history exists, suggesting that epigenetic processes are significant modifiers of disease pathobiology. Our results show widespread alterations in DNA methylation in MPNs and show that PMF is epigenetically distinct from polycythemia vera and essential thrombocytosis. These differences were not only seen at the involvement of different sets of aberrantly methylated genes but were also seen as a global level by the LUMA assay that showed increased hypomethylation in PMF.

Analysis of differentially methylated genes in MPNs revealed that they were enriched for binding sites for various transcription factors that have important roles in hematopoiesis. GATA1 has been shown to play a role in the pathogenesis of MPNs and myeloid leukemias (25) and was significantly associated with differentially methylated genes in polycythemia vera and essential thrombocytosis. NF1 was found to associated with differentially methylated genes in PMF and was recently found to be deleted and mutated in PMF samples in a large SNP array study (26). These changes show that the changes in methylome occur at specific genomic loci and may be driven by altered levels of transcription factors that are dysregulated in these diseases. In addition to these known transcription factors, our study reveals numerous other differentially methylated loci that are associated with binding sites of transcription factors that have not been studied in MPNs and can potentially be involved in its pathobiology.

We also determined that even though PMF was an epigenetically heterogeneous subgroup, cases of PMF with ASXL1 deletions/mutations were epigenetically similar and clustered tightly together. The ASXL1 gene has been shown to be mutated/deleted in 36% of cases of PMFs and is associated with a more severe clinical presentation (18). Thus, even though it appears that this gene is important in pathobiology of PMF, it is not clear how the disruption of ASXL1 leads to altered gene transcription. ASXL1 is a member of the polycomb repressor complex and is a chromatin-modifying protein (27). Our data show for the first time that deletion/mutations of this protein are associated with distinct signature of DNA methylation involving many important gene promoters. Relative hypermethylation seen in the cases of PMF with ASXL1 disruption also raise the possibility of potential therapeutic benefit of DNMT inhibitors in this subgroup. In fact, our data from MPN-derived cells lines showed that the cell line with ASXL1 deletion (HEL cell line) was most sensitive to growth inhibition by DNMT inhibitors.

TET2 mutations are also seen frequently in MPNs (28), and recent data have shown that this protein is important for the conversion of methyl cytosine to hydroxymethyl cytosines (29). Hydroxymethyl cytosines can further be removed by base excision repair, thus leading to demethylation. Thus, impairment of this process due to mutated TET2 protein has been associated with decreased levels of hydroxymethylation in experimental models and in primary samples from patients with myelodysplastia (23). Our data now show that decreased hydroxymethylation and increased cytosine methylation can be seen in MPN samples with the TET2 mutation also. The increasing incidence of mutations in proteins involved in the epigenetic machinery (TET, ASXL1, IDH, and others) suggest that the resulting epigenetic alterations may influence gene expression that may contribute to disease pathogenesis. Our demonstration of widespread epigenetic changes in TET2-mutant cases further supports this hypothesis. One might have anticipated that the methylation patterns in neutrophils would be very similar in controls compared with patients with MPNs, due to the intact myeloid differentiation program in the MPNs. Thus, any differences between controls and patients with MPNs or between MPN subclasses may reflect specific genomic lesions unique to the MPN, whether the methylation aberrancies are reflection of the MPN stem cell, progenitor, or differentiated myeloid cell.

Finally, we show that there are numerous significantly and uniformly hypermethylated loci in polycythemia vera, essential thrombocytosis, and PMF that may be targeted by epigenetic modifiers in future clinical trials. Epigenetic modifiers such as histone deacetylase inhibitors are being tried in MPNs including both polycythemia vera and PMF (30). Decitabine has also been tested in PMF with response in 37% of cases (31). These responses do point to the role of epigenetic alterations in disease pathobiology. The heterogeneity in DNA methylation seen by our study in MPNs raises the possibility of specific epigenetic clusters that may be responsive to these agents. Further similar correlative studies are needed to uncover epigenetic signatures of response to these agents.

No potential conflicts of interest were disclosed.

Conception and design: J. Greally, A. Moliterno, A. Verma

Development of methodology: S. Bhattacharyya, D. Sohal, Y. Mo, J.P. Maciejewski, A.M. Melnick, A. Moliterno, A. Verma

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S. Nischal, M. Christopeit, B. Will, A. Pardanani, M. McDevitt, J.P. Maciejewski, A. Moliterno, A. Verma

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S. Nischal, M. Christopeit, Y. Yu, D. Sohal, U. Steidl, A. Verma

Writing, review, and/or revision of the manuscript: M. Christopeit, D. Sohal, A. Pardanani, M. McDevitt, J.P. Maciejewski, A.M. Melnick, U. Steidl, A. Moliterno, A. Verma

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): L. Zhou, T. Bhagat, B. Will, A. Verma

Study supervision: M. Suzuki, J. Greally, U. Steidl, A. Verma

The study was supported by NIH grants R01HL082946 and RO1HL082995, Gabrielle Angel Foundation, Leukemia Lymphoma Society, Hershaft family Foundation, American Cancer Society, German Research Foundation, DFG (CH945/1-1; M. Christopeit), and Immunooncology Training Program T32 CA009173.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.