Diverse and Targetable Kinase Alterations Drive Histiocytic Neoplasms

Whole Exome and Transcriptome Sequencing across Histiocytic Neoplasms

To address the above, we performed whole exome sequencing (WES) of frozen tumor biopsies from 24 patients with LCH (n=10) or ECD (n=14) paired with peripheral blood mononuclear cells. Thirteen of 24 patients also underwent RNA sequencing (RNA-seq) (Supplementary Table S1). Both adult (n=15; n=2 with LCH)) and pediatric cases (n=9; n=8 with LCH)) were included. All mutations in activating kinases identified by WES or RNA-seq were validated by droplet-digital PCR (ddPCR), while all others were validated by targeted-capture next-generation sequencing (Supplementary Tables S2–S3 and Supplementary Methods).

Overall, a mean of 7 non-synonymous mutations per adult patient was identified (range 1–22) compared with 5 mutations per pediatric patient (range 4–9; p=ns; Supplementary Fig. S1). The median variant allele frequency (VAF) for the 23 individual known activating kinase mutations identified by WES was 11% while the median VAF for all other somatic mutations was 10% (Supplementary Table S3). Combined WES and RNA-seq revealed mutations or fusions activating MAP kinase signaling in 100% of patients (Fig. 1A; Supplementary Fig. S2A–B). In addition, mutations affecting diverse biological processes co-existed with activating kinase mutations, including recurrent mutations affecting the p38/MAPK and epigenetic regulatory pathways (Supplementary Table S3); however, inspection of the co-occurring genes demonstrated no clear differences in the frequencies of mutated specific genes or pathways of genes affected by mutations between the 24 LCH and non-LCH cases. An analysis of WES data from peripheral blood DNA used as the germline reference for pathologic mutations in a total of 565 cancer-associated genes (including 60 autosomal dominant cancer predisposition genes, MAP kinase pathway members, as well as all genes reported as mutated in the recent studies of clonal hematopoiesis (13–15)) failed to reveal mutations in any of these genes in blood DNA (Supplementary Methods).

Recurrent MAP2K1 and ARAF Mutations in Non-LCH Neoplasms

WES and RNA-seq revealed MAP2K1 mutations in both LCH and non-LCH patients and ARAF mutations in non-LCH patients (Fig. 1; Supplementary Fig. S3A–B). Since neither MAP2K1 nor ARAF mutations have previously been described in ECD, we next interrogated a validation cohort of 37 BRAFV600E-wildtype, non-LCH, FFPE tissue cases (Fig. 1B). Exons 2–3 of MAP2K1 and all coding regions of ARAF were sequenced, along with regions of recurrent mutations in NRAS, KRAS, and PIK3CA, since rare mutations in these genes have previously been noted in BRAFV600E-wildtype, non-LCH patients (16). This revealed recurrent activating mutations in MAP2K1 (32%; n=12), NRAS (16%; n=6), KRAS (11%; n=4), PIK3CA (8%; n=3), and ARAF (3%; n=1) (Fig. 1B–D). Expression of 10 of these MAP2K1 mutants identified activation of MAP kinase signaling over wildtype MAP2K1 (Fig. 1E). Two non-LCH cases demonstrated concurrent activating NRAS and ARAF mutations. Eleven cases (30%; n=11) did not demonstrate known activating mutations in these kinase genes based on targeted gDNA sequencing (Fig. 1A–D; Supplementary Fig. S2 and Supplementary Fig. S3C–D).

Recurrent Kinase Fusions in Non-Langerhans Histiocytic Neoplasms

