Targeted deep sequencing in primary myelofibrosis

Mutations in primary myelofibrosis (PMF) are operationally classified into 2 categories: drivers and others. The drivers include Janus kinase 2 (JAK2), calreticulin (CALR), and myeloproliferative leukemia virus oncogene (MPL). Others include additional sex combs such as transcriptional regulator 1 (ASXL1), serine/arginine-rich splicing factor 2 (SRSF2), isocitrate dehydrogenase 1 or 2 (IDH1/2), and enhancer of zeste 2 polycomb repressive complex 2 subunit (EZH2) among others.^1,2  Driver mutations in PMF are often mutually exclusive, and their frequencies are 50% to 60% for JAK2^V617F, 20% to 25% for CALR, and 6% to 7% for MPL^3,4 ; approximately 10% to 15% of patients with PMF do not express any of the 3 driver mutations and are referred to as being triple-negative.^3,5  Recently, whole-exome sequencing has suggested the presence of somatic or germline MPL or JAK2 variants in some of the patients with triple-negative PMF.⁶  Regardless, driver mutations in PMF might be accompanied by other mutations whose pathogenetic relevance is even less clear; they include ASXL1, SRSF2, IDH1/2, EZH2, TET2, DNMT3A, and CBL, with respective frequencies of approximately 22%, 9%, 3%, 5%, 10%, 6%, and 4%.⁷ 

Driver mutations in myeloproliferative neoplasms (MPNs) have now been included in the World Health Organization diagnostic criteria for both prefibrotic and overtly fibrotic PMF.^8,9  Furthermore, the World Health Organization system encourages screening for mutations other than JAK2/CALR/MPL to complement morphologic diagnosis in triple-negative cases.⁸  In terms of phenotypic correlations in PMF, JAK2 mutations are associated with older age, higher hemoglobin level, leukocytosis, and lower platelet count, and mutant CALR is associated with younger age, higher platelet count, and lower frequencies of anemia, leukocytosis, and spliceosome mutations.^10,11  Furthermore, type 2 CALR mutations have been associated with higher risk category, circulating blast percentage, leukocyte count, and inferior survival.¹²  Genotype-phenotype correlations regarding other mutations include clustering of ASXL1 and SRSF2 mutations with older age and high-risk disease and clustering of ASXL1 mutations with leukocytosis, anemia, and constitutional symptoms.⁷  Most recently, the order of mutation acquisition was suggested as an additional determinant of phenotype in MPNs.¹³ 

In terms of disease prognostication in PMF, type 1 or type 1-like CALR mutations have been associated with superior survival,^14,15  and ASXL1/SRSF2/EZH2 mutations have been associated with inferior survival⁷ ; in addition, leukemia-free survival was affected by IDH1/2 mutations.⁷  In a subsequent study,¹⁶  we showed that patients with 2 or more adverse mutations displayed significantly shorter survival compared with patients with 1 adverse mutation who in turn had worse survival than those with no adverse mutations. In this study, we applied a broader mutation screening strategy to identify the occurrence and prognostic value of additional mutations as well as their co-segregation patterns and phenotypic correlations.

This study was approved by the Mayo Clinic Institutional Review Board. Diagnoses and treatment approaches were in accordance with what was considered standard of care at the time of initial diagnosis.¹⁷  The study population was selected on the basis of the availability of sufficient DNA and was not otherwise biased by any other selection criteria. Targeted capture assays were carried out on bone marrow or whole blood DNA for the following genes: TET2 (exons 3, 9, 10, and 11), DNMT3A (exons 4, 8, 13, 15, 16, 18, 19, 20, 22, and 23), IDH1 (exon 4), IDH2 (exon 4), ASXL1 (exon 12), EZH2 (exons 8, 17, and 18), SUZ12 (all exons), SRSF2 (exon 1), SF3B1 (exons 13, 15, and 17), ZRSR2 (all exons), U2AF1 (exons 2 and 6), PTPN11 (all exons), TP53 (exons 5, 6, 7, and 8), SH2B3 (all exons), RUNX1 (exons 3, 4, and 8), CBL (exons 8 and 9), NRAS (exons 2 and 3), JAK2 (exons 12 and 14), CSF3R (exons 14 and 17), FLT3 (exons 14 and 20), KIT (all exons), CALR (all exons), MPL (exon 10), NPM1 (exon 11), CEBPA (exon 1), IKZF1 (all exons), and SETBP1 (exon 3).

