Patients with triple-negative, JAK2V617F- and CALR-mutated essential thrombocythemia share a unique gene expression signature

Classical Philadelphia-negative myeloproliferative neoplasms (MPNs) are multipotent, clonal hematopoietic stem cell neoplasms encompassing essential thrombocythemia (ET), polycythemia vera (PV), and primary myelofibrosis. These disorders have an overlapping yet heterogeneous clinical spectrum, with both ET and PV having the potential to transform to myelofibrosis.^1,2  Furthermore, MPNs predispose to an increased risk of thromboembolic events and have an inherent risk of transformation to acute myeloid leukemia (AML; 1% to 2% for ET/PV and up to 10% to 20% in primary myelofibrosis at 10 years).^2,3 

ET is specifically characterized by increased numbers of abnormal, clonal megakaryocytes and a clinical phenotype of thrombocytosis. Although rare cases of hereditary ET have been observed, the majority of ET cases occur sporadically, most commonly associated with somatic mutations in either JAK2 exon 14, CALR exon 9, or MPL exon 10, all of which lead to upregulated JAK-STAT signaling.⁴  Mutations in JAK2, CALR, and MPL are usually mutually exclusive, with an incidence of 55%, 25%, and 3%, respectively.⁵  However, 10% to 15% of patients with ET do not harbor detectable “driver” mutations, at least according to conventional diagnostic assays, and are termed “triple negative” (TN). Previous research using whole-exome sequencing (WES) has reported novel MPL and JAK2 mutations in TN MPN, highlighting disease heterogeneity.⁴ 

Limited data exist concerning both extended gene expression analyses and detailed methylation profiles of ET progenitor cells, particularly regarding the TN ET population. Using genome-wide methylation arrays, Pérez et al⁶  evaluated the epigenetic landscape of 71 patients with chronic MPN and 13 with transformed MPN compared with healthy volunteers (HVs). Although a distinct MPN-associated methylation signature was evident compared with HVs, no differential pattern existed between the chronic disease categories (primary myelofibrosis vs PV/ET). The differentially methylated regions appear enriched in a gene network centered on inflammation and the NFκB pathway. In a cohort of patients with MPN, Rice et al⁷  described aberrations affecting components of a noncanonical polycomb repressive complex 2 (EZH1, SUZ12, and JARID2) and, interestingly, upregulation of the transcription factor (TF) NFIB, resulting in cellular protection from apoptosis. More recently, Wong et al⁸  reported significant pro-inflammatory gene upregulation in an MPN cohort.

We hypothesized that patients with ET have a unique gene expression and methylation signature that may serve as a diagnostic tool in those who lack classical driver mutations. We therefore evaluated 46 TN ET patients, aiming to delineate if there are differential epigenetic and gene expression profiles compared with ET patients with known driver mutations.

A total of 56 patients with ET, 2 patients with post-ET myelofibrosis as defined by the 2008 World Health Organization criteria, and 4 healthy HVs were enrolled into the study. Baseline mutation status to exclude the presence of JAK2V617F, CALR, or MPL mutations was determined per local clinical laboratory protocols. Written informed consent was obtained in accordance with the Declaration of Helsinki and with approval from the local ethical review committee.

Using peripheral blood mononuclear cells (PBMCs), WES libraries were prepared using the SureSelect^XT Human All Exon V4 (Agilent Technologies) target enrichment system as per manufacturer’s protocols. Details of the bioinformatics analysis pipeline are provided in the supplemental Methods. Briefly, genomic DNA libraries prepared by using SureSelect^XT Human All Exon V4 were sequenced on the Illumina HiSeq 2000 platform (Illumina) with 100 bp paired-end reads. Reads were aligned to the reference human genome (hg19, b37) using NovoAlign (Novocraft Technologies). Alignments underwent quality control and filtering. In each sample, >87% of the GENCODE-defined coding bases of the exome was represented by at least 20 reads. Pairwise variant calling was performed by using SAM tools and filtered by using VarScan2.