RNA-seq identified kinase fusions involving BRAF, as well as ALK, exclusively in BRAFV600E-wildtype non-LCH (Fig. 1A; Supplementary Fig. S2). An RNF11-BRAF fusion was detected in an infiltrative, non-LCH brain lesion of a 14-year-old child (Fig. 2A; Supplementary Fig. S4A). Since kinase fusions have not been described in histiocytic neoplasms, we first confirmed expression by reverse transcription polymerase chain reaction (RT-PCR) with breakpoint-flanking primers, as well as with interphase FISH (Fig. 2A–B). Sequencing of the RNF11-BRAF fusion revealed that this transcript developed from an in-frame fusion of exon 1 of RNF11 to exons 11–18 of BRAF. This results in loss of the N-terminal regulatory, RAS-binding domain in BRAF and placement of the BRAF kinase domain under the aberrant regulation of the RNF11 promoter. Although RNF11 is a novel BRAF fusion partner, the RNF11-BRAF fusion has a similar configuration to previously described BRAF fusions (17). Stable expression of RNF11-BRAF in cytokine-dependent, murine pro-B cell Ba/F3 cells, along with an empty vector, wildtype BRAF, or BRAFV600E revealed activation of ERK and MEK phosphorylation by both RNF11-BRAF and BRAFV600E proteins (Fig. 2C). Likewise, expression of RNF11-BRAF or BRAFV600E resulted in cytokine-independent growth of Ba/F3 cells (Fig. 2D), indicating that the RNF11-BRAF fusion is an activating event with pathway activity similar to the BRAFV600E mutation. Moreover, RNF11-BRAF expression sensitized Ba/F3 cells to MEK inhibition but not to RAF inhibition by vemurafenib similar to previous reports analyzing other BRAF fusion genes (18, 19) (Fig. 2E).

Targeted RNA-seq analysis (20) of 9 BRAFV600E-wildtype, non-LCH patients using a panel of 265 genes known to be translocated in cancer detected an additional BRAF fusion. This was a CLIP2-BRAF fusion identified in retroperitoneal lesions from a non-LCH patient, resulting in juxtaposition of the kinase domain of BRAF (exons 11–18) to the N-terminal domain of the protease CLIP2 (Fig. 2F). Expression was confirmed by RT-PCR/Sanger sequencing. As with the RNF11-BRAF fusion, this is a previously unreported fusion of BRAF.

In addition to BRAF fusions, an in-frame fusion of KIF5B (exons 1–24) to the kinase domain of ALK (exons 19–29) was identified in skin lesions from a 25-year-old with BRAFV600E-wildtype ECD (Fig. 2G; Supplementary Fig. S4B). Expression was confirmed by RT-PCR/Sanger sequencing, interphase FISH, and immunohistochemistry (IHC) (Fig. 2G–H; Supplementary Fig. S4C). Although KIF5B-ALK fusions have been described in non-small cell lung cancer (NSCLC) (21), this particular case harbored a unique breakpoint. In both tumor types, the kinase domain of ALK is fused to the N-terminal coiled-coil domain of KIF5B resulting in inappropriate ALK expression, as well as constitutive ALK activation. Introduction of this patient’s ALK fusion into Ba/F3 cells resulted in clear ALK, MAP kinase, STAT3, and PI3K-AKT pathway activation (Fig. 2I), as well as cytokine-independent growth (Fig. 2J). KIF5B-ALK-expressing cells, but not empty vector controls, were exquisitely sensitive to ALK inhibition (Fig. 2K). Interestingly, targeted RNA-seq analysis identified a second KIF5B-ALK fusion in the liver lesions of a 50-year-old ECD patient (Fig. 2L). Finally, an LMNA-NTRK1 fusion was identified in the skin lesions of a 27-year-old ECD patient (Fig. 2M; Supplementary Fig. S4D). LMNA-NTRK1 fusions have been identified in Spitzoid neoplasms (22) and result in aberrant NTRK1 expression with consequent MAP kinase and PI3K-AKT pathway activation.

Gene Expression Analysis of Langerhans and Non-Langerhans Histiocytic Neoplasms