Paired-end indexed libraries were prepared from individual patient DNA in the Mayo Genomic Sequencing Core Laboratory by using the NEBNext Ultra Library Prep protocol on the Agilent Bravo liquid handling platform (New England BioLabs; Ipswich, MA; Agilent Technologies Inc., Santa Clara, CA). Capture libraries were assembled according to NimbleGen standard library protocol (Roche NimbleGen, Inc., Basel, Switzerland). A panel including the regions of the aforementioned 27 heme-related genes was selected for custom target capture using Agilent SureSelect Target Enrichment Kit (see below for gene regions covered). Capture libraries were pooled at equimolar concentrations and loaded onto paired-end flow cells at concentrations of 7 to 8 pM to generate cluster densities of 600 000 to 800 000/mm² following Illumina’s standard protocol using the Illumina cBot and HiSeq PE Cluster Kit Version 3 in batches of 48 samples per lane (Illumina, Inc., San Diego, CA). The flow cells were sequenced as 101 × 2 paired-end reads on an Illumina 2000 HiSequencing System using TruSeq SBS Sequencing Kit Version 3 and HiSeq Data Collection Version 2.0.12.0 software. Base-calling was performed by using Illumina’s Real-Time Analysis version 1.17.21.3.

GeneSifter software was used (PerkinElmer, Danvers, MA) to analyze targeted sequence data. Reads from the sequencing in fastq format were aligned by using Burrows-Wheeler Aligner against the genomic reference sequence for Homo sapiens according to Build 37.2 of the National Center for Biotechnology Information (Bethesda, MD). An additional alignment postprocessing set of tools was then used to perform local realignment, duplicate marking, and score recalibration to generate a final genomic aligned set of reads. Nucleotide variants were called by using the Genome Analysis Toolkit (Broad Institute, Cambridge, MA), which identified single nucleotide and small insertion/deletion events using default settings. AlamutVisual mutation analysis software (Interactive Biosoftware, Rouen, France) was used to help filter variations through public genomic databases. Variants were not further analyzed if they did not have a sequencing read depth of ≥50 reads and/or were present with ≤5% variant allele frequency. Minor allele frequency (MAF) was annotated by using both the National Center for Biotechnology Information Single Nucleotide Polymorphism public archive database (dbSNP), and the Exome Aggregation Consortium from the Broad Institute (Cambridge, MA). Remaining variants were further filtered through the Wellcome Trust Sanger Institute Catalogue of Somatic Mutations in Cancer (COSMIC) public database and characterized into 3 categories: VAR: variants not previously associated with a hematologic malignancy (by COSMIC) and either not present or present with ≤1% MAF in dbSNP; HEME: variants previously associated with a hematologic malignancy (by COSMIC or previously by the Mayo Clinic research laboratory; data not shown) and present with ≤1% MAF in dbSNP; and MUT: variants previously associated with a hematologic malignancy and also identified as being somatic (by COSMIC) and present with ≤1% MAF in dbSNP. Each case was assigned a unique patient number and corresponding variants/mutations can be viewed in supplemental Table 1. In cases with multiple variants/mutations in a single gene, the most significant type was illustrated in supplemental Table 1 (ie, MUT>HEME>VAR).

Prognostic evaluation of sequence variants/mutations considered both the number of sequence variants/mutations and the specific genes affected. Adverse sequence variants/mutations were identified by age-adjusted multivariable analysis of their impact on overall or leukemia-free survival. All statistical analyses considered clinical and laboratory parameters obtained at time of first referral at the Mayo Clinic. Differences in the distribution of continuous variables between categories were analyzed by either the Mann-Whitney U test or the Kruskal-Wallis test. Patient groups with nominal variables were compared by using the χ² test. Survival analysis was considered from the date of referral to date of death or last contact. Leukemia-free survival calculations considered the transformation event as the uncensored variable. Survival curves were prepared by the Kaplan-Meier method and compared by the log-rank test. A Cox proportional hazards regression model was used for multivariable analysis. P values less than .05 were considered significant. All analyses were conducted by using StatView (SAS Institute, Cary, NC).

A total of 182 Mayo Clinic patients with PMF were evaluated (median age, 63 years; range, 22-87 years; 65% males). Table 1 outlines the presenting clinical and laboratory details. Driver mutation distribution was 60% JAK2, 22% CALR, 6% MPL, and 12% triple negative. According to the Dynamic International Prognostic Scoring System Plus (DIPSS-Plus) model,¹⁸  survival risk distribution at time of referral was 31% high (n = 57), 38% intermediate-2 (n = 70), 17% intermediate-1 (n = 30), and 14% low (n = 25). Median follow-up was 4 years (range, 0.12-22.2 years). The numbers of documented deaths and leukemic transformations were 112 (62%) and 18 (10%), respectively.