Targeted deep sequencing of genes known to be mutated in myeloid malignancies was performed by using the HaloPlex^HS platform (Agilent Technologies) on peripheral blood PBMCs. The mean sequencing coverage was 95%, and 90% of target sequences could be analyzed with a coverage >90%. A list of the genes covered and details of the analysis pipeline are provided in supplemental Table 1.

RNA-sequencing was performed on a selected subset of individuals based on sample availability, which included 12 “true” TN patients, 3 known JAK2V617F and 3 CALR mutated patients, and 4 HVs. RNA was extracted from PBMCs by using TRIzol reagent (Life Technologies). Libraries were prepared by using the Illumina TruSeq Stranded RNA kit following the manufacturer’s instructions and sequenced by using an Illumina HiSeq 2500 instrument. Differentially expressed genes were identified by using limma.⁹  Full details of the bioinformatics analysis are provided in the supplemental Methods.

Methylation profiles were performed in a selected cohort of enrolled individuals based on sample availability: 21 TN patients (12 of 21 had gene expression performed as described earlier) and 4 JAK2V617F mutated, 2 MPL mutated, and 3 CALR mutated (one JAK2V617F and 2 MPL patients were identified in retrospect from the WES data) patients. Genomic DNA extracted from PBMCs was bisulfite converted by using the EZ DNA Methylation kit (Zymo Research), as per the manufacturer’s protocol and processed on an Infinium Human-Methylation 850K BeadChip platform. Further details are provided in the supplemental Methods.

To investigate the abundance of the gene-specific transcript of 20 candidate genes identified with altered expression, we prepared complementary DNA from total RNA in an additional 20 peripheral blood samples of ET patients (9 JAK2V617F positive, 7 CALR positive, and 4 TN) and 6 HVs¹⁰  (supplemental Table 2). Gene expression was assessed with reverse-transcriptase quantitative polymerase chain reaction (RT-qPCR) using gene-specific probes for the following 20 candidate genes: CEBPβ, ELK1, FAXDC2, FOSL2, GP1BB, GTPBP1, ITGA2B, ITGA7, ITGB3, JUN, JUND, KLF1, MAPKAPK2, MARCKS, NFκB2, PFKFB3, RELA, SPI1, TAL1, and TGFB1 (supplemental Table 3). Assays were run in triplicate on the ABI7900 platform (Applied Biosystems).

Assays were performed as previously described¹¹ ; ABL1 cycle threshold values ≤30 were used for quality assessment. We quantified differences in the expression level (fold change) of the 20 candidate genes between the ET cohort and HVs by using a comparative cycle threshold values platform (Applied Biosystems), and gene expression in the ET samples was compared with HVs, with ABL1 as a reference gene. No-template and no-reverse-transcriptase controls were used to ensure specific amplification of target transcripts.

The clinical and hematologic parameters and treatment regimens for the TN patients enrolled in this study are shown in Table 1.

Targeted, error-corrected sequencing and WES identified somatic or germline mutations in 10 patients (Table 2). Six patients were identified as having mutations in MPL: 4 patients with mutations in exon 10 (W515R, W5151G/S, and S505N), 1 patient with a mutation in exon 11 (R537W), and 1 patient with a mutation in exon 9 (MPL P453R). Four patients had the JAK2V617F mutation. Among the 6 patients identified as having an MPL mutation, this was a germline mutation in 3 patients, including an MPL P435R mutation in a female patient with post-ET myelofibrosis, and an MPL S505N mutation shared by a father and son pair. In the latter cases, the MPL S505N mutation was confirmed to be a germline mutation as it was detected by WES in both PBMCs and skin.

In total, 22% of clinically diagnosed TN ET patients in our cohort harbored driver mutations detected by next-generation sequencing that were not detected by conventional testing. However, most patients displaying a clinical ET phenotype still had no detectable driver mutations via next-generation sequencing, and hence we proceeded to investigate the gene expression profiles of these so-called “true” TN patients.

