Integrated molecular profiling of juvenile myelomonocytic leukemia

Juvenile myelomonocytic leukemia (JMML) is a rare myelodysplastic/myeloproliferative neoplasm that occurs in infants and during early childhood, characterized by excessive myelomonocytic cell proliferation and granulocyte-macrophage colony-stimulating factor hypersensitivity.¹  More than 80% of JMML patients harbor mutually exclusive somatic and/or germ line mutations in canonical RAS pathway genes, such as PTPN11, NF1, NRAS, KRAS, and CBL.^2-4  However, the molecular pathogenesis of JMML in 10% to 20% of patients remains unknown. Secondary genetic events such as mutations in SETBP1^2,5  and other genes,^3,4  gene expression profiles,^6,7  and selected gene promoter methylation profiles^8,9  correlate significantly with clinical outcomes. However, the interrelationships among these biomarkers and genome-wide methylation status of JMML are poorly understood. To address these issues, we conducted an integrated analysis of genome-wide methylation, resequencing, and RNA-sequencing (RNA-seq) to detect the molecular pathogenesis of JMML in patients without RAS pathway mutations and to investigate the selective co-occurrence of clinical and biological prognostic markers.

We studied 150 children (98 boys and 52 girls) with JMML, including 15 with Noonan syndrome-associated myeloproliferative disorder (NS/MPD), who were diagnosed with JMML in institutions throughout Japan. Of these, 92 were included in our previous publication.²  Written informed consent was obtained from the parents of the patients before sample collection. This study was approved by the ethics committee of the Nagoya University Graduate School of Medicine and was conducted in accordance with the principals of the Declaration of Helsinki. Diagnosis of JMML was made on the basis of internationally accepted criteria.¹⁰  Characteristics of the 150 patients are summarized in supplemental Table 1 available on the Blood Web site. The median age at diagnosis was 15 months (range, 1-160 months). Karyotypic abnormalities were detected in 28 patients, including 14 with monosomy 7. Of the 150 patients, 86 (57%) received allogeneic hematopoietic stem cell transplantation.

Genomic DNA and RNA were extracted from mononuclear cells derived from peripheral blood or bone marrow, CD3⁺ T cells, and umbilical cord tissue by using a QIAamp DNA Blood Mini Kit, a QIAamp DNA Investigator Kit, and an RNeasy Mini Kit (QIAGEN, Hilden, Germany), according to the manufacturer’s instructions. Human T-Cell Activation/Expansion Kits (Miltenyi Biotec, Bergisch Gladbach, Germany) were used for the expansion of CD3⁺ T cells from patients’ peripheral blood or bone marrow mononuclear cells.

