Epigenomic profiling of myelofibrosis reveals widespread DNA methylation changes in enhancer elements and ZFP36L1 as a potential tumor suppressor gene that is epigenetically regulated

Abstract

In this study we interrogated the DNA methylome of myelofibrosis patients using high-density DNA methylation arrays. We detected 35,215 differentially methylated CpG, corresponding to 10,253 genes, between myelofibrosis patients and healthy controls. These changes were present both in primary and secondary myelofibrosis, which showed no differences between them. Remarkably, most differentially methylated CpG were located outside gene promoter regions and showed significant association with enhancer regions. This aberrant enhancer hypermethylation was negatively correlated with the expression of 27 genes in the myelofibrosis cohort. Of these, we focused on the ZFP36L1 gene and validated its decreased expression and enhancer DNA hypermethylation in an independent cohort of patients and myeloid cell-lines. In vitro reporter assay and 5’-azacitidine treatment confirmed the functional relevance of hyper-methylation of ZFP36L1 enhancer. Furthermore, in vitro rescue of ZFP36L1 expression had an impact on cell proliferation and induced apoptosis in SET-2 cell line indicating a possible role of ZFP36L1 as a tumor suppressor gene in myelofibrosis. Collectively, we describe the DNA methylation profile of myelofibrosis, identifying extensive changes in enhancer elements and revealing ZFP36L1 as a novel candidate tumor suppressor gene.

Introduction

Philadelphia chromosome-negative myeloproliferative neoplasms, namely polycythemia vera (PV), essential thrombocythemia (ET) and primary myelofibrosis (MF), are characterized by a clonal transformation of hematopoietic progenitors leading to expansion of fully differentiated myeloid cells.1 Primary MF carries the worst prognosis of all Philadelphia chromosome-negative myeloproliferative neoplasms, with progressive marrow fibrosis, extramedullary hematopoiesis, mild to severe splenomegaly and an increased risk of transformation into leukemia.2 Secondary MF can also arise from PV and ET (hereafter referred to as post-PV and post-ET MF, respectively) by mechanisms that are still poorly understood and are clinically and morphologically indistinguishable from primary MF.3

MF has been intensively studied from the genetic perspective;54 in fact, the modified World Health Organization (WHO) diagnostic criteria for Philadelphia chromosome-negative myeloproliferative neoplasms require the demonstration of a genetic marker of clonal hematopoiesis (JAK2, CALR or MPL mutations).6 The frequency of mutations on relevant epigenetic genes (i.e., DNMT3A, EZH2 and ASXL1) suggests that MF might have an epigenetic component that, to our knowledge, remains poorly characterized.5 So far, epigenetic changes such as DNA methylation have been scarcely addressed in MF7 partly due to the limited changes in promoter DNA methylation compared to those in other hematologic malignancies, as previously published by our group.8 DNA methylation of CpG islands (CGI) (mostly on putative promoter regions) has been traditionally studied in both normal and neoplastic hematopoiesis.109 However, high-throughput platforms offer a wider coverage of the genome, allowing a better understanding of DNA methylation dynamics in regions distant from CGI.11 In this regard, enhancer regions have been characterized as potentially relevant sites of DNA methylation outside CGI.1412 Chromatin immunoprecipitation-sequencing studies have enabled reliable mapping of genome-wide active enhancer regions based on histone modifications (e.g., H3K4me1 and H3K27Ac),1615 allowing the identification of enhancers playing a role in dynamic transcriptional regulation during hematopoiesis.17

The present work describes a comprehensive genome-wide analysis of DNA methylation in MF patients, coupled with a gene expression analysis and information on functional chromatin states, compared with those of samples from healthy donors.16 Focusing on potential epigenetic alterations in enhancer regions, we identified ZFP36L1 as a potential tumor suppressor gene with relevance for the pathogenesis of MF.

Methods

Patients’ samples and clinical data