RNA-sequencing was performed on 12 true TN patients (ie, no driver mutations identified by next-generation sequencing), three JAK2V617F-mutated patients, three CALR-mutated patients, and four HVs. First, the data were interrogated to identify if any gene expression differences existed between the TN cohort and patients with known CALR and JAK2V617F mutations. Notably, gene expression profiles correlated strongly with each other across all ET patients regardless of mutational status (supplemental Figure 1A). Moreover, direct comparisons between the 12 TN patients and 6 patients with known mutations yielded limited significant differences (supplemental Figure 1B). In contrast, compared with that of HVs, the gene expression profiles of each ET subset showed consistent significant differences, regardless of disease subtype, platelet counts, treatment, and age (Figure 1A). Because there seemed to be a set of unique genes in the TN comparison, we further analyzed these genes; however, no significant differences in expression were determined comparing TN samples vs the rest of the ET cohort. Moreover, these genes followed a similar expression pattern in the JAK2V617F and CALR samples compared with the HV samples (supplemental Figure 1C-D).

From our adjusted analysis comparing all ET patients vs HVs, a total of 838 differentially expressed genes (DEGs; adjusted P <.05) (Figure 1B) were identified. Clearly, these genes define a gene set that could distinguish ET from HV in this cohort (Figure 1C). We also recognized that the disease-modifying nature of interferon treatment may affect gene expression patterns; however, sample clustering highlighted no anomalies (supplemental Figure 1A,H). We performed an additional differential analysis comparing patients with and without interferon and could not show any significant differences in gene expression profiles (supplemental Figure 1I-J).

Ingenuity pathway analysis (IPA) of the 838 DEGs highlighted upregulated genes were enriched for genes critical for myeloid cell differentiation and proliferation, specifically activation of megakaryopoiesis and oncogenic signatures (Figure 1D; supplemental Figure 2A). The converse was also true (ie, enrichment for thrombocytopenia among downregulated genes). Notable upregulated DEGs included ITGB3, ITGA2B, ITGA7, GP1BB, FAXDC2, SPI1 (PU.1), TAL1, TGFB1, and CSF2RB; their mean expression levels are displayed in Figure 1E.

To understand the potential molecular pathways that may be driving the expression of these proliferative genes critical for megakaryopoiesis, functional enrichment was performed using the KEGG database. This analysis revealed that upregulated DEGs were enriched for several proliferative and inflammatory pathways, including JAK-STAT, tumor necrosis factor (TNF), NFκB, and MAPK pathways (P < .001) (supplemental Figure 2B-C). Regarding proliferative pathways, key downstream targets of these pathways, including the TFs JUN, JUND, and CCAAT enhancer binding protein β (CEBPβ), were also upregulated.

These data suggest that regardless of disease subtype, ET converges onto similar gene expression profiles, which potentially drive the disease phenotype. To test this theory in an independent cohort, we processed previously published gene expression array data (described in the supplemental Methods) from a mixed ET cohort (n = 20) with HVs (n = 20).¹²  Directly comparing the ET samples vs the HV samples, we identified that 32 of the upregulated genes delineated in our analysis were also significantly upregulated in this cohort. Among the overlapping genes were key inflammatory genes: FAXDC2, ITGA2B, ITGA7, FOSL2, and ITGA2B (supplemental Figure 2D).¹³ 

To gain further validation of the gene expression profile differences observed between the ET patients and HVs, RT-qPCR was performed on an additional 20 independent patients with ET (9 JAK2V617F positive, 7 CALR positive, and 4 TN) and six HVs (supplemental Table 1) on 20 candidate genes identified to have altered expression. The cycle threshold values for the 20 candidate genes (CEBPβ, ELK1, FAXDC2, FOSL2, GP1BB, GTPBP1, ITGA2B, ITGA7, ITGB3, JUN, JUND, KLF1, MAPKAPK2, MARCKS, NFκB2, PFKFB3, RELA, SPI1, TAL1, and TGFB1) are listed in supplemental Table 3. As shown in Figure 1F, all of the candidate genes had significantly higher expression in the patients with ET compared with the HVs, irrespective of mutation subtype.