Given the presence of BRAFV600E, MAP2K1, and ARAF mutations in both LCH and non-LCH, we next performed gene expression analysis to understand potential biological differences in these clinically distinct disorders. Unsupervised hierarchical clustering of the top 1% of differentially expressed genes (Supplementary Table S4; Fig. 3A) revealed that samples clustered first by clinical diagnosis (LCH or ECD) then by kinase alteration. Genes differentially expressed between LCH and ECD included all the genes known to encode proteins used to discriminate LCH from ECD histologically (1) (Fig. 3B). Moreover, restriction of gene expression analyses to those LCH and non-LCH samples harboring BRAF mutations revealed the presence of a core set of genes discrepant between LCH and ECD (Fig. 3C–D; Supplementary Fig. S5A–B). Firstly, this included strong enrichment of genes expressed in LCH samples from this study with prior gene expression data from purified CD207+ LCH tumor cells (23). In addition, LCH samples harbored strong enrichment of gene sets upregulated in late-stage myeloid progenitor cells and granulocyte-monocyte progenitors (GMPs), as well as multiple gene sets from classical dendritic cells (cDCs) (Fig. 3C and Supplementary Fig. S5). In contrast, non-LCH demonstrated enrichment of multiple gene sets upregulated in hematopoietic stem cells (HSCs), common myeloid progenitors (CMPs), and monocytes. These differences in gene set enrichments were largely due to upregulation in LCH samples of genes known to be highly expressed in cDCs (24) (including IRF7, RUNX3, GPR82, and CCR7) (Fig. 3D). In contrast, genes known to be upregulated in HSCs such as BAALC (25) and ID1 (26), as well as core macrophage-associated genes (such as CEBPA (27)) were upregulated in non-LCH samples compared with LCH (Fig. 3D). In addition to gene sets related to hematopoietic ontogeny, gene sets involved in cell cycle regulation and IL1 signaling were found to be enriched in LCH but not ECD. Conversely, genes involved in lipid metabolism and adipogenesis were uniquely enriched in ECD but not LCH (Supplementary Fig. S5). Although this result is consistent with the histologic appearance of ECD lesions, these data highlight a potential role of adipogenesis not previously recognized in this disorder.

Therapeutic Targeting of MAP2K1- and ARAF-mutant Refractory Histiocytoses

Three of the investigated non-LCH patients with refractory disease and progressive organ dysfunction were treated with targeted therapies based on the discovery of novel kinase alterations described above. A MAP2K1K57N mutation was noted from WES of peri-renal lesions from a 53-year-old ECD patient with progressive disease following treatment with both interferon-alpha and anakinra, who was symptomatic by way of inflammatory ascites and renal failure. Treatment was initiated with a MEK inhibitor (single-agent trametinib 2mg daily), resulting in abrupt cessation of ascitic accumulation and normalization of creatinine. Additionally, there was metabolic resolution of FDG-avid infiltrates in the retroperitoneum (SUV 8.2) and spermatic cords (SUV 6.2) to background SUV (Fig. 4A–B). The patient has been maintained on trametinib single-agent therapy with a sustained clinical response for >180 days. Similarly, treatment of a second ECD patient with MAP2K1Q56P-mutant disease refractory to 4 lines of prior therapy led to resolution of PET-avid disease in renal (SUV 8.4), aortic (SUV 3.9), and maxillary sinus (SUV 9.5) infiltrations to background SUV within a single month of administration of the MEK inhibitor cobimetinib (Fig. 4C).

Further evidence of effective targeted inhibition was found in an ECD patient, with disease in the bones, orbits, cavernous sinuses, and choroid, whose tumor was found to harbor an ARAFS214A mutation. This patient’s disease had progressed following cladribine, clofarabine, and anakinra. The patient was also symptomatic with visual impairment from macular edema that required ongoing corticosteroids (Fig. 4D; Supplementary Fig. S3A–B). Given a recent report of complete response to sorafenib in a NSCLC patient with an ARAFS214C mutation (28), we initiated sorafenib therapy, titrated up to 600mg daily. Within 12 weeks, there was regression of lesions in the cavernous sinuses (by post-gadolinium MRI) and retina (as visualized directly by fundoscopy), and the patient was able to taper her steroid dose. This coincided with a >50% decrease in mutant ARAF DNA in plasma cell-free DNA (Fig. 4E).