Samples from MF patients (n=39) were bone marrow, granulocytes or total peripheral blood cells. The MF cohort comprised cases of primary MF (n=22), post-ET MF (n=7) and post-PV MF (n=10). Peripheral blood cells from healthy donors (n=6) were used as control samples in this study. All patients were diagnosed using the 2008 version of the WHO classification system of hematologic malignancies.18 Data on JAK2 mutation status were retrospectively available for all patients, whereas no data on CALR and MPL mutations were available. The patients’ data are accessible from the Gene Expression Omnibus (GSE118241).

Samples and patients’ data were provided by the Biobank of the University of Navarra and were processed following standard operating procedures approved by the local Ethics & Scientific Committee. Prior to the collection of samples, all patients consented to the use of their data and to the use of stored material for research purposes.

DNA methylation profiling

DNA methylation was assessed using a Human-Methylation 450K Bead-Chip kit (Illumina, Inc., San Diego, CA, USA) and the data were analyzed by Bioconductor open source software. The analytical pipeline implemented several filters to exclude technical and biological biases and take into account the performance characteristics of Infinium I and Infinium II assays.19 Differentially methylated CpG were defined as previously described.1913 Details on the experimental procedures, annotation of CpG sites, detection of differentially methylated regions, and Gene Ontology analysis20 are described in the Online Supplementary Methods.

Identification of candidate genes targeted by aberrant DNA methylation in enhancers

Data on gene expression profiling from primary MF and healthy peripheral blood samples were obtained from the publicly available Gene Expression Omnibus accession bank number GSE26049.21 Data were further processed using R and the open source Limma package.22 Further details are described in the Online Supplementary Methods.

Luciferase reporter assays

The CpG-free vector (pCPG-L), kindly provided by Dr. Michael Rehli,23 was used to clone the ZFP36L1 enhancer region. Luciferase experiments were performed in triplicate and the details are described in the Online Supplementary Methods. Primer sequences are available in Online Supplementary Table S1.

ZFP36L1 binding motif search

To further validate the potential relevance of the ZFP36L1 gene in MF, the DREME motif discovery algorithm24 was used to assess enrichment of genes with the ZFP36L1 consensus binding sequence among those genes differentially expressed in MF [false discovery rate (FDR)≤ 0.05].

Overexpression of ZFP36L1

A vector containing the ZFP36L1 open reading frame was kindly provided by Dr. Murphy and subcloned into a PL-SIN-GK vector.25 Further details are described in the Online Supplementary Methods.

Statistical analysis

For parametric group comparisons one-way analysis of variance (ANOVA) with the Dunnet correction was used, whereas for nonparametric group comparisons the Kruskall-Wallis test with the Dunn correction was employed. Paired data were analyzed with a Friedman non-parametric test with the Dunn correction for multiple comparisons, for the data with single measurements. Two-way ANOVA with the Tukey correction was used for data with multiple paired measurements. All tests were performed using Prism 7™ software (GraphPad, La Jolla, CA, USA).

Details of other experimental procedures are given in the Online Supplementary Methods.

Results

Myelofibrosis is characterized by a specific DNA methylation pattern enriched in enhancer regions

In order to provide an exhaustive analysis of the DNA methylation profile in patients with MF, we analyzed the DNA methylome of patients with primary MF, secondary MF (including post-ET/post-PV MF) and healthy donors as controls, using the Human- Methylation 450K array. The first result worth highlighting was the epigenetic similarity between primary and post-ET/post-PV MF. Interestingly, with a FDR<0.05, no differentially methylated CpGs were found between primary and secondary MF. Furthermore, we did not identify any differentially methylated CpG between post-ET and post-PV MF. However, both unsupervised principal component analysis (PCA) (Figure 1A) and hierarchical clustering studies (Online Supplementary Figure S1A) using all CpG analyzed confirmed an explicit segregation and a clear epigenetic difference between samples from patients with MF and those from healthy controls. These results allowed us hereafter to consider all MF samples as a single sample cohort.