We performed whole-exome sequencing (WES) of paired tumor-normal samples from 69 patients. Exome capture from paired tumor reference DNA was performed by using a SureSelect XT target enrichment system and SureSelect Human All Exon v3 or v5 bait (Agilent, Santa Clara, CA), according to the manufacturer’s instructions. The prepared libraries were run on a HiSeq 2000/2500 sequencing system (Illumina, San Diego, CA). Candidate somatic mutations were detected using our in-house pipeline for WES (Genomon: http://genomon.hgc.jp/exome/) as previously described.¹¹  All candidates were validated by using polymerase chain reaction (PCR)–based deep resequencing (primer sequences are available upon request).

PCR-amplicon–based target deep sequencing covering 8 genes (PTPN11, NRAS, KRAS, CBL, NF1, SETBP1, JAK3, and SH3BP1) had been performed for 92 of the patients in a previous study.²  In addition, we performed capture-based target deep sequencing in 14 of the 58 newly included patients covering 184 genes as previously described.¹² 

RNA-seq was performed for 129 patients. The quality of the extracted RNA was assessed by using an RNA ScreenTape and a 2200 TapeStation system (Agilent). Sequencing libraries were prepared using an NEBNext Ultra RNA Prep Kit for Illumina (New England Biolabs, Ipswich, MA) with the NEBNext Poly(A) mRNA Magnetic Isolation Module (New England Biolabs) (n = 114) or an NEBNext rRNA Depletion Kit (New England Biolabs) (n = 15), according to the manufacturer’s instructions. Prepared libraries were run on a HiSeq 2500 high-throughput sequencing system. Obtained reads were analyzed using TopHat-Fusion (for gene fusions),¹³  Cufflinks (for expression analysis),¹⁴  HTSeq (for expression analysis),¹⁵  DESeq2 (for expression analysis and differential expression),¹⁶  and VarScan 2 (for nucleotide variations).¹⁷  Candidate gene fusions were validated by reverse transcription-PCR (RT-PCR) using a ThermoScript reverse transcription system (Life Technologies, Carlsbad, CA) and PrimeSTAR GXL DNA Polymerase (TaKaRa Bio, Ohtsu, Japan). We used Cluster 3.0 to group samples on the basis of hierarchical clustering using the average linkage method.¹⁸  Briefly, raw count data for each gene obtained by RNA-seq were normalized by variance stabilizing transformation using DESeq2, and the values were adjusted to center the genes relative to the medians. Finally, clustering was performed using Pearson’s correlation coefficient analysis. Results were visualized using Java TreeView.¹⁹ 

The presence of LIN28B expression was defined by testing the number of reads on LIN28B using Fisher’s exact test against the total number of properly mapped reads. A P value of < .05 indicates a positive result. Thus, 44 of 105 patients with at least 1 read on LIN28B were considered to express LIN28B (supplemental Figure 1).

We analyzed global DNA methylation profiles in 106 patients using the Infinium HumanMethylation450 BeadChip (Illumina) according to the manufacturer’s instructions. We used Cluster 3.0 to group samples and repository data from normal CD34⁺ (n = 5) and acute myeloid leukemia (AML) CD34⁺CD38⁻ samples (n = 14)²⁰  by hierarchical clustering using the average linkage method. We used promoter-associated probes on autosomal chromosomes. β values were adjusted to center the probes relative to medians. Results were visualized using Java TreeView.¹⁹ 

We assessed normalized expression data obtained from RNA-seq and β values obtained from global DNA methylation using gene set enrichment analysis (GSEA) software and the Molecular Signature Database (http://www.broad.mit.edu/gsea/) as previously described.²¹  We used C2 curated gene sets (4729 gene sets), and a false discovery rate of <0.05 was considered to be statistically significant.

Detection of the breakpoint in the RANBP2-ALK fusion gene from genomic and complementary DNA was performed by PCR using PrimeSTAR GXL DNA Polymerase (TaKaRa Bio) and agarose gel electrophoresis. For minimal residual disease detection, we performed real-time quantitative PCR for RANBP2-ALK and RPP30 (ribonuclease P p30 subunit) using an SYBR Premix Ex Taq II reagent (TaKaRa Bio) and an ABI PRISM 7000 sequence detection system (Life Technologies). Primer sequences are listed in supplemental Table 2.

Complementary DNA synthesized by using a ThermoScript reverse transcription system (Life Technologies) was subjected to ALK and ROS1 gene expression quantification using a TaqMan Gene Expression Assays kit (catalog No. Hs01058318_m1 and Hs01090625_m1, Applied Biosystems, Foster City, CA) and ABI PRISM 7000 sequence detection system (Life Technologies), according to the manufacturer’s instructions.

To monitor colony formation ability, 1 × 10³ CD34⁺ cells obtained from patients were cultured for 2 weeks using MethoCult H4434 classic (with cytokines) or MethoCult H4230 (without cytokines) (STEMCELL Technologies, Vancouver, BC, Canada) in the presence or absence of crizotinib, alectinib, ceritinib, TAE684, ruxolitinib, and azacitidine (all from Selleck Chemicals, Houston, TX) dissolved in dimethyl sulfoxide (Sigma-Aldrich, St. Louis, MO), according to manufacturers’ instructions.

Immunohistochemical staining was performed using an anti-ALK antibody (ab17127; Abcam, Cambridge, United Kingdom), an anti-phospho-ALK antibody (PA5-40168; Invitrogen, Carlsbad, CA), and BOND-MAX Fully Automated immunohistochemical and in situ hybridization assays (Leica, Wetzlar, Germany), according to manufacturers’ instructions.

To test the probability of co-occurrence and mutual exclusivity of 2 parameters, we performed a Monte Carlo simulation using the total number of patients and the number of patients who were positive for each parameter, and we documented the number of co-occurrences. We performed a total of 100 000 cycles of simulation using a Mersenne Twister pseudorandom number generator and calculated the P value. The q values were obtained by adjusting P values using Benjamini-Hochberg correction, and a q value of <0.1 was considered to be significant.

To compare the frequency of mutations or other clinical features between the disease groups, categorical variables were analyzed using the χ² test, and continuous variables were tested using the Mann-Whitney U test. Overall survival (OS) and transplantation-free survival (TFS) were calculated using the Kaplan-Meier method; in the TFS analysis, transplantation and death as a result of any cause were censored as events. Hazard ratios for survivals with 95% confidence intervals (CIs) were estimated according to the Cox proportional hazard model. The parameters with P values of < .1 in univariable analysis for OS were incorporated into the multivariable model. The differences in survival were tested using the log-rank test. All statistical analyses were performed using EZR (Saitama Medical Center, Jichi Medical University), which is a graphical user interface for R (The R Foundation for Statistical Computing, Vienna, Austria).²² 

We performed mutational analysis in 150 patients using WES (n = 69) and targeted sequencing (n = 81). Ninety-two patients (13 with WES and 79 with targeted sequencing) were also included in our previous cohort, and each had a unique patient number (UPN) from 1 to 92.²  WES identified a total of 94 somatic mutations, corresponding to an average of 1.3 mutations per patient (supplemental Figure 2). We detected canonical RAS pathway mutations in 134 (89%) of the 150 patients with JMML (Figure 1; supplemental Tables 1 and 3). Although most of the canonical RAS pathway gene mutations were mutually exclusive, coexisting secondary RAS pathway mutations were found in 6 patients (PTPN11 and NF1 in 4 patients, PTPN11 and CBL in 1 patient, and NRAS and CBL in 1 patient). In addition, we identified 4 mutations in other RAS pathway genes: somatic mutations of RRAS2 (which encodes Ras-related protein R-Ras2) in 2 patients (p.Q37H in UPN133 and p.G24D in UPN142) and somatic mutations of SOS1 (which encodes son of sevenless homolog 1) in 2 patients (p.V171E in UPN7 and p.E191K in UPN145) (supplemental Figure 3). Of the 150 patients, 23 (15%) harbored alterations in other genes, including JAK3 (n = 14), SETBP1 (n = 9), and ASXL1 (n = 1). These mutations were more frequent in patients with PTPN11 and NF1 mutations than in other patients (18 of 68 patients vs 4 of 82 patients; P = .002). In a combined analysis of genes within and outside the RAS pathway, 29 patients (19%) who harbored 2 or more gene mutations in recurrently mutated genes had a poor prognosis (supplemental Figure 4), consistent with previous reports.^3,4 

During RNA-seq analysis, we identified 3 in-frame fusions involving receptor tyrosine kinases (DCTN1-ALK in UPN5, RANBP2-ALK in UPN168, and TBL1XR1-ROS1 in UPN106; Figure 2A-C). None of the patients with tyrosine kinase fusion genes harbored canonical RAS pathway mutations; in addition, all patients had monosomy 7 and were relatively older (56, 59, and 153 months) at the time of diagnosis. High levels of ALK and ROS1 expression in patients carrying corresponding gene fusions were confirmed by real-time quantitative PCR (supplemental Figure 5). In addition, the protein expression of ALK and phosphorylated ALK in patients with ALK fusions were also confirmed by immunohistochemical assays (supplemental Figure 6).

Crizotinib, a potent ALK/ROS1/MET inhibitor, yielded an obvious clinical response in patients with ALK/ROS1 aberrant non–small cell lung cancer (NSCLC),²³  inflammatory myofibroblastic tumor,²⁴  AML,²⁵  or ALK-positive anaplastic large cell lymphoma (ALCL).^26,27  To assess the in vitro effect of crizotinib on inhibition of RANBP2-ALK–positive and DCTN1-ALK–positive JMML cells, we performed colony formation assays in the presence of the drug (Figure 3). As expected, crizotinib significantly suppressed granulocyte-macrophage (GM) colony formation of JMML cells; this effect was dose-dependent within the therapeutic concentration range. We translated these findings into clinical application in a patient with RANBP2-ALK (UPN168).

UPN168 (RANBP2-ALK positive) was refractory to conventional AML-type chemotherapy. On the basis of the results of a phase 1 study indicating that crizotinib at 280 mg/m²/day was well tolerated by pediatric patients with refractory solid tumors or ALCL,²⁷  we administered crizotinib in addition to conventional chemotherapy. The patient achieved complete molecular remission indicated by negative ALK break-apart fluorescence in situ hybridization, and negative RANBP2-ALK–specific RT-PCR results. She was successfully bridged to allogeneic bone marrow transplantation from an HLA-matched sibling donor, and she survived without disease recurrence for 15 months after transplantation (Figure 4; supplemental Table 4). Two other patients harboring the ALK or ROS1 fusion did not receive crizotinib and died of tumor progression.

Moreover, given the similar clinical manifestation of RAS pathway mutation–positive JMML and ALK/ROS1-positive JMML, we assessed the effect of crizotinib on JMML cells from 4 patients with canonical RAS pathway gene mutations (CBL, n = 2; PTPN11, n = 1; KRAS, n = 1). Notably, crizotinib inhibited GM colony formation in a dose-dependent manner in all 4 patients (supplemental Figure 7). We also assessed the effect of other ALK inhibitors (alectinib, ceritinib, TAE684) and JAK2 inhibitor (ruxolitinib) on JMML cells from a patient with PTPN11 mutation (c.182A>T, p.D61V). As a result, all tested drugs inhibited GM colony formation with efficacy similar to that of crizotinib (supplemental Figure 8).

We performed a genome-wide methylation analysis of 106 patients and combined these data with repository data from normal CD34⁺ samples (n = 5) and AML CD34⁺CD38⁻ samples (n = 14) for an unsupervised consensus clustering of CpG methylation profiles, which yielded 2 distinct subgroups: a hypermethylation profile and a hypomethylation profile (Figure 5A; supplemental Figure 9). All AML samples fell within the hypermethylation profile cluster, whereas all normal CD34⁺ and NS/MPD samples had the hypomethylation profile. A GSEA revealed marked enrichment of AML-associated genes in the hypermethylation cluster (Figure 5B-C; supplemental Table 5),²⁸  as well as association of this cluster with PTPN11, NF1, 2 or more mutations, older age, higher fetal hemoglobin (HbF) levels, and lower platelet count at diagnosis (supplemental Table 6). The 5-year OS rate was 46.0% (95% CI, 31.0%-59.8%) among patients with the hypermethylation profile and 73.4% (95% CI, 57.6%-84.1%) among patients with the hypomethylation profile (P = .0028). The hypermethylation group also had a significantly poorer 5-year TFS rate than the hypomethylation group (2.2% [95% CI, 0.2%-10.1%] vs 41.2% [95% CI, 27.1%-54.8%]; P = 3.1 × 10⁻¹⁰) (Figure 5D-E).

We assessed the effect of a hypomethylating agent (azacitidine) on JMML cells from patients with hyper- or hypomethylation profiles. As a result, azacitidine was found to inhibit GM colony formation in both group of patients (supplemental Figure 10).

RNA-seq expression profiles of 105 patients were subjected to all-gene-based hierarchical clustering. This identified 4 clusters (supplemental Figure 11), with cluster A1 correlating significantly with PTPN11 mutations and poor TFS. Cluster A4 was associated with NRAS mutations and better survival compared with other clusters (5-year OS, 76.9% [95% CI, 44.2%-91.9%] vs 52.7% [95% CI, 40.3%-63.6%]; P = .12; 5-year TFS, 53.8% [95% CI, 24.8%-76.0%] vs 12.1% [95% CI, 5.2%-22.1%]; P = .003). Hierarchical clustering was also performed using 435 genes associated with a diagnostic classification model originally designed for the classification of myelodysplastic syndrome and leukemia²⁹ ; this identified 3 distinguishable clusters (clusters 1-3) (supplemental Figure 12A). Cluster 1 corresponded to an AML-like expression profile observed in the preceding study, as indicated by correlations of differentially expressed genes⁶  (supplemental Figure 12B), and it correlated strongly with the hypermethylation profile, PTPN11 mutations, NF1 mutations, the existence of 2 or more mutations, cluster A1 identified in all-gene based hierarchical clustering, and poor prognosis (5-year OS, 47.2% [95% CI, 28.2%-64.1%] vs 61.4% [95% CI, 47.2%-72.7%]; P = .20; 5-year TFS, 0% vs 28.5% [95% CI, 17.3%-40.7%]; P = 2.6 × 10⁻⁶] (supplemental Figure 12C-D; supplemental Table 7).

We analyzed differential expression of genes between the hypermethylation and the hypomethylation profiles and identified LIN28B as exhibiting the most substantially elevated expression level in the former (supplemental Figure 13A). LIN28B expression was frequently observed in patients with the hypermethylation profile (P = 2.1 × 10⁻¹³; supplemental Figure 13B) and AML-like expression profile (P = .015). Hypomethylation of a recently identified promoter region in medulloblastoma³⁰  was observed in patients with LIN28B expression (supplemental Figure 14).

We integrated these molecular analyses of JMML, including genome-wide methylation profiling, RNA-seq, and resequencing data. We assessed correlations among methylation and gene expression profiles, genetic alterations, and clinical and other laboratory findings. We found that the hypermethylation profile was highly correlated with most established risk factors for JMML, including an older age, higher HbF, lower platelet count, PTPN11/NF1 mutations, the presence of 2 or more mutations, LIN28B overexpression, and an AML-like expression profile, as well as poor survival (Figure 6; supplemental Tables 8 and 9).

To the best of our knowledge, this is the first integrated molecular analysis of JMML that includes genome-wide methylation profiling, RNA-seq, and resequencing data. Unsupervised consensus clustering of CpG methylation profiles identified patients with the hypermethylation profile, which was significantly associated with an AML-type molecular signature. A recent study showed that AML with wild-type DNMT3A displayed CpG island hypermethylation as a consequence of the progression of AML.³¹  Analogous to AML, JMML with the hypermethylation profile was associated with disease progression and significantly poorer survival.

We analyzed differential gene expression in the hypermethylation and the hypomethylation profiles and detected overexpression of LIN28B and its downstream target HMGA2 with the hypermethylation profile. One previous study on medulloblastoma identified a novel LIN28B promoter whose hypomethylation resulted in high LIN28B expression.³⁰  Hypomethylation of this region was also observed in a JMML patient with LIN28B overexpression. Notably, LIN28B overexpression was recently identified as a poor prognostic marker in JMML.⁷ 

We detected recurrent fusion genes associated with receptor tyrosine kinases (RANBP2-ALK, DCTN1-ALK, and TBL1XR1-ROS1) in 3 (18%) of 16 patients without canonical RAS pathway gene mutations. ALK and ROS1 fusion genes have been reported in other nonhematologic and hematologic malignancies, including NSCLC,^32,33  and RANBP2-ALK and DCTN1-ALK were previously identified as gain-of-function fusions in other malignancies.^34-37  Although the TBL1XR1-ROS1 fusion gene is novel, ROS1 rearrangement with various fusion partners has been reported in several malignancies, including NSCLC³³  and gastric adenocarcinoma.³⁸  Moreover, a TBL1XR1-RET (receptor tyrosine kinase) fusion was reported in thyroid cancer.³⁹  The predicted TBL1XR1-ROS1 fusion protein retained the ROS1 kinase domain, and the ability of TBL1XR1 to form protein dimers strongly suggests that the fusion yielded a gain of function. ALK and ROS1 are RAS/MAPK pathway-activating receptor tyrosine kinases,⁴⁰  and our data suggest that the receptor tyrosine kinase aberration represented an alternative mechanism of molecular RAS/MAPK activation and JMML progression.

To the best of our knowledge, this is the first report of a JMML patient with RANBP2-ALK who achieved a molecular complete response after receiving crizotinib with chemotherapy and was successfully bridged to hematopoietic stem cell transplantation. In previous studies, crizotinib, a potent ALK/ROS1/MET inhibitor, yielded obvious clinical responses in patients with ALK/ROS1-aberrated NSCLC,²³  ALK-positive ALCL,^26,27  and other malignancies. As expected, we found that crizotinib inhibited progression of RANBP2-ALK–positive JMML cells in vitro and in a patient. Furthermore, crizotinib significantly suppressed colony formation by JMML cells that harbored canonical RAS pathway gene mutations. However, we could not determine how ALK inhibitors inhibited colony formation of these cells, and further investigation is necessary.

We detected 4 mutations in other RAS pathway genes including RRAS2 and SOS1. All 4 patients harbored concomitant PTPN11 mutations, suggesting that these newly identified mutations were secondary events. A previous study searched for SOS1 mutations in JMML but found none.⁴¹  SOS1 regulates RAS proteins by facilitating the exchange of guanosine triphosphate for guanosine diphosphate, and germ line SOS1 mutations cause hereditary gingival fibromatosis⁴²  and NS type 4.⁴³  An identical somatic SOS1 missense mutation in UPN145 (c.571G>A, p.E191K) was previously reported in a patient with gastric adenocarcinoma (http://cancer.sanger.ac.uk/cosmic/mutation/overview?id=4094177). However, the functional evaluation of RRAS2 and SOS1 mutations and solitary mutations of JAK3 (p.I688F, p.V765D, and p.L857P, found in a single case each) is desirable to confirm their pathogenicity.

The co-occurrence and mutual exclusivity test identified that various clinical biomarkers may reflect different aspects of the same poor prognostic subgroup of patients with JMML. The hypermethylation profile was also significantly associated with PTPN11/NF1 mutations, clinical features (age, HbF level, and platelet count), and LIN28B expression. These could potentially be used as surrogate biomarkers for this distinct subgroup, which shares several molecular signatures with AML.

In summary, our molecular analyses integrated previously identified and novel biomarkers and defined 2 distinct JMML subgroups with different biological consequences. Our findings should contribute to diagnostics and to the development of therapeutic approaches by providing a model of precise risk stratification.

The authors acknowledge all clinicians, patients, and their families, and thank Yoshie Miura, Yuko Imanishi, and Hiroe Namizaki for their valuable assistance. They also acknowledge the Division for Medical Research Engineering, Nagoya University Graduate School of Medicine, for providing technical support, and the Human Genome Center, Institute of Medical Science, the University of Tokyo (http://sc.hgc.jp/shirokane.html), for providing supercomputing resources.

This work was supported by Japan Society for the Promotion of Science KAKENHI grant No. 22790975, and the Project for Development of Innovative Research on Cancer Therapeutics from the Ministry of Education, Culture, Sports, Science, and Technology of Japan, and from Japan Agency for Medical Research and Development.

Contribution: N.M. performed research, analyzed data, and wrote the paper; Y.O. and H. Muramatsu designed and performed the research, led the project, and wrote the paper; K.Y., Y.S., G.N., K.C., H.T., M.S., and T.U. performed bioinformatic analyses of the resequencing and methylome data; K.S., A.N., H.S., N.K., A.H., M.H., and A.W. collected specimens and performed research; X.W. and Y.X. performed Sanger sequencing; M.I. performed immunohistochemical assays; and S.K., H.A., H. Mano, S.M., S.O., and Y.T. designed the research and analyzed data.

Conflict-of-interest disclosure: N.M., Y.O., and H. Muramatsu have a patent application related to this study. H. Mano served as a scientific adviser for Pfizer Inc. The remaining authors declare no competing financial interests.

Correspondence: Hideki Muramatsu, Department of Pediatrics, Nagoya University Graduate School of Medicine, 65 Tsurumai-cho, Showa-ku, Nagoya, Aichi 466-8560, Japan; e-mail: hideki-muramatsu@med.nagoya-u.ac.jp; and Yoshiyuki Takahashi, Department of Pediatrics, Nagoya University Graduate School of Medicine, 65 Tsurumai-cho, Showa-ku, Nagoya, Aichi 466-8560, Japan; e-mail: ytakaha@med.nagoya-u.ac.jp.