DNA sequence variants/mutations other than JAK2/CALR/MPL were detected in 147 patients (81%), with no difference among driver mutational categories: 83% in JAK2, 73% CALR, 91% MPL, and 82% in triple-negative cases (P = .43) (Tables 1 and 2). Figure 1 depicts all sequence variants/mutations detected and their frequencies whereas supplemental Table 1 provides additional information on specific variants/mutations; the most frequent genes involved were ASXL1 (36%), TET2 (18%), SRSF2 (18%), and U2AF1 (16%). Approximately 35%, 26%, 10%, and 9% of the patients harbored 1, 2, 3, or 4 or more sequence variants/mutations other than JAK2/CALR/MPL, respectively (Table 1). The distribution in the number of mutations was similar among the different driver mutational categories (Table 2; P = .35); however, significant associations were noted for JAK2 with U2AF1 variants/mutations and infrequent occurrence of SRSF2 variants/mutations were noted in CALR-mutated cases (Table 2). Furthermore, among the most frequent sequence variants/mutations, significant associations were apparent for ASXL1 with U2AF1 (P = .002), CBL (P = .007), and SETBP1 (P = .04); for TET2 with CBL (P = .04), SUZ12 (P = .003), and PTPN11 (P < .001); for SRSF2 with IDH1/2 (P < .001) and MPL (P = .047); and for U2AF1 with JAK2 (P = .02), PTPN11 (P = .02), and ASXL1 (P = .002). Conversely, ASXL1 was less likely to coexist with DNMT3A (P = .03) and SRSF2 was less likely to coexist with U2AF1 (P = .006).

In general, older patients were more likely to display U2AF1, SF3B1, TET2, SETBP1, RUNX1, CEBPA, and PTPN11 variants/mutations, whereas female sex was associated with DNMT3A and male sex was associated with CSF3R variants/mutations (P < .05 for all). Additional phenotypic correlations were examined for sequence variants/mutations with recurrence rate of above 10% (Figure 1). Parameters examined for possible association included DIPSS-plus, karyotype, leukocyte count, hemoglobin level, platelet count, and splenomegaly. Anemia was associated with U2AF1 and SRSF2, thrombocytopenia with U2AF1, and constitutional symptoms and higher DIPSS-plus score with ASXL1 and SRSF2 variants/mutations (P < .05 for all, including P = .003 for U2AF1 and anemia association and P = .006 for U2AF1 and thrombocytopenia association). No association with karyotype, spleen size, or leukocyte count was evident.

Age-adjusted univariable analysis identified ASXL1, SRSF2, IDH1/2, EZH2, CBL, KIT, and RUNX1 variants/mutations as being significantly associated with inferior survival; ASXL1, SRSF2, CBL, and KIT remained significant on age-adjusted multivariable analysis, with respective hazard ratios (HRs) of 2.1 (95% CI, 1.4-3.1), 2.0 (95% CI, 1.3-3.0), 3.0 (95% CI, 1.4-6.5), and 38.1 (95% CI, 8.6-168.2), whereas EZH2 and RUNX1 maintained borderline significance. For leukemia-free survival, univariable analysis showed significant associations with SRSF2, RUNX1, CEBPA, SH2B3, and IDH1/2; the first 4 remained significant during multivariable analysis with respective HRs of 4.9 (95% CI, 1.9-12.8), 8.7 (95% CI, 1.8-42.5), 5.4 (95% CI, 1.6-17.6), and 5.8 (95% CI, 1.6-21.7).

On the basis of their independent impact on overall or leukemia-free survival outlined above, ASXL1, SRSF2, CBL, KIT, RUNX1, CEBPA, and SH2B3 variants/mutations were identified as being adverse and at least 1 of them was present in 56% of the 182 study patients. Adverse variants/mutations were associated with inferior overall survival (median, 3.6 vs 8.5 years; P < .001) and leukemia-free survival (7-year risk, 25% vs 4%; P < .001), and the effect on survival was independent of both DIPSS-Plus and JAK2/CALR/MPL mutational status (HR, 2.0; 95% CI, 1.3-3.1). Figure 2 depicts survival data of patients with adverse vs nonadverse vs no detectable variants/mutations; the difference between the latter 2 became insignificant when adjusted for age (P = .15). Additional prognostic information was obtained by considering the number of adverse variants/mutations; median survivals in patients with zero (n = 80), 1 or 2 (n = 93), or 3 or more (n = 9) adverse variants/mutations were 8.5, 4, and 0.7 years, respectively (Figure 3; P < .001), and the difference was independent of DIPSS-Plus with HRs of 7.3 (95% CI, 3.2-16.8) for 3 or more and 1.9 (95% CI, 1.2-3.0) for 1 or 2 adverse variants/mutations.