Figure 1.Patients with myelofibrosis have a different DNA methylation profile from controls, with changes located primarily in enhancer regions. (A) Unsupervised principal component analysis (PCA) showing a differential DNA methylation profile of myelofibrosis (MF) patients and healthy controls with no differences between primary and secondary MF. (B) Distribution of differentially methylated CpG according to CpG island mapping (left graph) or functional chromatin analysis (right graph) grouped by DNA methylation status of the probes (legend). *P≤0.05. (C) Heatmap of DNA methylation levels of differentially methylated CpG sites located in enhancer regions in MF patients and healthy controls. (D) GO-PANTHER analysis of genes adjacent to differenatially methylated CpG located in enhancer regions. Analysis of hypermethylated and hypomethylated genes is shown in the left and right panels, respectively. PC1: principal component 1; PC2: principal component 2; PMF: primary myelofibrosis; PV: polycythemia vera; MF: myelofibrosis; ET: essential thrombocythemia; CGI: CpG islands; DMC: differentially methylated CpG.

Next, we sought to interrogate differences in DNA methylation between MF samples and healthy controls. In this supervised analysis, we detected 35,215 differentially methylated CpG (FDR≤0.05) corresponding to 10,253 coding genes. Among all of these differentially methylated CpG, 65.3% were hypomethylated (corresponding to 22,998 CpG) and the remaining 34.7% were hypermethylated (a total of 12,217 CpG), suggesting that loss of DNA methylation is the predominant alteration in MF. Global DNA hypomethylation has also been a common finding in other hematologic malignancies such as chronic lymphocytic leukemia, multiple myeloma and acute myeloid leukemia.272613

Analysis of the genomic location of differentially methylated CpG showed that both hyper- and hypomethylated CpG were underrepresented in classical CGI and significantly enriched outside CpG islands (Figure 1B). This is an interesting finding, because traditionally, neoplasms acquire hypomethylation outside CGI and hyper-methylation inside the islands,2713 and suggests that patterns of methylation gain in MF might differ from those of other neoplasms. To shed light onto the specific function of the differentially methylated CpG, the chromatin state of each CpG was categorized adapting a publicly available annotation of chromatin immunoprecipitation-sequencing data from CD34 hematopoietic progenitor cells, in which four distinct states were defined: promoter (with H3K4me3), active enhancer (with H3K4me1 and H3K27ac), transcribed regions (showing H3K36me3) and heterochromatin (including H3K9me3 and H3K27me3).16 Both hyper- and hypo-methylated CpG showed significant enrichment in enhancer regions, together with a striking underrepresentation in promoter regions (Figure 1B). Unsupervised clustering of differentially methylated CpG located exclusively in enhancer regions (Online Supplementary Table S2) displayed a clear segregation of the majority of MF patients from healthy controls (Figure 1C) identifying 4,182 hypermethylated and 10,935 hypomethylated probes. These results suggest that patients with MF show an intrinsic aberrant pattern of DNA methylation preferentially located in enhancer regions of the genome.

To further characterize the aberrant DNA methylation of enhancer regions in MF, GO-PANTHER enrichment analysis was performed separately in differentially methylated genes. GO terms with an adjusted FDR<0.05 were selected, showing in the case of hypermethylated enhancers relevant cellular processes such as cellular defense response or induction of apoptosis (Figure 1D).

DNA methylation of enhancer regions is associated with gene expression profile in myelofibrosis

DNA methylation levels of enhancer regions were correlated with the expression of host and adjacent coding genes using publicly available gene expression data of an independent cohort of MF patients and healthy donors (GSE26049).21 Fold increases in gene expression values were grouped according to the hypermethylated (Δβ>0.4) or hypomethylated (Δβ<-0.4) enhancer status in MF versus controls. This analysis showed that enhancer DNA hyper-methylation was associated with decreased gene expression of host/adjacent coding genes. In contrast, hypomethylated enhancer regions were not associated with increased gene expression (Figure 2A).