Patients

The study was conducted according to the Declaration of Helsinki, and human tissues were obtained with patient-informed consent under approval by the Institutional Review Boards of Memorial Sloan Kettering Cancer Center, St. Jude’s Children’s Research Hospital, MD Anderson Cancer Center, the National Human Genome Research Institute, and Pitié-Salpêtrière Hospital.

Excised lesions were either flash-frozen for DNA/RNA extraction and/or fixed in 4% neutral-buffered formalin, embedded in paraffin, and processed by routine histological methods. For patients undergoing whole exome sequencing (WES), DNA extracted from peripheral blood mononuclear cells was utilized as a paired germline control. In total, specimens from 64 patients were analyzed, and the clinical, histological, and genetic characteristics are summarized in Supplementary Tables S1 and S3.

Whole Exome Sequencing

Analysis of exome sequencing data, which includes mapping, coverage and quality assessment, SNV/Indel detection, tier annotation for sequence mutations, and prediction of deleterious effects of missense mutations have been described previously (40, 41). Approximately 250ng of DNA from each sample was sheared to an average of 150 bp in a Covaris instrument for 360 seconds (Duty cycle - 10%; intensity - 5; cycles/Burst - 200). Barcoded libraries were prepared using the Kapa Low-Throughput Library Preparation Kit Standard (Kapa Biosystems), amplified using the KAPA HiFi Library Amplification kit (Kapa Biosystems) (8 cycles), and quantified using Qubit Fluorimetric Quantitation (Invitrogen) and Agilent Bioanalyzer. An equimolar pool of the 4 barcoded libraries (300ng each) was used as input to capture the exome using one reaction tube of the Nimblegen SeqCap EZ Human Exome Library v3.0 (Roche, cat # 06465684001), according to the manufacturer’s protocol. The pooled capture library was quantified by Qubit (Invitrogen) and Bioanalyzer (Agilent) and sequenced on an Illumina HiSeq 2500 using a paired end, 100 nucleotide in length run mode, to achieve an average of 100X coverage.

Confirmation of Mutations

The following mutations in genes coding tyrosine kinases were confirmed using ddPCR: ARAF S214A, S186R, A225V, and P539H; BRAF V600E and R603Q; MAP2K1 K57N, F53_Q58del, Q58_E62del, and F68L; NRAS Q61R and PIK3CA E542K, using BioRad probes with FAM for the mutants and HEX for the wildtype sequences, on a BioRad QX200 ddPCR system, following the manufacturer’s instructions, starting from 30ng of template DNA. For the rest of the mutations, we used a custom-designed, TruSeq Custom Amplicon probe to confirm the mutations detected by exome sequencing. Design Studio (Illumina, San Diego, CA) was used to design amplicons covering the regions of interest. The regions were amplified using 250ng of template genomic DNA, using the manufacturer’s instructions, with 25 cycles of amplification, and were run on an Illumina MiSeq 2 × 250 cartridge.

RNA Sequencing

RNA library construction for transcriptome sequencing was done as per manufacturer’s instructions using the Illumina TruSeq RNA sample preparation V2. Sequencing was completed on the Illumina HiSeq 2000 as per manufacturer’s instructions. Analysis of transcriptome sequencing data, which includes mapping, coverage, and quality assessment, SNV/Indel detection, tier annotation for sequence mutations, and prediction of deleterious effects of missense mutations have been described previously (40, 41). For gene expression analyses, transcript expression levels were estimated as Fragments Per Kilobase of transcript per Million mapped reads (FPKM); gene FPKMs were computed by summing the transcript FPKMs for each gene using Cuffdiff2 (42, 43). A gene was considered “expressed” if the FPKM value was ≥ 0.5 based on the distribution of FPKM gene expression levels. Genes that were not expressed in any sample group were excluded from the final data matrix for downstream analysis. Gene set enrichment analysis (GSEA) was performed as described previously (44), and the methodology and approach are described in the Supplementary Methods.

Targeted DNA Sequencing

We sequenced the regions of known mutations in MAP2K1, NRAS, KRAS, and PIK3CA, as well as all coding exons of ARAF using Sanger sequencing (primer sequences listed in Supplementary Table S5) and/or hybrid-capture, next generation sequencing using the MSKCC IMPACT assay as previously described (45) or the Foundation One Assay (Foundation Medicine, Inc Cambridge, MA) as previously described (46). Prior to DNA extraction, FFPE samples from all cases were reviewed to confirm that the tissue was of sufficient size to generate a minimum of 50ng of 20% histiocyte nucleic acid. DNA was isolated from 40-μm-thick sections of FFPE tissue.

Targeted RNA Sequencing

Total RNA extracted from 40-μm-thick sections of FFPE tumor was reverse transcribed with random hexamer primers using the Super-Script III First-Strand Synthesis System (Invitrogen). Double-stranded cDNA was synthesized with the NEBNext mRNA Second Strand Synthesis Module (New England Biolabs). Hybrid selection of indexed, adaptor-ligated libraries was performed using the cDNA Kinome hybridization kit with 612 transcripts of kinases and kinase-related genes (Agilent SureSelect Human Kinome Kit). Selected libraries were sequenced on the HiSeq-2000 instrument (Illumina) with 49 × 49 paired reads. For RNA sequencing, we used a sequencing approach targeting 612 transcripts of kinases and kinase-related genes. We aimed for a high number of unique read pairs (~50,000,000) per sample (Supplementary Table S2).

Gene Fusion Confirmation

All gene fusions were validated with RT-PCR followed by direct sequencing, interphase fluorescence in situ hybridization (FISH), and/or immunohistochemical analysis. RT-PCR was performed from cDNA (primer sequences available upon request), followed by analysis on a Bioanalyzer (Agilent). Specific PCR amplicons were only detected with the appropriate combination of primer and template and not with negative controls. The nucleotide sequence at the fusion site was confirmed with Sanger sequencing. Details of interphase fluorescence in situ hybridization (FISH) are included in the Supplementary Methods.

Plasma Cell-Free DNA (cfDNA) Collection and Analysis by Droplet Digital PCR (ddPCR)

10mL of blood was collected into Streck tubes, and plasma was then separated from blood using standard techniques. Plasma cfDNA was isolated using the QIAamp Circulating Nucleic Acid Kit (QIAGEN; Germantown, MD) according to the manufacturer’s instructions, and the concentration was assessed using a BioAnalyzer.

The S214A (c.640T>G) mutation for ARAF was evaluated by ddPCR (BioRad QX200 Hercules, CA) in a custom-designed, allele-specific assay (primer sequences available upon request). The experiments were performed using the following protocol: 1 cycle at 95°C for 10 minutes, 40 cycles at 94°C for 30 seconds and 55°C for 1 minute, 1 cycle at 98°C for 10 minutes, then 1 cycle at 4°C infinite, all at a ramp rate of 2°C/second. Bio-Rad’s T100 thermal cycler was used for the PCR step. When available, 5ng of DNA was assessed in a 20μl PCR reaction, partitioned into approximately 20,000 droplets. A total of two replicates were used per sample. Droplets were quantified using the BioRad Quantasoft Software (version 7.0) to identify the concentration of ARAFS214A mutant copies/mL DNA and of ARAF wild type copies/mL of DNA. The ratio of the concentration of mutant to wild type ARAF ([Mu]/[WT]) was then calculated for each patient pretreatment and following 1 and 3 months of sorafenib therapy.

Germline DNA Mutational Analysis

An automated analysis of mutations from the WES of germline DNA was performed to search for pathologic mutations in 565 cancer associated genes, including 60 autosomal dominant cancer predisposition genes (ALK, APC, BAP1, BMPR1A, BRAF, BRCA1, BRCA2, CBL, CDC73, CDH1, CDK4, CDKN1C, CDKN2A, CEBPA, DICER1, EPCAM, FH, GATA2, HRAS, KRAS, MAP2K1, MAP2K2, MAX, MEN1, MLH, MSH2, MSH6, NF1, NF2, NRAS, PALB2, PAX5, PHOX2B, PMS2, PRKAR1A, PTCH1, PTEN, PTPN11, RAF1, RB1, RET, RUNX1, SDHA, SDHAF2, SDHB, SDHC, SDHD, SHOC2, SMAD4, SMARCA4, SMARCB1, SOS1, STK11, SUFU, TMEM127, TP53, TSC1, TSC2, VHL, and WT1), as well as all genes previously reported to be mutated in individuals with clonal hematopoiesis associated with a somatic mutation (ASXL1, DNMT3A, FLT3, GNAS, GNB1, IDH1, IDH2, JAK2, KIT, NPM1, SF3B1, SRSF2, TET2, U2AF1, and others (13–15, 47)). With the exception of variants that matched those contained within the IARC TP53 database, only novel, non-silent single nucleotide variants (SNVs) or SNVs with <0.1% population frequency in the NHLBI ESP database received an automatic classification. In-frame indels that matched dbSNP did not receive a classification. Indels present in 1000 Genomes or with multiple submissions in dbSNP were also excluded.

To identify pathologic lesions, variants were classified into three tiers, named gold, silver, and bronze. Gold is reserved for truncation mutations in tumor suppressor genes, matches to truncation mutations and hotspot mutations in somatic mutation database, and perfect matches to highly-curated locus-specific databases, including the IARC TP53 database, NHGRI BRCA1 and BRCA2 database for records marked clinically important, ARUP MEN2 database for mutations in RET, ASU database for TERT mutations, LOVD for mutations in APC and MSH2, and the RB1 mutation database. Silver is used for variants of less certain significance, and bronze for variants more likely to be benign. Additional databases that were utilized for classification included HGMD, ClinVar, and UMD, matches to which generally received a silver classification. Truncation mutations located close to the C-terminus of the protein that do not map to a functional domain were examined manually and consulted with external gene experts for their pathogenicity.

Drug Studies

Vemurafenib, sorafenib, GDC-0973 (cobimetinib), crizotinib, and alectinib were purchased from Selleckchem (Houston, TX, USA). Drug studies were conducted in vitro using FACS-sorted, DAPI⁻ eGFP⁺ Ba/F3 cells that stably expressed the MIGII-EV, MIGII-BRAF V600E, MIGII-RNF11-BRAF, and MIGII-KIF5B-ALK constructs using the CellTiter-Glo® Luminescent Cell Viability Assay from Promega Corporation (Madison, WI, USA), according to the manufacturer’s instructions. The MIGII-BRAF V600E, MIGII-RNF11-BRAF, and MIGII-KIF5B-ALK FACS-sorted Ba/F3 cells were maintained in RPMI + 10% Fetal Bovine Serum (FBS) + penicillin and streptomycin media without murine IL-3 while MIGII-EV was maintained in RPMI + 10% FBS + penicillin and streptomycin with recombinant murine IL-3 (1 ng/mL).

Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) was performed as described previously³. The dataset was converted to gene symbols, and the gene expression signatures were analyzed using the java GSEA package. The most differentially expressed genes for each comparison were used to generate a signature for GSEA analysis. The input motif gene sets were extracted from the Molecular Signature Database, version 4 (MSigDBv4). Gene sets with a p<0.05 and a false discovery rate (FDR) q-value < 25% were considered to be significantly enriched in LCH and non-LCH cases evaluated by RNA-seq.

Accession Codes

Whole exome (n=24 patients) and transcriptome (n=13) sequencing data for the patients from the discovery cohort have been deposited in the Sequence Read Archive (SRA) and Gene Expression Omnibus (GEO), respectively. The SRA accession number is SRP065600, and the GEO accession number is GSE74442.