Somatic activating mutations in MAP2K1 cause melorheostosis

Melorheostosis is a sporadic disease of uncertain etiology characterized by asymmetric bone overgrowth and functional impairment. Using whole exome sequencing, we identify somatic mosaic MAP2K1 mutations in affected, but not unaffected, bone of eight unrelated patients with melorheostosis. The activating mutations (Q56P, K57E and K57N) cluster tightly in the MEK1 negative regulatory domain. Affected bone displays a mosaic pattern of increased p-ERK1/2 in osteoblast immunohistochemistry. Osteoblasts cultured from affected bone comprise two populations with distinct p-ERK1/2 levels by flow cytometry, enhanced ERK1/2 activation, and increased cell proliferation. However, these MAP2K1 mutations inhibit BMP2-mediated osteoblast mineralization and differentiation in vitro, underlying the markedly increased osteoid detected in affected bone histology. Mosaicism is also detected in the skin overlying bone lesions in four of five patients tested. Our data show that the MAP2K1 oncogene is important in human bone formation and implicate MEK1 inhibition as a potential treatment avenue for melorheostosis.

M elorheostosis (OMIM%155950) is a rare dysostosis characterized by excess bone formation in the classic "dripping candle wax" pattern on the surface of bone, revealed radiographically (Fig. 1a, b and Supplementary Fig. 1) 1 . The disease occurs in males and females equally, without familial clustering, and often affects multiple contiguous bones unilaterally in an asymmetric distribution 2 . Symptoms of melorheostosis may begin in childhood; the diagnosis is usually apparent by the age of 20 years 3 . The bone overgrowth lesions are associated with pain, functional impairment, joint contracture, and deformity. There is a wide spectrum of clinical appearance 2 with skin lesions overlying the affected bone noted in up to 30% of patients. Bone lesions do not metastasize, but conversion to osteosarcoma has rarely been reported 4 . Increased uptake on bone scintigraphy suggests that melorheostotic bone is metabolically active (Fig. 1c-e). Patients are diagnosed by a combination of clinical findings, radiographs and bone scan. There is no definitive diagnostic test or specific treatment for melorheostosis 2 .
At present, two patterns of occurrence of melorheostosis have been differentiated by genetics. Patients with the autosomal dominant "spotted bone disease" osteopoikilosis and the related condition Buschke-Ollendorff syndrome (OMIM #166700) may display features of melorheostosis, and have been shown to have germline loss-of-function mutations in LEMD3, encoding the inner nuclear membrane protein MAN1 5 . LEMD3 mutations have not been found in patients with the more common sporadic melorheostosis without osteopoikilosis 6 ; the genetic cause of their condition is undefined 7 . Recently, a patient with melorheostosis and familial osteopoikilosis was found to have a germline LEMD3 mutation and a somatic KRAS mutation in an overlying area of scleroderma-like skin 8 . Whether the KRAS mutation was causative for the bone findings in this one patient remains undetermined. The lack of vertical genetic transmission and asymmetric involvement suggests that sporadic melorheostosis may be caused by somatic mutations in bone-forming cells 9 , but bone tissue has not been previously investigated for somatic mutations in bone forming cells.
We biopsied affected and unaffected bone of 15 patients with melorheostosis and compared the tissue whole exome sequences for each patient. We identified mutations in the MAP2K1 gene in affected, but not unaffected, bone of eight patients. MAP2K1 encodes the protein MEK1, a kinase whose activity is modulated by its negative regulatory domain. The MAP2K1 mutations identified in melorheostosis patients cause substitutions in two residues of the MEK1 negative regulatory domain. We present genetic, functional and histological data supporting the enhancement of MEK1 activity predicted by the location of the mutations and the causative role of the MAP2K1 mutations in melorheostosis. This study provides evidence of a bone disease caused by mutations in MAP2K1.

Results
Identification of mosaicism for MAP2K1. Fifteen patients with melorheostosis underwent biopsies of affected and unaffected bone (Supplementary Table 1 and 2). A consistent intraoperative finding was extremely dense, rigid affected bone that often dulled the osteotomes and drill bits. Genomic DNA was extracted from both tissues and subjected to high-depth whole exome sequencing. We identified 8365 sequence variants that were present in a subject's affected bone but not unaffected bone. After restricting these variants to those that were found in more than 1% of sequences (but less than germline frequency), rare in the general population, and predicted to alter protein coding or mRNA splicing, we were left with 284 variants (Fig. 2a and Table 1). Only five genes contained putative somatic mutations in more than one affected bone sample (MAP2K1, USH2A, CCDC13, SIRT5 and RAB44). Among these five genes, it was striking that sequences from five affected bone samples contained one of three tightly clustered missense mutations in the coding region of MAP2K1, c.167 A > C (p.Q56P), c.169 A > G (p.K57E), and c.171 G > T (p.K57N), encoding substitutions in two adjacent amino acids of the protein MEK1 ( Fig. 2b and Supplementary Data 1). These coding variants are not present in the exomes of 123,136 healthy individuals searchable with the ExAC 10 . Both substituted residues are located in a MEK1 α-helix that negatively regulates kinase function ( Fig. 2b and Supplementary Fig. 2) 11 . All three of these mutations have been previously shown to lead to gain-of-function through loss of the inhibitory role of the negative regulatory domain and have been identified in malignancies, including lung 12 , melanoma 13 and hairy cell leukemia 14 (Supplementary  Table 3). All other 281 variants were unique to the affected bone of 1 of the 15 subjects (Supplementary Data 1). No germline or somatic LEMD3, MAP2K2 or MAP2K3 variants were identified in automated or manual re-examination of exome sequencing data Fig. 1 Clinical findings in melorheostosis. a Clinical appearance of Melo-10 with melorheostosis of right lower extremity. Note irregular thickening of the right leg. The affected bone was found to harbor a MAP2K1 mutation. b Radiograph of the right tibia/fibula of Melo-10 with classic candle-wax appearance. c Maximal intensity projection (MIP) 18 F-NaF PET image of the lower extremities showing intensely increased 18 F-NaF activity in the bones of the right leg, primarily along the medial and distal femur, and along the tibia extending into the foot. d, e Axial CT and fused 18 F-NaF PET/CT images at the level of mid-calf demonstrating dense cortical bone formation in the anterior and posterior aspect of the right tibia (red arrows in d with corresponding abnormally increased 18 F-NaF uptake (white arrows, in e)  from any patient. A KRAS p.Q61R mutation was found in both the affected and unaffected bone of one MAP2K1 mutationnegative patient with skin lesions consistent with RASopathy, who likely has melorheostosis as part of a more complex early embryonic mosaic RASopathy, which will be reported separately. Thus 5 of 15 patients were identified to have MAP2K1 mutations by whole exome sequencing. The allele frequency of these somatic heterozygous MAP2K1 mutations identified by automated somatic mutation calling was relatively high (10-25%). We performed ampliconbased targeted sequencing to validate the exome sequencing findings and screen the remaining patients for MAP2K1 mutations at higher sensitivity. This approach revealed three additional patients with MAP2K1 p.K57N mutations at lower mutant allele frequencies (Table 1). Amplicon-seq polymerase errors were less than 1% at the base positions in question (Supplementary Tables 4, 5 and 6). Their mutations were missed by WES due to low allele frequencies in patients Melo-16 and Melo-6 (2.7% and 4.8%, respectively) or low coverage in patient Melo-2 (58×). All eight MAP2K1 mutations were then confirmed by droplet digital PCR (ddPCR) of DNA from affected bone (Fig. 2c). The MAP2K1 mutant allele frequency in affected bone by ddPCR ranged from 9 to 28% for p.Q56P, 3-34% for p.K57N, and 18% for p.K57E and correlates well with amplicon and WES data.
The ddPCR mutant allele frequency in osteoblasts cultured from melorheostotic bone of seven subjects, varied from 0.05 to 47% and was similar in different samples from the same cell culture, but often did not match the frequency found in bone tissue, consistent with tissue mosaicism. For example, there was a markedly higher frequency of the p.K57N mutation in cultured osteoblasts than bone tissue from patient Melo-2 (40-47% in cell passages 1-5 vs. 6.9% in bone) and a lower frequency of p.Q56P in Melo-4 (0.05-0.2% in cell passages 1-5 vs. 10% in bone tissue).
Mosaicism for their respective MAP2K1 mutations was also detected in skin tissue overlying affected bone of patients Melo-2, Melo-4, Melo-6 and Melo-18, but not detected in skin from patient Melo-16, who has a lower disease burden. To look for germline mutations, we tested blood samples available from seven patients. MAP2K1 mutations were not identified in blood samples including those with mutation-positive skin tissue, further confirming mosaicism (Table 1).
Melorheostosis histology reveals intense remodeling. To gain insights into the effects of MAP2K1 mutations on bone structure, remodeling and mineralization, we compared histomorphometry of tissue sections from affected and unaffected bone in six patients. Melorheostotic bone was characterized in the outer regions by distinctive parallel layers of primary lamellar bone, an organization that underlies the surgical hardness of the bone (Fig. 3a). This newly formed compact tissue was intensely remodeled into a highly porous osteonal-like bone at greater depth from the surface (Fig. 3b). Compared to unaffected bone, affected bone showed increased active bone-resorbing osteoclasts, elevated eroded surfaces and a higher number of bone-forming osteoblasts (Fig. 3c). There was an approximately six-fold increase in the thickness of osteoid (unmineralized bone matrix) and a greater than 50-fold increase in osteoid surface/bone surface in the affected bone samples compared to their respective unaffected counterparts (Fig. 3c, d).

MEK1 signaling increased in melorheostosis.
To determine how the MAP2K1 mutations identified in our cohort affected the MEK1 target kinases ERK1 and ERK2, we examined MEK1 activity by immunohistochemical analysis of phosphorylated ERK1/2 (p-ERK1/2) in sections of affected bone tissue from two patients with MAP2K1 mutations. Osteocytes and osteoblasts in affected bone displayed stronger p-ERK1/2 signals than control bone, but no p-ERK1/2 signal was detected in osteoclasts ( Fig. 4a). Consistent with somatic mosaicism, the p-ERK1/2 positive cells occurred in patches in affected bone tissue. Furthermore, osteoblasts cultured from affected bone and having a high frequency of mutant MAP2K1 cells (approximately 80%) exhibited a bimodal ERK1/2 activation after serum-stimulation when examined by flow cytometry for intracellular p-ERK1/2. Bimodal ERK1/2 activation was not observed in cells from unaffected bone. Both peaks of ERK1/2 activation in osteoblasts from affected bone, likely representing wild-type and mutant cells, were reduced by the MEK-specific inhibitor U0126. Baseline p-ERK1/2 was also higher in a proportion of affected cells than in unaffected cells (Fig. 4b). These findings would be expected from a mixed cell population and confirm MAP2K1 mosaicism at the cellular level. We also investigated the MEK1-ERK1/2 pathway in cultured osteoblasts harboring MAP2K1 mutations by western blot analysis in seven patients (Fig. 4c). Osteoblasts from the same patient and passage as those analyzed by flow cytometry showed increased p-ERK1/2 compared to cells from unaffected bone, indicating MEK1 gain-of-function ( Supplementary Fig. 4). FBS-induced ERK1/2 activation could more readily be seen in unaffected osteoblasts by western blotting, likely due to the increased sensitivity in comparison to flow cytometry. Pre-treatment of cells with U0126 significantly diminished serum-stimulated p-ERK1/2, supporting enhanced MEK1 activity as mediating the increased ERK1/ 2 activation. The level of ERK1/2 activation by MEK1 generally correlated with mutant allele frequency. Activation was clearly seen in Melo-2 and Melo-18, with the highest percentage of mutant cells. In cells cultured from affected bone of Melo-10 with lower mutant allele frequency of 10%, the activation of ERK1/2 by MEK1 was equivocal ( Supplementary Fig. 4).
MEK1 activation increases cell proliferation. MAP2K1 mutations did not affect MAP2K1 transcripts (Fig. 5a) or MEK1 protein levels (Fig. 5b), supporting a gain-of-function mechanism intrinsic to the protein itself. ERK1/2 activation promotes G1 to S-phase cell cycle progression and accelerates cell proliferation. The enhanced MEK1-ERK1/2 activation in affected osteoblasts increased cell proliferation in vitro (Fig. 5c). For example, the doubling time of osteoblasts from affected bone of patient Melo-2 (mutant allele frequency 45%) was 18 h, while osteoblasts from this patient's unaffected bone had a doubling time of 54 h. Osteoblasts from affected bone did not lose contact inhibition. Furthermore, affected osteoblasts showed increased levels of cyclin D3 expression compared to unaffected cells, consistent with accelerated G1-S cell cycle transition (Fig. 5d).

MEK1 activation decreases osteoblast differentiation in vitro.
To further investigate the functional consequences of the MAP2K1 mutations causing melorheostosis, osteoblast differentiation and mineralization were assessed in vitro. The presence of the MAP2K1 mutation strikingly inhibited BMP2-stimulated mineralization in affected compared to unaffected osteoblasts (Fig. 5e). Expression of marker genes for osteoblast differentiation was also reduced in cells with the MAP2K1 mutation during BMP2-stimulated differentiation (Fig. 5f). Reduced expression of RUNX2 (runt-related transcription factor 2), ALPL(alkaline phosphatase) and COL1A1 (the α1 chain of type I collagen) indicates the prevalence of immature bone forming cells. Osteoblasts early in differentiation secrete proteins to stimulate osteoclastogenesis 15 . In fact, we observed a markedly increased ratio of RANKL/OPG transcripts in osteoblasts from affected bone, which would generate a strongly positive context for osteoclastogenic stimulus (Fig. 5f) 16 . Together, these data support key histological findings in affected bone. Inhibition of mineralization by overactive MEK1 accounts for the massive accumulation of unmineralized osteoid in affected bone tissue (Fig. 3d). Furthermore, the elevated RANKL/OPG ratio stimulates the increased osteoclast number and intense remodeling that occurs in melorheostotic bone (Fig. 3c).

Discussion
Genetic osteosclerotic disorders are a rare but interesting subset of bone diseases which can provide important insights into bone biology. Identification of a germline mutation causing sclerosteosis led to the development of a potential treatment for osteoporosis [17][18][19] . Investigations into fibrodysplasia ossificans progressiva elucidated the role of inflammation in heterotopic ossification and exposed targets for treatment 20 . Here, our investigation has revealed a role for MEK1 in a disorder characterized by excessive bone formation.
We identified mutations in MAP2K1 in the affected, but not unaffected, bone of 8 of 15 patients diagnosed clinically and radiologically with melorheostosis. The mutations cluster at two residues in the MEK1 negative regulatory domain, p.Q56 and p. K57, where identical substitutions have been identified in multiple malignancies, but not primary bone tumors. The occurrence of the mutations in a MEK1 "hot spot," and the evidence of enhanced MEK1-ERK1/2 signaling in melorheostotic bone and osteoblasts, strongly support MAP2K1 somatic mutations as causing about half of cases of sporadic melorheostosis. In our cohort, we found no occurrences of somatic or germline mutations in LEMD3, which have previously been reported in patients with familial osteopoikilosis and features of melorheostosis 6,8 .
Mosaicism for cells with increased ERK1/2 activation was confirmed directly in melorheostotic bone by immunohistochemistry. Furthermore, two populations of cells with distinct levels of ERK1/2 activation were detected by flow cytometry. Increased ERK1/2 activation was also detected in cultured osteoblasts from affected, but not unaffected bone. In individual patients, the proportion of cultured mutant osteoblasts often differed from the proportion of mutant allele in genomic DNA extracted from bone chips from the same lesion, consistent with clustering of affected cells in mosaic conditions.
The three MAP2K1 mutations we identified in melorheostosis have been shown to increase MEK1 activity in transfection studies 21 , and likely act by destabilizing the α-helix comprising the MEK1 negative regulatory region, which normally keeps the ATPbinding site in an inactive conformation 11 . MEK1 is a downstream activator in the RAS pathway. Mutations in RAS pathway components are associated with developmental defects falling under the umbrella designation RASopathies 22 and involve variable cardiac, facial and neurodevelopmental defects. Interestingly, germline mutations identified in leukocytes from patients with the RASopathy cardio-facio-cutaneous syndrome 23, 24 include two MAP2K1 mutations (p.G128V/p.Y130C) in the catalytic core, and putatively activating mutations in the negative regulatory region of MAP2K1 (p.F53S) or MAP2K2 (p.F57V, analogous to MAP2K1 p.F53). Bone overgrowth is not reported in these RAS-MAPK syndromes, implying that the somatic mutations we report in osteoblasts produce a different phenotype than the same mutation in the germline. Mosaicism for mutations at the same residues (p.Q56P and p.K57E) we identified in the MEK1 negative regulatory domain in While the MAPK pathway is known to be involved in skeletogenesis and MEK1 is essential in embryogenesis, our study documents the importance of the MEK1-ERK1/2 pathway in human bone cell biology. While somatic mutations in MAP2K1 have also been reported in Langerhans cell histiocytosis and   . d Western blot analysis shows that affected osteoblasts (lanes 2, 4, and 6) expressed higher level of cyclin D3 compared to unaffected (lanes 1, 3, and 5), correlating with the increased p-ERK1/2 level in affected cells shown in Fig. 4c. Note that levels of cyclin D3 and p-ERK1/2 decreased on day 3 (D3) and day 6 (D6) compared to day 1 (D1), because culture media was not refreshed after day 1, similar to conditions used during the live-cell imaging shown in Fig. 5c. e Alizarin Red S staining of mineralization in osteoblast cultures from patient Melo-18 (MAP2K1 p.K57N, VAF 46%). After 7 weeks of BMP2-stimulated mineralization in vitro, significantly inhibited mineralization was observed in affected cells compared to unaffected. f Realtime qPCR analysis of expression of osteogenic marker genes, RUNX2, COL1A1, and ALPL in osteoblasts from melorheostosis patient Melo-2 (MAP2K1 p. K57N, VAF 45%) after two weeks of osteogenic stimulation. Expression level of RUNX2, COL1A1, and ALPL was significantly lower in affected osteoblasts compared to unaffected. g The RANKL/OPG transcript ratio, an index of osteoclastogenic stimulus, was assessed by real-time qPCR with Melo-2 patient osteoblasts as in Fig. 5f. Significantly higher ratio of RANKL/OPG in affected osteoblasts indicates increased osteoclastogenesis compared to unaffected Erdheim-Chester disease (ECD) 26 , the mutated cells in ECD are derived from monocytes rather than osteoblasts. Bone lesions in ECD are typically bilateral and symmetrical and occur in trabecular bone of the epiphysis 27 . In contrast, activation of ERK1/2 MAP kinases in melorheostosis occurs in osteoblast lineage and the lesions are primarily cortical. Transgenic murine models, in which the entire MEK1 negative regulatory domain was deleted, share the elevated ERK1/2 signaling seen in melorheostosis 28 . However, the mice displayed accelerated bone remodeling rather than bone overgrowth. The mouse osteoblasts also showed increased differentiation and mineralization, but unchanged proliferation in vitro.
In melorheostosis, activation of the MEK1-ERK1/2 pathway in mature bone results in an increase in osteoblast surface and increased production of unmineralized bone matrix (osteoid). This is consistent with the faster proliferation of osteoblasts cultured from melorheostotic bone and with the less differentiated osteoblast phenotype revealed in vitro in BMP2stimulated cells, in which transcripts of osteogenic genes RUNX2, COL1A1, and ALPL were decreased. Histologically, we observed a significant increase in osteoid formation and osteoclast numbers in affected bone from the MAP2K1 variant cohort. Interestingly, these features were previously noted in a single case report 29 with unknown genetics. The increased osteoid seen in vivo is consistent with the decreased mineralization deposited by affected osteoblasts in vitro with BMP2 stimulation. Increased osteoclast numbers and bone turnover in affected bone is consistent with the increased RANKL secretion by cultured cells, which would stimulate osteoclast development and intense bone remodeling. Furthermore, it is known that chronic activation of ERK1/2 strongly inhibits BMP2-mediated mineralization, as previously shown in a murine model for neurofibromatosis 30 , and provides an underlying mechanism for the increased unmineralized osteoid found in melorheostotic bone. The phenotypic distinction between murine models and melorheostosis patients may be ascribed to overexpression of the MEK1 transgene in all murine osteoblasts as opposed to the mosaic mutations in melorheostosis, or greater disruption of function by deletion of the full negative regulatory domain 28,31 . Furthermore, studies in fruit flies and zebrafish indicate that activating MAP2K1 mutations affect ERK1/2 activation differently depending on pathways intrinsic to the target tissue and developmental state 22 .
Identification of a somatic mosaic mutation in osteoblasts accounts for the unique clinical features of melorheostosis. The "dripping candle wax" pattern observed on radiographs is due to marked appositional bone growth followed by aggressive remodeling. The radiographic pattern of multiple contiguous bones affected while others are spared may be due to longitudinal expansion of mutant osteoblasts precursors along the axis of development ( Supplementary Fig. 5) 32 . The marked variability in the pattern of bone lesions between patients is likely related to mutations arising during different points during development, with earlier mutations resulting in an increased burden of disease. Further studies will delineate the roles of specific MAP2K1 mutations in different tissues at various developmental stages.
Melorheostosis can now be viewed as a genetically heterogeneous group of dysostoses. The melorheostosis seen in patients who have "spotted bone" osteopoikilosis is associated with a mutation in germline LEMD3 6 . Whyte et al. have shown that additional somatic mutations in KRAS (and possibly other genes) may play a role in the development of "dripping candle wax" bone lesions on radiographs in some patients in an osteopoikilosis kindred 8 . Patients who have sporadic melorheostosis are a distinct group and do not have germline 6,33 , or as we now show, somatic mutations in LEMD3. It remains to be seen if the radiographic similarities in melorheostosis caused by LEMD3 mutations and MAP2K1 mutations are due to involvement of a common pathway with activation of SMADs, or are simply a similar clinical outcome.
Our finding that three of four patients demonstrated mutations in the overlying skin raises the possibility that testing of skin for MAP2K1 mutations may be a diagnostic test for melorheostosis. However, if the skin sample is negative, bone biopsy may still be necessary to detect the mutation, especially in patients with a low mutant allele frequency. Although the causative gene is an oncogene, the clinical data from many years is reassuring to patients that malignant degeneration is rare 2,4 .
Development of a diagnostic test is especially valuable because identification of somatic MAP2K1 mutations in melorheostosis raises the possibility of inhibiting MEK1 to treat melorheostotic bone lesions. Several MEK1 inhibitors have already been developed for malignancies with MAP2K1 mutations identical to those we identified in melorheostosis 34 . Because melorheostosis is one of few osteosclerotic bone-forming diseases, study of the MEK1-ERK1/2 pathway may lead to insights relevant to diseases where bone formation is impaired, such as delayed fracture healing and osteoporosis.

Methods
Patient cohort and bone sampling. We recruited 20 unrelated adults with a diagnosis of melorheostosis to a NIAMS Institutional Review Board-approved protocol (NCT02504879). Diagnoses were confirmed by radiographs and increased uptake on 18 F-NaF bone scan (Fig. 1). Patients were excluded if pregnant or lactating, if under the age of 18, if actively infected or if unable to provide informed consent. Study subjects were not compensated for their participation. Fifteen patients consented to bone biopsies of affected and contralateral unaffected bone (see Supplementary Table 1

for clinical characteristics). Patients consented to publication of photographs.
Patients underwent open surgical biopsy of the affected bone under fluoroscopic guidance. The site of the biopsy of affected bone was chosen to minimize risk. Thus the most superficial area of melorheostotic bone was chosen. Under general anesthesia, the limb was exsanguinated and a tourniquet elevated (in all but one patient with proximal femur lesion). The bone was exposed and any thickened periosteum divided. Fluoroscopy was used to confirm the site of affected bone. Four drill holes were used to create a rectangle, and an osteotome used to connect the drill holes. The rectangular piece of cortical bone was then removed and sectioned in the operating room. The sections were placed immediately in cell culture media (for DNA extraction and osteoblast culture) or 70% ethanol for histology. Hemostasis was obtained with electrocautery and the wound closed in layers. A contralateral control sample of bone with no radiographic evidence of melorheostosis was also obtained.
For some patients, such as Melo-2 and Melo-18, abundant melorheostotic bone was available for evaluation. For other patients, only a small amount of affected bone could be safely sampled; priority was given to bone for DNA extraction followed by osteoblast culture over other studies. In all patients, the smallest informative sample of normal bone was taken to minimize morbidity. Thus, normal bone was often a limiting reagent. See Supplementary Table 2 for details on which patient samples underwent particular experiments.
Additional tissue sampling. In our initial study design, we were concerned that surgical biopsy of melorheostosis (which had never been done before) could lead to worsening of the disease, infection, or other complication. Therefore, only the smallest amount of usable bone was collected. After a shared mutation was identified and the absence of surgical complications assured, four patients returning for a follow-up visit underwent 6 mm punch biopsy of skin overlying an area of affected bone, and a control biopsy of unaffected skin from the contralateral side. Seven of the eight patients with MAP2K1 mutations agreed to have blood tested for mutations. DNA from peripheral blood leukocytes was tested by ddPCR for presence of the mutation.
Genomic DNA extraction and mutation detection in bone tissue. Soft tissues and cartilage were manually dissected from affected and unaffected bone of each patient. Following repeated washes in phosphate-buffered saline (PBS) containing no Ca 2+ or Mg 2+ to remove marrow, mineralized tissues were minced using a rongeur and spring scissors in a sterile conical bottom glass tube (Roboz, Gaithersburg, MD).
Genomic DNA was isolated from minced bone chips following the manufacturer's instructions (Gentra Puregene, Qiagen). In brief, minced bone samples were resuspended in cell lysis solution and digested with 0.2 mg/ml proteinase K at 55°C for 12 h. Digested proteins were precipitated and removed by centrifugation at 14,000×g for 3 min. The cleared supernatant was transferred to a NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-03720-z ARTICLE NATURE COMMUNICATIONS | (2018) 9:1390 | DOI: 10.1038/s41467-018-03720-z | www.nature.com/naturecommunications new tube and genomic DNA was precipitated with isopropanol and collected by centrifugation at 14,000×g for 5 min. The DNA pellet was washed once with 70% ethanol and air-dried at room temperature for 5 min. DNA was resuspended in DNA hydration buffer (10 mM Tris, 1 mM EDTA, pH 8.0) and its concentration was measured using a NanoDrop spectrophotometer (ThermoFischer Scientific). Genomic DNA samples underwent deep Whole Exome Sequencing (Otogenetics), with 107-270× coverage for each sample. Most exome targets (82%) had at least 50× coverage in each sample. Agilent SureSelect Human 51 Mbp All Exon kit was used for exome capture. Paired-end sequencing was performed on Illumina HiSeq sequencers. Sequencing reads were returned in FASTQ format and were aligned to Human Reference Genome Build 37 using Burrows-Wheeler alignment Tool. The standard PICARD-GATK pipeline was used to remove duplicate reads, refine alignment around indels, and recalibrate base quality scores. The resulting BAM files (one per sample) served as inputs to somatic mutation callers (Fig. 2a). Two somatic mutation callers (muTect v1.1.7 and Strelka v1.0.14 35 ) were used, with default parameters for WES data. Somatic variants called by either method were annotated with functional impact and population frequency using ANNOVAR. We filtered somatic variants for those that were not present in the unaffected bone, caused protein or ncRNA sequence changes or splicing changes, had less than 1% frequency in ExAC and 1000 Genomes Project databases, were not in a duplicated genomic region, and were called by both muTect and Strelka. Using this method, 284 putative somatic variants were identified in 15 affected bone samples (Supplementary Data 1 and Supplementary Fig. 2).
Exome sequencing was complemented by high-coverage amplicon sequencing with an average sequencing coverage of more than 10,000 × per sample. A targeted NGS amplicon library prep kit from Swift Biosciences (56 G Oncology) was used to screen patients for MAP2K1 somatic mutations. The panel targets cancerassociated mutations in 56 genes including 5 mutation hotspot regions in MAP2K1 (exons 2, 3, 6, 7, 11), as well as the RAS codon 12, 13 and 61 hotspot mutations. The amplicon libraries were sequenced on Illumina MiSeq and variants were called with LoFreq 2.1.2.
To validate and quantitate the allele frequency of somatic mutations in other tissues and cells cultured ex vivo, we utilized ddPCR (Fig. 2c). The approach determines the abundance of each allele using allele-specific hydrolysis probes (Supplementary Table 7) complementary to mutant allele (FAM-labeled) and wildtype allele (HEX-labeled) and is based on partitioning of each sample into >10,000 nanoliter-sized droplets, in which PCR is carried out. gBlock oligos of 200 bp bearing the mutation of interest were designed by Integrated DNA Technologies, Inc. (Coralville, IA) and used as mutation-positive controls. Reactions were performed on the QX200 Droplet Digital PCR System (Bio-Rad), using standard PCR cycling conditions based on manufacturer's instructions (10-min enzyme activation at 95°C, followed by 40 cycles of 30-s 94°C denaturation and 1-min 55°C annealing/extension, finishing with a 10-min enzyme deactivation at 98°C). QuantaSoft software was used to quantitate the concentration of each allele as the number of copies per microliter, using a Poisson distribution model to account for the number of target copies per individual droplet. Fractional abundance of mutant allele was calculated as the concentration (copies per microliter) obtained in FAM channel divided by sum of concentrations of FAM and HEX channels.
Bone histology and histomorphometry. Bone samples from a subset of 6 patients with MAP2K1 mutations were embedded in polymethylmethacrylate and prepared for histomorphometric analyses using standard procedures 36 . For histological examinations, thin sections (4 µm) were cut from the tissue blocks with a hard tissue microtome (Leica SM2500, Leica Microsystems Nußloch GmbH), deplasticized with 2-methoxyethyl-acetate and stained with modified Goldner's Trichrome 36 (Fig. 3d). A light microscope equipped with a video camera (Zeiss Axiophot, Zeiss AxioCam, Oberkochen, Germany) was used to obtain digital images of the sections that were analyzed using NIH ImageJ software (version 1.63) 37 . Bone histomorphometric analyses were performed on 4 randomly chosen areas throughout each bone section (Fig. 3c) Subsequently, the residual blocks were prepared for backscattered electron imaging by grinding and polishing to obtain plane parallel surfaces, then carbon coated (Fig. 3a, b). The entire cross-sectioned bone sample area was imaged with a spatial resolution of 1.8 µm per pixel using a field emission scanning electron microscope (FESEM) (Zeiss Supra 40, Oberkochen, Germany) equipped with a four-quadrant semiconductor backscatter electron detector. The FESEM was operated with electron energy of 20 keV. The gray levels reflecting the mineral/ calcium content were calibrated by the material contrast of pure carbon and aluminum.
Statistical evaluation was performed with GraphPad Prism 6.0 (GraphPad Software, Inc., La Jolla, CA, USA). Comparison of histomorphometric indices were based on paired t-test comparing affected vs. unaffected tissue in each patient. Statistical significance was considered as p < 0.05.
Cell culture. Primary osteoblasts were grown from freshly minced bone chips according to the method of Robey and Termine 38 . Bone chips (~1 mm) were incubated in alpha-minimal essential medium (α-MEM, Life Technologies) containing type II collagenase (200 U/ml, Worthington Biochemical Corp., Lakewood, NJ) for 2 h at 37°C with gentle rocking. Bone chips were allowed to settle briefly, then collagenase solution was removed and washed twice with PBS. Bone chips were placed in T75 tissue culture flasks with α-MEM supplemented with 10% fetal bovine serum (Gemini Bio-Products, West Sacramento, CA) and antibiotics (penicillin (100 U/ml) and streptomycin (100 µg/ml), Thermo Fisher Scientific Inc.). Bone chips were incubated for 2-3 weeks at 37°C/8% CO 2 , refreshing culture media every 3 days. Osteoblasts of passages 1-4 were used in this study.
Genomic DNA was extracted from cultured osteoblasts using the Gentra Puregene kit (Qiagen), as described above. The proportion of mutant MAP2K1 allele in DNA from multiple cell passages and revivals was determined by ddPCR as above.
Cell proliferation assay by live-cell imaging. Osteoblasts were plated at the indicated densities in 96-well tissue culture dish and imaged at 2 h intervals for 120 h using the IncuCyte ZOOM Kinetic Imaging System (Essen Bioscience) at 10× magnification (Fig. 5c). During the live-cell imaging, cells were incubated in 37°C/ 8% CO 2 . Percent confluence was calculated from the percentage of the well area occupied by cells. Doubling time was calculated using GraphPad Prism 6.0 using non-linear regression. To analyze comparable populations of cells between groups, only wells with confluency at the first time-point of image acquisition within one half standard deviation of the mean were included in the analysis. Statistical comparison was done by two-way ANOVA.
Gene expression in osteoblasts by real-time quantitative PCR. Total RNA was extracted from osteoblasts using the RNeasy mini kit (Qiagen) following the manufacturer's instruction. RNA concentration was measured by the NanoDrop spectrophotometer. Synthesis of cDNA was performed using 1 µg of total RNA and the High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scientific), following the manufacturer's instructions. Comparative real-time PCR was performed in triplicate with an Applied Biosystems Prism 7500 Fast Sequence Detection System using TaqMan universal PCR master mix, according to the manufacture's protocol (Applied Biosystems Inc., Foster City, CA) (Fig. 5a, f, g). The TaqMan probes and primers were purchased from Applied Biosystems: Human RUNX2 (Hs01047973_m1), COL1A1 (Hs00164004_m1), ALPL (Hs01029144_m1), RANKL (TNFSF11, Hs00243522_m1), OPG (TNFRSF11B, Hs00900358_m1). Human GAPDH (Hs02786624_g1) and TBP (Hs00427620_m1) were used as endogenous controls for normalization. Levels of MAP2K1, RUNX2, ALPL and COL1A1, RANKL (TNFSF11), and OPG (TNFRSF11B) transcripts were determined using the 7500 Fast System SDS software version 1.3.1 (Applied Biosystems). Relative expression was calculated using the comparative ΔΔCt method.
Immunohistochemistry for p-ERK1/2. Fresh bone tissue was fixed in 10% neutral buffered formalin, decalcified by formic acid, and embedded in paraffin. Thin sections (5 µm) were mounted onto tape adhesive slides. Following deparaffinization and rehydration with ethanol, endogenous peroxidase was blocked using 3% hydrogen peroxide/distilled water. Antigen retrieval was performed using a citrate buffer (BioGenex) for 10 s at 100°C, then slides were cooled and rinsed in distilled water and Tris-buffered saline/Tween-20 (TBS-T). After blocking in 5% normal goat serum (Vector Labs), sections were incubated with anti-phospho-p44/42 MAPK ERK1/2 rabbit monoclonal antibody (Thr202/Tyr204; diluted 1:200, Cell Signaling Technology, USA #4376) overnight at 4°C or monoclonal rabbit IgG isotype control. Subsequently, biotinylated goat anti-rabbit IgG (1:100, Vector Labs) was applied for 30 min, before slides were incubated with ABC Elite Standard (Vector Labs) for 30 min. p-ERK1/2 was visualized by a 5-min incubation in DAB Reagent (Sigma) and subsequent hematoxylin counterstain. Slides were dehydrated, cleared in xylene and mounted with coverslips with Permount. Paraffin embedded colon cancer cells (SW48, carrying a p.Q56P mutation in MAP2K1) was used as positive control 13 . Normal bone from patients undergoing surgical resection for unrelated indications served as negative control (Fig. 4a).
To visualize mineralization in osteoblast culture, cells were washed once with PBS and fixed with 4% paraformaldehyde for 30 min at room temperature. Fixed cells were washed three times with PBS and stained with 2% Alizarin Red S solution (pH 4.2) for 30 min at room temperature. Excess Alizarin Red S stain was removed with three distilled water washes. Cells were air-dried in the dark before imaging (Fig. 5e). During osteoblast mineralization in osteogenic medium, total RNA was collected at 2 weeks' post-induction and used to examine gene expression by real-time quantitative PCR analysis.
Data availabilty. The data that support the findings of this study are available from the corresponding author upon reasonable request.