We next explored whether the epigenetic landscape at the level of DNA methylation profiles differed among the patients with ET according to their mutation subtypes. The DNA methylation profiles of 25 TN ET patients, 4 JAK2V617F-mutated ET patients, and 3 CALR-mutated ET patients were determined by using the human methylation 850K array. Similar to their gene expression profiles, all ET samples, regardless of mutation subtype, shared highly correlated methylation profiles, with a median Pearson correlation of 0.98 (supplemental Figure 3A). Methylation profiles did not cluster on the basis of cytoreductive agent, age, sex, or white blood cell composition (supplemental Figure 3B). On average, 43% of cytosine guanine dinucleotides (CpGs) were hypermethylated and 25% were hypomethylated. As expected, gene promoters were predominantly hypomethylated, whereas gene bodies showed a methylation distribution similar to that of total CpGs (supplemental Figure 3C).¹⁴ 

For each patient sample, we aggregated promoter- or gene body–associated CpGs for each gene and explored regions that had consistently hypomethylated/hypermethylated CpGs across all ET patient samples for that gene. We found that 9611 genes had hypomethylated promoters and 71 genes had hypermethylated promoters; 33 genes had hypomethylated gene bodies and 159 genes had hypermethylated gene bodies (Figure 2A). As expected, correlation with gene expression profiling showed that genes with hypermethylated promoters had significantly lower expression than genes with either hypomethylated promotors or hypomethylated/hypermethylated gene bodies (Figure 2B; supplemental Figure 3D). Interestingly, however, genes identified as having a consistent hypomethylated/hypermethylated pattern at either their promoter of gene body had a significantly higher expression than genes without a clear hypomethylation/hypermethylation pattern. This was consistent for all patient samples in the RNA-sequencing cohort. Although the exact significance of gene body methylation remains unclear, gene body hypermethylation has been reported to be associated with increased gene expression.^15,16 

We next questioned if there were any methylation changes compared with HV PBMCs (n = 2). We identified 652 CpGs with significantly different methylation (adjusted P ≤ .01), of which 572 (80%) CpGs had increased methylation and 80 (20%) had decreased methylation in the ET cohort (supplemental Figure 4A). There was no significant overlap between genes with differential methylation and DEGs (χ² > 0.5), with only 15 overlapping genes. This finding suggests that gene expression changes between ET and HV result primarily from DNA methylation–independent mechanisms.

We subsequently investigated the promoter methylation patterns of DEGs. A significant difference was observed in the median methylation between upregulated and downregulated gene promoters (Figure 2C). Interestingly, however, both upregulated and downregulated promoters were predominantly hypomethylated. Overall, 58% of the DEGs had consistently hypomethylated promoters, and almost none were hypermethylated (Figure 2D). Furthermore, the 58% of DEGs with hypomethylated promoters maintained enrichment for key genes involved in upregulated pathways uncovered in our RNA-sequencing analysis, including the TNF, MAPK, and platelet pathways. Overall, although promoter hypomethylation is consistently observed at upregulated DEGs, it does not appear that changes in DNA methylation alone can account for these changes in gene expression.

A fundamental aspect of gene regulation is by TFs and repressors. Thus, to explore the regulatory influence of TFs on the DEGs, the promoters of all upregulated DEGs were assessed for TF-binding motifs. There was a significant enrichment for many TFs involved in myeloid cell development, including RARA, PU.1, and several KLF genes (Figure 3A). Twenty-five of the enriched TFs were highly expressed in the ET cohort, with expression levels in the top quartile. A significant enrichment for TFs key in the NFκB, MAPK, and TNF pathways such as JUN, JUNB, FOSL2, RELA, and BCL6 were identified, all of which were upregulated in the ET cohort. Interestingly we also observed enrichment for RARA and its cofactor RXRA. Using this approach, several TFs were identified that could be driving the expression of genes leading to the ET phenotype.