Figure 2.Aberrant enhancer DNA methylation regulates gene expression in myelofibrosis. (A) Violin density plots of expression of genes with differentially methylated CpG located in enhancer regions. The vertical axis represents log fold change in gene expression. The horizontal width of the plot represents density of data along the y axis. (B) Candidate genes with substantial changes in DNA methylation (FDR <0.01 and Δβ>0.4) and differential gene expression (log_FC<0). Red bars represent the average DNA methylation of all enhancer-mapped probes, the black bars represent the average expression of all probes, and the error bars represent the standard deviation (SD) (C) ZFP36L1 downregulation validation by real-time quantitative polymerase chain reaction analysis of myelofibrosis (MF) patients and three myeloid cell lines (including SET-2) compared to healthy controls (n=3). (D,E) Bisulfite sequencing of the ZFP36L1 enhancer region (D) and promoter region (E) in healthy controls, cell lines and primary MF samples. For each sample, the graph shows the mean ± SD of ten CpG dinucleotides for enhancer regions and 15 CpG dinucleotides for promoter regions. (F) pCpG-L luciferase reporter assay showing the inhibition of luciferase activity after treatment of the ZFP36L1 enhancer region with Sss-I methyltransferase. (G) DNA methylation levels of the enhancer region – the same ten CpG dinucleotides as in (D) after 5-azacytidine treatment of SET-2. (H) ZFP36L1 expression levels after 5-azacytidine treatment of SET-2. Plots/bars indicate mean ± SD. FC: fold change; DMC: differentially methylated CpG; FDR: false discovery rate; CONTROL: healthy controls. MF; myelofibrosis; AZA: 5-azacytidine.

Next, we designed a more stringent approach to identify the set of genes underlying the most significant and substantial changes in enhancer DNA methylation (FDR <0.01, Δβ>0.4), coupled with downregulation of their expression (log_FC<0) (Figure 2B). After identifying a number of potential candidates (27 genes), we focused on ZFP36L1, which codes for a RNA-binding protein that mediates the decay of unstable mRNA with AU rich elements in the 3′ untranslated region.2928 Interestingly, the enhancer region associated with this candidate gene was located in its intragenic region, presumably acting as a cis-regulatory element of ZFP36L1 transcription. It is worth noting that this regulatory element was consistently hypermethylated in the cohort of MF patients and showed the largest number of hypermethylated enhancer-related CpG probes among the final 27-gene list.

ZFP36L1 enhancer hypermethylation correlated with downregulation of expression in MF as compared to controls (Figure 2B and Online Supplementary Figure S1B), which was further confirmed in an independent cohort of MF patients and myeloid cell lines (Figure 2C). Bisulfite sequencing confirmed that DNA methylation of the enhancer region of ZFP36L1 was consistently higher in all MF samples and myeloid cell lines than in control samples, whereas the promoter region remained unmethylated (Figure 2D,E and Online Supplementary Figure S1C). Results obtained from luciferase-reporting assays demonstrated that the exogenous DNA methylation significantly reduced ZFP36L1 enhancer activity (Figure 2F). Moreover, 5′-azacytidine hypomethylating treatment was able to reverse the DNA methylation levels of the enhancer region in vitro, partially restoring the gene expression levels of ZFP36L1 in the SET-2 cell line (Figure 2G,H).

ZFP36L1 acts as a tumor suppressor gene and potentially affects the myelofibrosis transcription profile

We hypothesized that ZFP36L1 downregulation could lead to upregulation of its putative targets in MF. We used DREME, a motif discovery algorithm specifically designed to find short, core DNA-binding motifs enriched in the 3′ untranslated region of genes. We found that the GTATTTDT motif (E-value=4.5×10) was in fact overrepresented in transcripts upregulated in MF patients (Figure 3A). Subsequently, an analysis of motif enrichment was performed, revealing a significant enrichment of upregulated genes in MF patients among the group containing the mentioned motif (P=7.69×10; log_FC >1; P<0.05).