The frequent detection of multiple DNA sequence variants/mutations in PMF in this study is in line with that of a previous smaller study of 34 patients who underwent a 104-gene panel next-generation sequencing (NGS).¹⁹  The larger number of patients in this study has enabled us to examine the effect of driver mutational status on the particular phenomenon and its relevance to disease phenotype and prognostication. Frequencies of sequence variants/mutations were mostly similar among JAK2/CALR/MPL mutational categories, including triple-negative cases in which more than 80% of the patients displayed at least 1 sequence variant/mutation; however, significant clustering was noted between JAK2 and U2AF1 mutations, whereas SRSF2 mutations were infrequent in CALR-mutated patients (Table 2). In other words, multigene sequencing allows genetic confirmation of clonal myeloproliferation in virtually all patients with PMF, and it also discloses clustering of specific mutations that might signify pathogenetic relevance and provide an explanation for differences in phenotype and prognosis among patients with different driver mutations.

The phenotypic relevance of driver mutations in PMF is now well established and includes the association of JAK2 mutations with older age, higher hemoglobin level, leukocytosis, lower platelet count, and the association of CALR mutations with younger age, higher platelet count, and lower frequencies of anemia, leukocytosis, and spliceosome mutations.¹⁰  In this study, we show that older age was also associated with a number of other variants/mutations, including U2AF1, SF3B1, TET2, SETBP1, RUNX1, CEBPA, and PTPN11; the particular information underscores the need to adjust the analysis of survival impact for age. Additional phenotypic correlations included anemia with U2AF1 and SRSF2, thrombocytopenia with U2AF1, and constitutional symptoms and higher DIPSS-Plus score with ASXL1 and SRSF2 variants/mutations. Most of these associations were previously recognized^7,20,21 ; the phenotypic effect of U2AF1 is particularly noteworthy because of its significant association with both JAK2 and ASXL1 and because of a recent mechanistic report on its potential contribution to abnormal hematopoiesis in patients with myelodysplastic syndromes.²² 

The key objective of this study was to identify variants/mutations with prognostic relevance. The need for such information is dire, considering the increasing use of targeted NGS in myeloid malignancies and its potential to facilitate refinement of current prognostic models. In a previous study⁷  that used mutation analysis of a limited number of genes, we had identified ASXL1, SRSF2, EZH2, and IDH1/2 mutations as risk factors for overall or leukemia-free survival in PMF.⁷  We also showed that the number of such prognostically relevant genes provided additional prognostic value; patients with 2 or more prognostically relevant mutations fared worse, in terms of both overall and leukemia-free survival, compared with those with 1 such mutation, which in turn displayed inferior survival compared with those without such mutations.¹⁶  In this study, we used NGS and a broader panel of genes, which enabled us to add CBL, KIT, RUNX1, SH2B3, and CEBPA variants/mutations to the unfavorable list; whether or not EZH2 and IDH1/2 mutations are thus no longer relevant in the presence of the newly identified adverse mutations requires additional studies for confirmation. This study also confirms the additional prognostic value of a number of adverse variants/mutations but, unlike our previous study,¹⁶  patients with 1 or 2 adverse variants/mutations displayed similar survival, which was inferior to that of patients without adverse variants/mutations and superior to that of those with 3 or more such mutations.

Additional studies are needed to validate and/or refine our observations. For example, some of the variants/mutations we labeled as being adverse have not yet been shown to affect protein function. Similarly, the distinction between rare germline variants and acquired somatic mutations was not fully addressed in this study. Thus, the observations from this study should be viewed with caution and regarded as being preliminary. Regardless, this study provides general information that might help interpret results from multigene sequencing studies that are being implemented with increasing frequency in routine clinical practice. The demonstration of additional (ie, independent) prognostic value of adverse variants/mutations to that of both DIPSS-Plus and driver mutational status suggests the practical value of including such molecular information in currently available prognostic models, which would also allow for a more appropriate use of such molecular information.

This study was supported in part by the Mayo Clinic Harvey-Yulman Charitable Foundation for Myelofibrosis Tissue Bank, by the Clinical Database of Molecular and Biological Abnormalities, and by the Henry J. Predolin Foundation for Research in Leukemia, Mayo Clinic, Rochester, MN.

Contribution: A.T. designed and sponsored the study, contributed patients, collected clinical data and patient samples, reviewed the molecular data, performed statistical analysis, and wrote the paper; T.L.L. performed and analyzed the molecular analysis; C.M.F. performed the molecular analysis; Y.E. collected clinical data; C.A.H. performed pathology review; R.P.K. performed cytogenetics review; N.G. contributed patients and collected clinical data; A.P. contributed patients and reviewed the molecular data; and all authors reviewed and approved the paper.

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

Correspondence: Ayalew Tefferi, Mayo Clinic College of Medicine, 200 First St SW, Rochester, MN 55905; e-mail: tefferi.ayalew@mayo.edu.