To understand if these gene promoters are active and accessible to allow TF binding in an ET context, we interrogated publicly available RNA-sequencing and H3k27ac (acetylation) data in the SET-2 cell line.^17,18  We first explored expression of the upregulated genes previously uncovered. Most upregulated genes have a high expression in the SET-2 cell line, with 63% in the upper expression quartile, which is significantly higher than average (Figure 3B). We also found that 42% of gene promotors have significant H3K27ac deposition, including ITGA2B, ITGB3, FOSL2, TAL1, JUND, and MAPK3 (supplemental Figure 5). Those genes with this active mark had a considerably higher expression overall. This suggests these promoters are highly active in this cell line and could be regulated through the TF-binding sites we identified.

In keeping with the literature, WES and targeted sequencing facilitated the identification of canonical and noncanonical mutations in MPL and JAK2 in 22% (10 of 46) of our TN cohort. This included 2 of 10 patients harboring a somatic MPL mutation in exon 10 at positions W515R/G/S, all previously reported in patients with MPN.¹⁹  One patient had a mutation on exon 11 position R537W. The MPLR537W mutation is not widely described in ET; this missense mutation in MPL located within the cytoplasmic domain results in failure of JAK binding and signal transduction, and has been associated with congenital amegakaryocytic thrombocytopenia.²⁰  We found MPL germline mutations in 3 patients: MPL P453R (n = 1) and MPL S505N (n = 2). The mutation on exon 9, position P453R is associated with the cytokine receptor domain and has not been previously detailed in ET. The germline mutation in S505N, a point mutation that leads to the loss of the inhibitory effect on the TPO genes, was identified in 2 patients in the study cohort and is well recognized in several kindreds.²¹  Four of 10 patients harbored the classic mutation JAK2V617F. No germline JAK2 mutations were identified in our cohort.

It is important to note that the absence of mutations does not exclude alternative disease drivers, as it is well established that WES cannot detect fusion oncogenes, mutations in regulatory sequences, or alterations in gene expression. To further characterize the remaining true TN cohort, we used RNA-sequencing to directly compare them with ET patients with known driver mutations and HVs. Significant correlations in gene expression were noted regardless of mutation status, and consistent differences were found between the ET patients and HVs. Consistent with the well-established “inflammatory” phenotype associated with MPN, we found upregulation of several pathways involved in hematopoiesis, inflammation, and immune regulation, specifically NFκB, MAPK, and TNF signaling, in parallel with upregulation of genes that drive megakaryopoesis.²² 

We identified 838 genes in ET patients that are consistently differentially expressed compared with HVs. Contrary to our expectation, however, ET gene expression profiles correlated highly with one another and did not cluster according to mutation subtypes. Indeed, this was also observed at an epigenetic level in terms of DNA methylation profiles, with no clustering based on mutation subtypes. This was validated further with RT-qPCR in an independent ET cohort, which exhibited significant upregulation of 20 key candidate genes in the ET population compared with HVs, irrespective of mutation subtype. Taken together, our systematic genetic and epigenetic profiling indicates that different driver mutations may activate similar molecular pathways that converge on a similar gene expression profile that ultimately results in the ET phenotype.

In our cohort, upregulation of genes associated with megakaryopoiesis, including ITGA2B and ITGB3, were identified. Such upregulation of ITGA2B and ITGB3 is in keeping with previous findings by Müller et al,²³  who also described altered expression in an MPN population. The ITGA2B gene encodes the αIIb subunit (GPIIb) and is restricted to the megakaryocyte lineage, whereas the ITGB3 gene encodes β3 (GPIIIa) and occurs in many cell types as a component of the vitronectin receptor (αvβ3).^24,25  Together, they form the GPIIb/IIIa complex, which is a platelet–platelet contact receptor for the clotting factors von Willebrand factor and fibrinogen in activated platelets, and initiates aggregation.^25,26  The synthesis of GPIIb/IIIa occurs primarily in the megakaryocytes and is detectable at an early stage of maturation.²⁷ 

As highlighted, upregulated DEGs displayed consistent promoter hypomethylation. The DEG promoter sequences were enriched for TF-binding motifs, including that of JUN family members. We identified that a large proportion of these gene promoters may be highly active and accessible to the TFs in an ET cell line (SET-2). The activating protein-1 TF family is made up of multiple Jun (cJun, JunB, and JunD) and FOS members. Activating protein-1, as with NFκB, regulates many of the inflammatory and immune pathways, functioning downstream of the MAPK pathway and resulting in essential cellular processes that include proliferation and differentiation.²⁸ 