Figure 3.ZFP36L1 rescue decreases cell viability in myelofibrosis. (A) Consensus binding motif for ZFP36L1 obtained by DREME motif discovery among transcripts with putative AU-rich motifs upregulated in myelofibrosis samples. (B) Efficiency of infection measured by the percentage of EGFP-positive cells after lentiviral infection. (C) Quantitative polymerase chain reaction validation of ZFP36L1 restoration in the SET-2 cell line after lentiviral infection. (D) ZFP36L1 protein restoration measured by western blot in the SET-2 cell line after lentiviral infection. (E,F) ZFP36L1 rescue with lentiviral vector infection in the SET-2 cell line decreased cell proliferation rate (E) and increased annexin V-positive cells (F).

To complement DREME analysis, we searched the AREsite30 database for AU-rich elements to determine whether we could detect, among the genes differentially expressed (B-value>10) between MF and controls, an enrichment of these sequences in the upregulated subset. Of all the possible AU motifs, we focused on the most restricted 9, 11 and 13-mer motifs. Interestingly, we were able to identify an enrichment of a 9-mer sequence WTATTTATW (P=0.01) and a 13-mer sequence WWW-TATTTATWWW (P=0.03) exclusively among the upregulated genes in MF patients (Figure 3A). Remarkably, both AU motifs strongly resemble the ZFP36L1 core-binding motif predicted by the DREME algorithm.

Re-expression of ZPF36L1 was achieved through lentiviral infection of the SET-2 cell line. Seventy-two hours after infection, the levels of EGFP-positive cells used as the positive control confirmed successful infection, and the level of expression of ZFP36L1 confirmed satisfactory overexpression of the gene (Figure 3B-D). Rescue of ZFP36L1 expression resulted in a decrease of more than 50% in cell proliferation, alongside an increase of annexin V-positive cells as measured by flow cell cytometry (Figure 3E,F).

Discussion

In the present study, we have extended previous knowledge regarding the DNA methylome in both primary and secondary MF, focusing specially on those CpG sites located in enhancer regions of the genome. A preliminary analysis of the global DNA methylome revealed the absence of DNA methylation differences between primary and secondary MF. This constitutes the first key finding of the present study and allowed us to use all the MF samples in a single cohort for further analysis. Primary and secondary MF are known to have very similar biological features, presenting symptoms and clinical course and in fact, both entities are treated indistinctively according to most published guidelines313 Nevertheless, some recent evidence from large retrospective trials has suggested that traditional prognostic factors may not be applicable to secondary MF as patients with post-ET MF seem to survive longer than those with post-PV MF or primary MF.33313

The remarkably homogenous epigenetic profile of all our MF samples supports a common biological origin of primary and secondary MF.34 The DNA methylomes of the novel MF subtypes defined by the new 2016 WHO classification (prefibrotic and overt MF) remain to be characterized and it will be interesting to establish whether these subtypes have different methylation profiles. This aspect exceeded the possibilities of our cohort (retrospective availability of histology samples) but warrants further investigation.

Although previous studies have already interrogated the DNA methylation landscape of MF,7 their findings are limited to small numbers of epigenetic abnormalities mainly focused on the study of promoter regions. Our genome-wide approach of DNA methylation analysis using the 450k array allowed us to interrogate regulatory regions outside traditional promoters and obtain a deeper insight into the aberrant DNA methylome of MF. Changes in DNA methylation levels are known to cooperate with the deposition of chromatin marks, particularly H3K4 methylation, to render the enhancers/promoters accessible/inaccessible to the transcription machinery.3735 Hence, the changes in DNA methylation observed in MF are expected to have an impact on the transcriptional profile of MF and potentially contribute to the MF malignant phenotype. Enhancer DNA methylation changes have been described to play a more prominent role in transcriptional regulation than promoter DNA methylation, governing processes such as hematopoietic differentiation and neoplastic transformation through the regulation of key transcription factors and genes.3837271312 Translated into the context of Philadelphia chromosome-negative myeloproliferative neoplasms, this evidence might support the involvement of aberrant enhancer DNA methylation in the abnormal pattern of differentiation leading to MF. Enhancer hyper-methylation has been reported in neutrophils,12 B cells,39 AML cells26 and myeloma13 adding evidence to dynamic enhancer DNA hypermethylation as a relevant regulatory mechanism of gene expression both in normal and neoplastic hematopoietic cells.

Although the potential involvement of ZFP36L1 in myeloid differentiation has been described previously,40 our results suggest that epigenetic downregulation of ZFP36L1 might be a prominent event in the pathobiology of MF; more importantly, hypermethylation of an enhancer regulatory element represents a novel mechanism of disrupted gene expression in the context of MF and ZFP36L1. ZFP36L1 has been previously implicated in normal hematopoiesis41 and specifically associated with erythroid and myeloid differentiation,40 suggesting a possible role of this gene in MF onset and progression. Moreover, ZFP36L1 is also known to mediate mRNA decay of genes relevant to cell proliferation, survival and differentiation such as CDK6, TNFα, BCL2, NOTCH1 and STAT5B.4342 Interestingly, the enhancer region associated with this candidate gene was consistently hypermethylated in the cohort of MF patients and was located in its intragenic region, presumably acting as a cis-regulatory element of ZFP36L1 transcription. The motif discovery experiments support our hypothesis of epigenetic deregulation of ZFP36L1, suggesting that MF samples with ZFP36L1 loss of expression experience upregulation of the gene’s putative targets. Consequently, when ZFP36L1 expression levels are restored with the lentiviral model, SET-2 cells lose their malignant proliferative phenotype, strengthening the tumor suppressor role of this gene in MF. Taken together, these results link ZFP36L1 to the pathobiology of MF, ultimately resulting in transcriptome deregulation of genes relevant to cell proliferation, survival and differentiation, as previously described.45444240

Conclusion

The DNA methylation landscape of patients with primary MF or post-ET/post-PV MF is consistently different from that of healthy individuals. The absence of differences between primary MF and post-ET/post-PV MF suggests that the changes seen in MF are founding epigenetic alterations occurring at the level of stem cells of this myeloproliferative neoplasm and maintained in differentiated myeloid cells. Aberrant DNA methylation in MF is predominantly located in enhancer regions and has a significant impact on the expression of their target genes. Combining DNA methylation and gene expression data, we identified ZFP36L1 as an attractive new possible therapeutic target that shows a decrease of gene expression mediated by enhancer hypermethylation. Our results also suggest a direct effect of ZFP36L1 downregulation on the gene expression profile of MF, through upregulation of mRNA harboring canonical sites with AU-rich elements. In vitro rescue of ZFP36L1 expression had an impact on cell proliferation and induced apoptosis in the SET-2 cell line, indicating a possible role of ZFP36L1 as a tumor suppressor gene in MF. Moreover, treatment with 5’-azacytidine further evidenced the plausibility of ZFP36L1 pharmacological manipulation. Taken together, these results provide evidence of an unexplored therapeutic target in MF patients, which remains to be properly evaluated in the pre-clinical setting.

Acknowledgments

We particularly acknowledge the patients for their participation and the Biobank of the University of Navarra for its collaboration. We thank John J Murphy and Amor Alcaraz for providing ZFP36L1 expression constructs. This research was funded by grants from Instituto de Salud Carlos III (ISCIII) PI14/01867, PI16/02024 and PI17/00701, TRASCAN (EPICA), CIBERONC (CB16/12/00489; co-financed with FEDER funds), RTICC (RD12/0036/0068) and the Departamento de Salud del Gobierno de Navarra 40/2016. NM is supported by a FEHH-Celgene research grant, MP was supported by a Sara Borrell fellowship CD12/00540 and RO was supported by the Ministerio de Ciencia, Innovación y Universidades of Spain, Subprograma de Formación de Profesorado Universitario (FPU) award number FPU14/04331.

Footnotes

• ↵# These authors share first authorship
• * These authors share senior authorship
• Check the online version for the most updated information on this article, online supplements, and information on authorship & disclosures: www.haematologica.org/content/104/8/1572

• Received August 23, 2018.
• Accepted January 15, 2019.