The JAK signal transducer, which has a central role in MPN pathogenesis, is critical for signaling of multiple cytokines and, in addition to STAT, leads to MAPK pathway activation with resultant dysregulation of cellular proliferation.²⁹  Across this ET cohort, we reported the upregulation of genes downstream of the MAPK pathway that included MAP3K11, ELK1, NFκB, JUND, and JUN (supplemental Figure 2). We also observed an upregulation of the FAXDC2 gene, which has been recently implicated as a driver critical for megakaryopoiesis, acting through the activation of MAPK/extracellular signal-regulated kinase signaling.³⁰  Moreover, we identified upregulation of NFκB, in parallel with the subunit RELA, which is the primary driver for inflammation and cell proliferation.³¹  In many cells, NFκB promotes a prosurvival state by enhancing the G1 phase to the S phase cell cycle progression and upregulating antiapoptotic proteins such as BCL2 and BCL-xl; however, in platelets, their functional relevance remains unclear.^32-35 

To our knowledge, this study is the first to use a systematic approach to TN ET patients in the evaluation of molecular, gene expression, and methylation profiles. We acknowledge that a potential weakness of our work is the absence of comparison with reactive thrombocytosis patients; this topic could be assessed in future studies. However, we are able to show that a proportion of TN ET patients have low-level canonical and noncanonical mutations in MPL and JAK2 genes. ET gene expression and DNA methylation profiles correlated highly with one another and did not cluster according to mutation subtypes. Thus, it is suggested that different driver mutations activate similar molecular pathways that converge on a similar gene expression profile (that includes upregulation of genes implicated in platelet proliferation such as ITGA2B and ITGB3), which ultimately results in the development of ET. We speculate that this gene expression profile may be regulated by JUN family TFs, downstream of NFκB, MAPK, and TNF signaling. We hence provide an insight into the transcriptional and epigenetic regulation of ET and identify putative underlying molecular pathways that should be investigated further for the potential for therapeutic targeting.

The authors thank all the patients and HVs who participated in this study. The authors extend a special thanks to the Cangen group for their support throughout the project.

Contribution: S.A. designed the study, collected and processed samples, conducted experiments, analyzed data, and wrote the manuscript; W.V. analyzed data and contributed to writing of the manuscript; R.D. oversaw experimental design, processed samples, and contributed to the manuscript; M.S., A.S., P.L., and K.M. performed data analysis; P.C. performed data analysis and contributed to the manuscript; M.R., A.K., J.K.-C., S.O., and B.B.P. processed samples; D.R. contributed patients to the study; C.W. and Y.F. performed sample collection; S.K. contributed to the manuscript; G.M. supervised the project; D.P.M. and C.N.H. supervised the project and contributed to writing the manuscript; and all authors read and approved the submitted manuscript.

Conflict-of-interest disclosure: C.N.H. reports funded research from Novartis; speaker fees from Novartis, Janssen, CTI, Celgene, and Medscape; and advisory board membership for Incyte, CTI, Sierra Oncology, Novartis, Celgene, Roche, AOP pharma, Geron, and AstraZeneca. D.P.M. reports speaker fees and advisory board memberships with Novartis, Celgene, and Jazz Pharmaceuticals. S.K. reports research grants from Celgene and Novartis; and speaker honorarium from Alexion. The remaining authors declare no competing financial interests.

Correspondence: Samah Alimam, Department of Haematology, 4th Floor, Southwark Wing, Guy’s and St. Thomas’ NHS Foundation Trust, Great Maze Pond, London SE1 9RT, United Kingdom; e-mail: samah.alimam@nhs.net; and Claire N. Harrison, Department of Haematology, 4th Floor, Southwark Wing, Guy’s and St. Thomas’ NHS Foundation Trust, Great Maze Pond, London SE1 9RT, United Kingdom; e-mail: claire.harrison@gstt.nhs.uk.