Somatic SMAD3-activating mutations cause melorheostosis by up-regulating the TGF-β/SMAD pathway

Melorheostosis is a rare, sporadic, sclerotic bone dysostosis with asymmetric distribution (Online Mendelian Inheritance in Man #155950). The bone lesions in this progressive hyperostotic disorder are associated with pain, functional impairment, joint contracture, and deformity. First described by Leri and Joanny (1922), the etiology of isolated sporadic melorheostosis was long elusive until recent findings of somatic mosaicism for the genes in the RAS-MEK-ERK pathway. Postzygotic mosaicism for a mutated KRAS proto-oncogene GTPase (p.Q61H) was identified in dermatoses from a boy with polyostotic melorheostosis, who also has familial osteopoikilosis (OPK; Whyte et al., 2017). The same KRAS mutation was found in the cervical lymphatic malformation and hyperpigmented skin from a girl with polyostotic melorheostosis and vascular anomalies (Seidel, V., E. Guillén, V.M. Martínez-Glez, Á.M. Lancharro Zapata, F. Ballesteros Tejerizo, V.A. Parra Blanco, A. García Martín, A. Salcedo Posadas, A., and M. Campos Domínguez. 2018. European Society of Human Genetics Conference. Abstr. P11.051C). The definitive connection of somatic mutations to bone lesions of patients with the classic form of melorheostosis, with a radiographical “dripping candle wax” appearance, involved mutations in MAP2K1. Half of 15 melorheostosis patients biopsied in our study had somatic mosaicism for activating mutations in MAP2K1 in affected, but not unaffected, bone and skin (Fratzl-Zelman et al., 2019; Jha et al., 2019; Kang et al., 2018).

Melorheostosis can also coexist with OPK (“spotted bone”) or Buschke–Ollendorff syndrome (BOS; dermatoosteopoikilosis), a benign autosomal dominant sclerosing bone dysplasia with symmetric distribution. Germline heterozygous loss-of-function mutations in LEMD3 (Hellemans et al., 2006; Hellemans et al., 2004; Mumm et al., 2007) cause OPK and BOS. LEMD3 encodes MAN1, an integral protein of the inner nuclear membrane, which interacts directly with phosphatase PPM1A and receptor-regulated SMADs (R-SMADs; SMAD2/3 or SMAD1/5) to antagonize TGF-β and bone morphogenetic protein (BMP) signaling (Bourgeois et al., 2013; Lin et al., 2006). Although somatic mutations in LEMD3 have not been reported in isolated sporadic melorheostosis (Mumm et al., 2007), as melorheostosis often coexists with OPK in the same individual or within a family with LEMD3 mutations, the TGF-β/SMAD pathway may contribute to melorheostosis pathogenesis.

TGF-β subfamily ligands transduce cellular signals through cell surface serine/threonine kinase receptors to intracellular SMADs. The TGFβRI kinase phosphorylates serine residues within the conserved SSxS motif at the carboxy-terminus of R-SMADs, SMAD2 and SMAD3. Phosphorylation of the carboxy-terminus of R-SMADs mediates oligomerization between R-SMADs and the coSMAD SMAD4. The complex translocates into the nucleus, where it functions to activate or repress expression of specific sets of TGF-β/SMAD target genes, depending on spatial and temporal tissue contexts (David and Massagué, 2018).

The TGF-β/SMAD pathway is crucial for skeletal embryonic development and postnatal homeostasis (Wu et al., 2016). Dysregulation of the TGF-β signaling pathway is associated with a spectrum of osseous defects as seen in several dominant genetic disorders. In Marfan syndrome, excessive SMAD-independent TGF-β signaling via ERK/MAPK activation induces aortic aneurysms (Holm et al., 2011), as well as tall stature and scoliosis. In Loeys-Dietz syndrome (LDS), skeletal involvement is caused by pathogenic variants of TGF-β signaling pathway components, including SMAD3 (Loeys and Dietz, 1993; Schepers et al., 2018). Activating germline mutations in the TGFB1 gene and uncontrolled activation of the TGF-β pathway cause Camurati-Engelmann disease, a progressive diaphyseal hyperostosis (Janssens et al., 2000; Kinoshita et al., 2000).

We noted distinctive features in the clinical, radiographical, and bone phenotypes between MAP2K1 mutation–positive versus –negative patients in the National Institutes of Health (NIH) cohort of 15 biopsied melorheostosis individuals (Jha et al., 2019). We hypothesized that different melorheostosis radiographical patterns may be associated with distinct molecular and cellular mechanisms caused by postzygotic mutations in different genes. In our cohort, some MAP2K1 mutation–negative patients showed an endosteal hyperostosis with asymmetric distribution. Here, we report four melorheostosis patients with an endosteal radiographical pattern whose affected, but not unaffected, bone samples contain novel somatic missense mutations causing substitution of serine 264 in the MH2 domain of SMAD3. The mutations increase SMAD3 activity and drive increased osteogenesis. Functional studies and RNA sequencing (RNA-Seq)–based transcriptome analysis implicate the TGF-β/SMAD pathway in the pathogenesis of endosteal melorheostosis and connect this form of melorheostosis to the TGF-β pathway activation that occurs with LEMD3 loss-of-function mutations causing OPK and BOS.

Seven patients in a cohort of 15 melorheostosis patients who underwent paired bone biopsies had an endosteal pattern of melorheostosis (Fig. 1, A–C; and Fig. S1). We are using the new term “endosteal melorheostosis” to describe the radiographical pattern seen in our patients, although a previous report used the term “osteoma-like” to categorize a similar radiographical pattern seen in 7 of 23 cases (Freyschmidt, 2001). Distinct from osteomas—which are small, have a radiolucent nidus, and tend to expand beyond the cortices—the abnormal bone tends to fill the intramedullary and intracortical space, so endosteal melorheostosis is a more descriptive term. This pattern is defined by hyperostosis that does not extend outside the bone cortical margin and has eccentric location at or in the bone. Lesions are oriented in the long axis of the involved bone and have a diameter ≥5 cm. One or more bones are involved, or if only one bone is involved, there is a circumscribed skin or subcutaneous fibrosis above the involved skeleton (Table S1; Freyschmidt, 2001). The lesions noted on radiographs corresponded to increased uptake on positron emission tomography/computed tomography using sodium fluoride tracer (Fig. 1 D), suggestive of increased bone turnover and vascularity of the lesions. This endosteal pattern of melorheostosis contrasts with the classic “dripping candle wax” radiographical appearance of melorheostosis seen in patients with MAP2K1 mutation–positive melorheostosis (Kang et al., 2018).

In four of the seven patients (Table 1), somatic mutations were detected in DNA from affected bone at the codon for the same residue in SMAD3 (Fig. 2). The remaining three patients with an endosteal pattern of melorheostosis were negative for mutations in either MAP2K1 or SMAD3. Relevant clinical findings from patients with SMAD3 variants are summarized (Table S1). Disease duration ranged from 6 to 43 yr; the two patients with higher variant allele frequency (VAF) in bone tissue (Melo-12 and Melo-17) presented in childhood. No cutaneous changes were noted in any of the four patients. Melo-12 also had a history of renal artery stenosis ipsilateral to her melorheostosis skeletal lesion. Due to uncontrolled hypertension, she required revascularization and bypass at 8 yr of age. In addition, Melo-12 incurred multiple stress fractures in her foot, the majority of which occurred on the side affected with melorheostosis. There was no history of prior fractures or vascular problems in the remaining three patients.

As germline mutations in SMAD3 have been described in LDS, characterized by vascular (tortuosity, aneurysms, and/or dissections), skeletal (pectus excavatum, scoliosis, or joint laxity), craniofacial (widely spaced eyes, strabismus, and cleft palate), and cutaneous findings (velvety and translucent skin, easy bruising, and dystrophic scars; Loeys and Dietz, 1993), the four patients were evaluated for features of LDS. No significant LDS findings were noted on physical exam or on computed tomography of the chest and abdomen with contrast (Table S1). Echocardiograms revealed normal aortic root and ascending aorta dimensions.

Biochemical markers of bone turnover were normal, except one patient (Melo-11) had decreased urinary N-telopeptide, a marker of bone resorption, and osteopenia. Melo-12 had osteopenia with uncertain etiology.

After applying the filtering strategy shown in Fig. 2 A, 284 potential somatic mutations were present in affected, but not unaffected, bone in 15 patients with melorheostosis. Among them was a novel SMAD3 p.S264Y (c.791C>A) variant at 24% VAF in patient Melo-12. While this variant was present in only one patient, the SMAD3 gene was prioritized based on the relatively high VAF and biological relevance of the TGF-β pathway to bone disease.

After validation of the initial SMAD3 p.S264Y variant by droplet digital PCR (ddPCR), we screened several MAP2K1 mutation–negative patients for this mutation by ddPCR and identified a second individual (Melo-8) with a SMAD3 p.S264Y mutation at 2.3% VAF. Prompted by these findings, we developed a custom amplicon-based targeted sequencing assay, which included all exons of SMAD3 as well as other genes potentially related to melorheostosis (KRAS, MAP2K1, MAP2K2, and MAP3K3). We sequenced the remaining MAP2K1 mutation–negative patients and identified two more patients (Melo-11 and Melo-17) with SMAD3 variants at the same residue: one, p.S264Y (c.791C>A), at 1.7% and one, p.S264F (c.791C>T), at 8.2% VAF (Fig. 2 B and Fig. S2, A and B). Both variants were confirmed by ddPCR (Table 1). Although skin overlying affected bone appeared normal in our patients, mosaicism for the SMAD3 mutation was also found in skin overlying affected bone in all four patients, at much lower allele frequencies than in bone (skin VAF ≈10% of bone VAF), but not in skin from the contralateral extremity. DNA from peripheral blood leukocytes tested negative for these mutations in all four patients (Table 1).

The SMAD3 serine264 residue is highly conserved across species (van de Laar et al., 2011). Both SMAD3 mutations (p.S264Y and p.S264F) are absent from population databases such as gnomAD and 1000G. Pathogenicity prediction software tools (SIFT, PolyPhen2, MutationTaster, MutationAssessor, and CADD [Combined Annotation Dependent Depletion]) assigned both mutations as damaging or disease-causing. The CADD scores for the two mutations are 32 for S264Y and 29.7 for S264F.

Histological evaluation revealed that the bone biopsy samples from melorheostotic lesions consisted of a mixture of osteonal-like remodeled bone, with concentric lamellae around vascular channels, and areas of aligned parallel lamellar bone (Fig. 3 A). Affected bone from Melo-12 showed mostly multiple layers of long parallel lamellae, indicating intense periosteal bone apposition associated with the endosteal radiographical pattern of melorheostosis. As in unaffected bone, osteonal bone surfaces in affected bone were characterized by low osteoid formation, sparse eroded surfaces, and absence of osteoblasts (Fig. 3 A and Table S2). Unaffected bone samples contained mostly trabecular bone (Table S2). Histomorphometric indices of bone formation and bone resorption were not significantly altered between affected and unaffected bone (Table S2). These results contrast markedly with our previous findings in melorheostotic lesions caused by somatic-activating mutations in MAP2K1 that showed intense bone remodeling and increased osteoid (Jha et al., 2019; Kang et al., 2018). In comparison, affected bone from patients with the SMAD3 mutations has significantly lower osteoid thickness (−62%, P = 0.003), osteoid surface per bone surface (−92%, P < 0.0001), and eroded surface per bone surface (−84%, P = 0.0016) than affected bone from patients with MAP2K1 mutations. Active bone-forming osteoblasts were not viewed in the selected areas; osteoclast indices are not statistically different between these groups.

Bone mineralization density distribution (BMDD) of biopsy samples from SMAD3 mutation–positive melorheostosis patients are shown in Fig. 3 B. Compared with unaffected bone, the BMDD curves from affected bone shifted toward higher mineral content in Melo-11 and Melo-12, while that of Melo-8 and Melo-17 were within normal limits. Statistical analysis supported an overall increase in mineralization of affected bone with the SMAD3 mutations, compared with both paired unaffected biopsies and MAP2K1 mutation–positive affected bone (Table S2). Of the BMDD variables, only CaHigh was significantly altered in the set of biopsies with SMAD3 mutations. CaHigh was nearly doubled in SMAD3 mutation–positive melorheostotic lesions compared with unaffected bone (P = 0.0240). There was also a strong trend to higher average mineral content in affected bone (CaMean: +6.7%, P = 0.0596).

Compared with MAP2K1 mutation–positive affected bone, lesions with SMAD3 mutations were more highly mineralized, with CaMean +7.8%, P = 0.002 and CaPeak +4.0%, P = 0.0048. CaHigh values were very variable in both groups; on average, however, the value was doubled in the SMAD3 mutation–positive group (P = 0.0005). In contrast, CaLow was markedly reduced (−65%, P = 0.0012) compared with MAP2K1 mutation–positive affected bone (+89%, P = 0.009; Fratzl-Zelman et al., 2019).

The identified SMAD3 missense mutations occurred in exon 6 of the gene, substituting serine at residue 264 with either tyrosine (p.S264Y) in three patients or phenylalanine (p.S264F) in one patient. To investigate the biological consequences of the novel SMAD3 mutations, primary osteoblasts isolated from affected and unaffected bone tissues of the same individual were subjected to Western blot analysis for SMAD3 activation upon TGF-β stimulation. Cells from affected bone showed higher SMAD3 carboxy-terminal SSxS motif phosphorylation relative to their unaffected counterparts (Fig. 4 A). This was detectable in affected cell cultures containing relatively high VAF (41%). In affected cells with relatively low VAF (∼1–2%), there were marginal differences between affected and unaffected cells (Fig. S3, A and B).

To further examine the effect of the mosaic SMAD3 mutation on the canonical TGF-β/SMAD pathway, cells were stimulated with TGF-β and target gene expression was examined by real-time quantitative PCR (qPCR). In agreement with the increased level of SSxS motif phosphorylated (activated) SMAD3, expression of TGF-β/SMAD pathway target genes, such as CDKN2B, SERPINE1, COL1A1, and COL1A2, was higher in affected than in unaffected cells (Fig. 4 B). These results indicate that the SMAD3 mutation is activating for the TGF-β/SMAD pathway.

Functional activity of the SMAD3 mutations (p.S264Y and p.S264F) was further investigated by transfecting plasmid vectors encoding wild-type and mutant SMAD3 into MC3T3-E1 (Fig. 4, C–E). In the absence of TGF-β stimulation, the basal level of SMAD3 phosphorylation was significantly higher in SMAD3 S264 mutants than wild-type SMAD3. Both SMAD3 mutants showed markedly higher levels of carboxy-terminal SSxS motif phosphorylation than wild-type SMAD3 upon TGF-β stimulation (Fig. 4 C). Under these conditions, the increased phosphorylation of SMAD3 S264 mutants appears to be constitutive, as it did not respond to further stimulation by TGF-β and was resistant to inhibition of the kinase activity of TGF-β receptor (Fig. S3 C). The hydroxyl group of tyrosine residues may be subject to post-translational modification, such as phosphorylation or glycosylation. However, substitution of the Tyr264 residue by alanine (p.S264Y to p.S264A) did not abolish the increased SMAD3 phosphorylation (Fig. S3 D), implying that loss of serine at residue 264 is critical to SMAD3 activation.

We noted that when SMAD3 was overexpressed, SMAD3 mutant had lower total SMAD3 than wild-type SMAD3 on Western blots, indicative of a higher rate of SMAD3 turnover or lower stability of mutant than wild-type SMAD3 (Fig. 4 C). Since activated SMAD3 undergoes rapid turnover by proteasome-mediated degradation as a part of the mechanism that regulates duration of TGF-β/SMAD signaling (Massagué, 2000), this indicates that levels of phosphorylated SMAD3 are increased in overexpressed mutants even in the context of lower total SMAD3.

Nuclear SMAD complexes (e.g., SMAD2/3/4) partner with transcription coactivators or corepressors to execute their TGF-β signaling responses (David and Massagué, 2018). Due to the mosaic nature of the SMAD3 mutations in cells of melorheostotic bones, we relied on a plasmid vector–based transfection system to assess nuclear translocation of the SMAD3 mutants. SMAD3 p.S264Y was expressed in MC3T3-E1, and nuclear translocation of SMAD3 was assessed by immunofluorescence. Consistent with its increased phosphorylation, SMAD3 p.S264Y displayed increased nuclear localization compared with wild-type SMAD3, even without TGF-β stimulation (Fig. 4 D). Further increase in nuclear translocation upon TGF-β stimulation appears to represent activated endogenous SMAD3.

The increased phosphorylation and nuclear translocation of mutant SMAD3 resulted in the expected increase in transcription of TGF-β/SMAD pathway target genes Ctgf, Serpine1, and Vegfa (Fig. 4 E). Up-regulation of target gene expression in SMAD3 mutant cells was significant, even without TGF-β stimulation. TGF-β stimulation further increased the target gene expression, but these enhanced levels likely represent the response of the endogenous wild-type SMAD3 rather than of plasmid vector–expressed SMAD3, since the extent of up-regulation was similar between wild-type and mutant SMAD3. These results suggest that the SMAD3 mutations (p.S264Y and p.S264F) are gain-of-function mutations augmenting signaling of the TGF-β/SMAD pathway.

SMAD3 plays a critical role in TGF-β–mediated regulation of cell growth and osteoblast differentiation (Sowa et al., 2002). Temporal cell cycle arrest in the G₁ phase is a prerequisite for cell differentiation (Kaji et al., 2006). Overexpression of SMAD3 has both a cytostatic effect and a stimulatory effect on bone matrix protein expression (Kaji et al., 2006). To examine the effect of the mosaic-activating SMAD3 mutations on cell growth, the cell proliferation rate of osteoblasts isolated from affected and unaffected bone of a melorheostosis patient was assessed by live-cell imaging. Osteoblast cultures from affected bone with relatively high VAF (37%) of SMAD3 p.S246Y displayed a significantly decreased cell growth rate compared with unaffected cells (Fig. 4 F). Proliferation of affected cells with very low VAF (1–2%) of SMAD3 p.S264Y was comparable to that of unaffected cells (Fig. 4 G).

To examine the impact of the mosaic SMAD3 mutation on osteogenesis in vitro, cells from affected and unaffected bone tissue of a melorheostosis patient (Melo-11, SMAD3 p.S264Y, VAF 41% at passage 3) were cultured in conventional osteogenic media. Osteogenic differentiation and ECM mineralization were significantly increased in affected compared with unaffected cells (Fig. 5 A). A concomitant increase in osteogenesis was confirmed by significantly higher expression of osteogenic marker genes in cells from affected bone (Fig. 5 B).

To further illuminate cellular pathways and biological processes influenced by the mosaic SMAD3 mutation, transcriptome profiles were acquired by RNA-Seq on cells harvested at 1-wk intervals during osteogenic differentiation in culture. Totals of 261, 172, and 250 genes were found to be differentially expressed with more than twofold changes (false discovery rate [FDR] < 1e-20) at 0-, 1-, and 2-wk time points, respectively (Fig. S4 A). Osteogenic markers and transcription factors, represented by SP7 (Osterix), and other skeletal-related genes (Ho et al., 2000) were among the up-regulated genes in affected cells (Fig. 5, C and D). Examination of the differentially expressed genes (DEGs) by Ingenuity Pathway Analysis identified β-catenin (CTNNB1), the signal transducer of the canonical Wnt pathway, as one of the top upstream regulators (KMT2A, POLG, KAT6A, CTNNB1, and HOXB3) of the DEGs, further accounting for the increased osteogenesis in affected cells (Table S3).

Since genes sharing similar expression patterns during osteogenic differentiation may also be involved in similar cellular processes, coexpression clustering of the DEGs identified 10 clusters, each of which displayed distinct expression trends over the time points (Fig. S4 B and Data S1). To gain insight into putative cellular pathways that were changed by the mosaic SMAD3 mutation, gene ontology and pathway enrichment analyses were performed on the DEGs residing in each coexpression cluster (Data S1 and Data S3). Further exhibiting functionality of the activating SMAD3 mutation, the TGF-β pathway was enriched in cluster 2, with gene transcript levels increased at baseline and subsequently (Fig. 5 E). Up-regulation of the TGF-β pathway is not attributable to LEMD3 (MAN1) action, as our RNA-Seq data indicated that transcript levels of LEMD3 and PPM1A, whose protein products form a complex to dephosphorylate R-SMADs and terminate TGF-β signaling, remained steady during osteoblast differentiation and were not altered by the SMAD3 mutations. Intriguingly, the interferon signaling pathway and cytokine signaling in the immune system were also enriched in cluster 2, with significantly higher expression in affected than in unaffected cells. The set of 21 genes whose expression was previously validated in systemic lupus erythematosus for interferon-inducible gene signature was found to have high- to medium-fold increase in affected cells of Melo-11 (Table S4; Yao et al., 2009). The genes in cluster 3, with increased transcripts in affected cells at baseline and week 1, were associated with adipogenesis and inflammatory response, in addition to ossification (Fig. S4 B and Data S3). Enrichment of genes for ECM organization was found in clusters 4, 5, 6, and 8, in which expression is down-regulated in affected cells. Genes related to translation and ribosome biogenesis were over-represented in cluster 7. Given that ribosomal proteins have a significant role in ribosome biogenesis, protein synthesis, and cell growth (Zimmermann, 2003), down-regulation of ribosomal protein subunit genes contributes to the reduced cell growth in affected cells (Fig. S4 C and Data S3).

Surprisingly, BMP2 supplementation in osteogenic culture media strikingly dampened the mutant SMAD3’s stimulatory effect on osteogenesis (Fig. 6 A). Concordantly, osteogenic marker genes elevated in untreated affected cells (Fig. 5 B) were significantly down-regulated with BMP2-treated compared with unaffected cells (Fig. 6 B). These results suggest that activated SMAD3 inhibited cellular signaling events activated by BMP2—including those associated with osteogenic differentiation and ECM mineralization.

Transcriptome profiling by RNA-Seq on Melo-11’s affected and unaffected cells, cultured in the presence of BMP2 and compared pairwise at each time point, identified 269, 206, and 624 DEGs with more than twofold changes at 0, 1, and 2 wk, respectively (FDR < 1e-20; Fig. S4 D). Supporting the antagonistic effect of BMP2, osteogenic marker genes, as exemplified by SP7 in MA plot, were significantly down-regulated in affected cells in the presence of BMP2 (Fig. 6, C and D).

Analysis of coexpressed genes identified nine clusters in the BMP2(+) condition (Data S2 and Data S4). Gene ontology analysis and pathway enrichment of each cluster revealed that genes in cluster 1 were associated with the TGF-β pathway, ECM organization, and ossification (Data S2 and Data S4). In agreement with the inhibitory effect of the SMAD3 mutation on osteogenic differentiation in the presence of BMP2, the transcriptional response of genes in cluster 1 was noticeably down-regulated at week 2 (Fig. 6 E). Genes in cluster 2 and cluster 8, with up-regulation at week 1 or at all time points, respectively, were implicated in interferon signaling and immune response (Fig. S4 E and Data S4). Genes for ribosome biogenesis and protein translation were over-represented in cluster 5, containing genes down-regulated in affected cells at baseline and subsequently. Interestingly, the transcriptional response of genes for interferon response (cluster 2, both BMP2[−] and BMP2[+]) and ribosomal proteins (cluster 7 in BMP2[−]; cluster 5 in BMP2[+]) appeared to be highly up-regulated and down-regulated, respectively, in affected cells regardless of BMP2 (Fig. S4, B and E; and Data S3 and Data S4). These results showed that BMP2 pathway activation in cells with mosaic-activating SMAD3 mutation decreased osteogenic differentiation and ECM mineralization (graphical summaries in Fig. 7). In contrast, the activating MAP2K1 mutations, which we identified in our previous report (Kang et al., 2018), significantly inhibited osteoblast differentiation and ECM mineralization regardless of the presence of BMP2 (Fig. S5). Along with the histology of melorheostotic bone with SMAD3 mutations, these results indicate that melorheostosis caused by SMAD3 mutations is clinically, pathologically, and mechanistically distinct from MAP2K1 mutation–positive melorheostosis.

Melorheostosis is an asymmetrically distributed sclerosing bone dysostosis caused by postzygotic mosaicism in the osteoblast lineage. In our previous report, we identified activating mutations in MAP2K1 as causative for classic melorheostosis with a radiographical “dripping candle wax” appearance (Kang et al., 2018). While this form is best known, its prevalence is lower than that of other radiographical patterns of melorheostosis (Freyschmidt, 2001). In this study, we report somatic mosaicism in SMAD3 for two novel mutations recurring at the same residue (p.S264Y or p.S264F) in four individuals whose melorheostosis has an endosteal radiographical pattern. The serine 264 residue is located in the region within the MH2 domain that is evolutionarily conserved (van de Laar et al., 2011). These patients had tested negative for mutations in MAP2K1, KRAS, or LEMD3.

The mosaic SMAD3 p.S264 substitutions were demonstrated to cause gain of function, with a constitutive increase in TGF-β signaling. Nuclear translocation of mutant SMAD3 was increased, and expression of TGF-β target genes was enhanced. Further, mosaicism for the SMAD3 mutation inhibited the proliferation of affected cells, consistent with the known action of TGF-β/SMAD signaling to modulate cell proliferation and growth in a cell context–dependent manner (Muñoz et al., 2008). TGF-β stimulates cell growth in fibroblasts and epithelial cells by positively regulating the mTOR pathway (Rahimi et al., 2009; Wu and Derynck, 2009). The growth inhibitory function of TGF-β was observed in muscle atrophy and hematopoietic stem cells, in which TGF-β/SMAD signaling counterbalanced the Akt/mTOR pathway (Chabanon et al., 2008; Sartori et al., 2009). Likewise, our RNA-Seq data revealed that osteoblasts suppressed expression of genes related to growth pathways regardless of BMP2 stimulation. The proximal cause of decreased cell growth in melorheostotic cells may be down-regulation of multiple ribosomal protein genes, with their significant role in protein synthesis (Maguire and Zimmermann, 2001; Naora, 1999).

In melorheostotic cells, the SMAD3 mutation promoted osteogenic differentiation and ECM mineralization. This was consistent with the findings of generally increased mineralization in affected tissue determined by quantitative back-scattered electron imaging. Specifically, CaHigh reflects the percentage of bone mineralized above 25.30 weight percent calcium, and thus older tissue age; in SMAD3 mutation–affected versus –unaffected bone, the threefold increase in CaHigh is thus indicative of slow remodeling rates, leading to increased mineral content over time.

TGF-β has bi-functional effects on osteogenesis depending on the stage of cell differentiation. TGF-β inhibits differentiation of mature osteoblasts and ECM mineralization (Kaji et al., 2006) but stimulates osteoprogenitor proliferation and early-stage osteoblast differentiation by enhancing ECM deposition and anti-apoptotic activity. The stimulated osteogenesis of SMAD3 mutation–positive melorheostosis reflected increased expression of osteogenic transcription factors (Chen et al., 2012), such as SP7 and MSX2 in affected cells. SP7 (Osterix) is the master regulator of osteoblast differentiation acting downstream of RUNX2 and regulates cell differentiation from preosteoblast to mature osteoblasts and osteocytogenesis (Nakashima et al., 2002). TGF-β signaling regulates SP7 expression through a SMAD-dependent pathway (Choi et al., 2016). MSX2, a homeobox transcription factor, is a direct transcriptional target of TGF-β and inhibits RUNX2 DNA-binding activity (Shirakabe et al., 2001); mutations in MSX2 are associated with craniosynostosis and enlarged parietal foramina (Ma et al., 1996; Wilkie et al., 2000). Our results suggest that the constitutively active form of SMAD3 in melorheostotic bone augmented the output of the TGF-β/SMAD pathway and overrode stage-specific osteogenesis regulation.

The SMAD3 serine 264 residue, which is substituted in melorheostosis, is located in the MH2 domain of SMAD3 on the L1 loop at the periphery of the protein-protein interface and forms a second network of hydrogen bonds with residues of the adjacent monomer (Fig. S2; Wu et al., 2001). There is an overrepresentation of hits in the SMAD3 MH2 domain among somatic missense mutations causing colorectal cancers (CRCs; Fleming et al., 2013). Most of these mutations reduce protein stability or hinder SMAD complex formation, undermining the role of SMAD2/3/4 as tumor suppressors in these cells (Fleming et al., 2013). Substitution in CRC at p.D258, adjacent to SMAD3 p.S264, interferes with oligomerization and decreases transcription in luciferase reporter assays, while p.S266 is critical for binding of coregulators to SMAD3. In contrast, our functional study of the melorheostosis-causing mosaic SMAD3 mutations in osteoblastic cells revealed that SMAD3 p.S264Y and p.S264F are gain-of-function substitutions that elicit stimulatory effects on osteogenesis and an inhibitory effect on cell growth.

SMAD3 germline mutations have also been demonstrated in connective tissue disorders associated with increased TGF-β signaling and skeletal, articular, and aortic findings. LDS is caused by defects in SMAD3 as well as five other genes (SMAD2, TGFB2, TGFB3, TGFBR1, or TGFBR2) in the TGF-β signaling pathway. LDS manifestations of vascular tortuosity and craniofacial and skeletal abnormalities are often associated with excess activation of TGF-β signaling, but LDS-causing genetic variants themselves are mostly loss-of-function mutations (Pezzini et al., 2012). Of the 67 SMAD3 mutations reported in LDS (Schepers et al., 2018), the majority were in the MH2 domain, but no hotspot was found. Substitutions in SMAD3 in LDS spanned exon 6, including the p.D258 and p.S266 residues associated with CRC and residues flanking the p.S264 residue substituted in melorheostosis (Fleming et al., 2013). Families with aneurysms and early onset of a distinctive osteoarthritis (originally AOS, now reclassified as LDS3) have SMAD3 p.T262I and p.R287W substitutions in the MH2 domain (van de Laar et al., 2011). As in our melorheostosis-causing SMAD3 mutations, these LDS3-causing SMAD3 mutations are predicted to interfere with heterotrimer formation. The LDS3-causing mutations were associated with increased immunohistochemical labeling of phosphorylated SMAD3, as well as upstream ligand TGF-β, in aortic tissue (van de Laar et al., 2011).

Sclerotic bone phenotypes in OPK and BOS are attributable to augmented SMAD-dependent TGF-β/BMP signaling caused by heterozygous loss-of-function mutations in LEMD3 (MAN1; Hellemans et al., 2004). Smad3 KO mice are osteopenic secondary to reduced bone formation, caused by increased apoptosis of osteoblasts and osteocytes (Borton et al., 2001). They have reduced TGF-β signaling, leading to abnormal chondrocyte development and osteoarthritis, and the majority die from ruptured aneurysms (van der Pluijm et al., 2016). Since patients with melorheostosis caused by SMAD3 p.S264 substitutions do not have findings of LDS, melorheostosis expands the type of bone abnormality associated with abnormal TGF-β signaling. Dynamic expression and stage-specific function of SMAD3 in osteoblast differentiation (Lin et al., 2019) may explain the phenotypic discrepancy between SMAD3 overexpression and melorheostosis-causing mosaic SMAD3 mutations.

Surprisingly, addition of BMP2, also a pro-osteogenic factor, to the conventional osteogenic media caused affected cells to display inhibition of osteogenic differentiation and ECM mineralization in vitro. Several factors are relevant to the mechanism of this adverse influence on osteogenesis. First, although the pro-osteogenic transcription factor RUNX2 is activated by BMP signaling, SMAD3 inhibits RUNX2 activity by interfering with its DNA binding (Soltanoff et al., 2009). Second, TGF-β stimulation results in a decreased level of SMAD1/5–SMAD4 complexes formed in response to BMP and an increase in mixed R-SMAD complexes (SMAD1/5–SMAD2/3), which bind the BMP-responsive element and exert inhibitory effects (Grönroos et al., 2012). This TGF-β action depends on SMAD3 rather than on SMAD2. Thus, the SMAD3 mutation’s constitutive activation of TGF-β may augment the inhibition of BMP-responsive gene expression and contribute to decreased osteogenesis in the presence of BMP2. It will be important to obtain further insight into the crosstalk between TGF-β/SMAD3 and BMP pathways in affected cells, since this may impact future therapeutics for SMAD3 mutation–positive melorheostosis. For example, the distinct signaling characteristics of BMP type I receptors could be exploited. BMP2 signals through the hetero-oligomerization of two type I receptor subtypes, BMPRIA (ALK3) or BMPRIB (ALK6), and the type II receptor, BMPRII. BMPRIA and BMPRIB transmit different signals to bone-derived mesenchymal progenitors, dictating whether adipogenesis or osteogenesis is induced (Chen et al., 1998). In addition, the BMPRIA expression level was ∼10 times higher than that of BMPRIB in our RNA-Seq data from melorheostotic osteoblasts. Although further studies are needed, we speculate that BMP2 muteins may be designed to signal preferentially through BMPRIA and exert an antagonistic function on the pro-osteogenic process in affected cells with the activating SMAD3 mutation.

Interestingly, up-regulation of genes regulated by the interferon signaling pathway was apparent in affected cells, regardless of BMP2 presence in culture. Osteoblasts from affected melorheostotic bone showed strong up-regulation of interferon-inducible signature genes (Yao et al., 2009) in early differentiation (Table S4). Although TGF-β generally suppresses type I interferon production, SMAD3 has been shown to enhance the effects of IRF7 and promotes interferon pathway target gene expression (Qing et al., 2004). The important physiological anabolic effects of IFN-γ signaling in bone formation were demonstrated by reduced bone volume in IFNγR1 KO mice with subsequent rescue of the ovariectomy-induced osteoporosis by IFN-γ administration (Duque et al., 2011). Given that the up-regulation of interferon signal-regulated genes was not altered by BMP2, the interferon signaling pathway does not seem to have an antagonistic crosstalk with the BMP pathway in affected cells. Whether or how the augmentation of the interferon signaling pathway and up-regulation of downstream target genes in affected cells may contribute to pathological phenotypes by inducing immune responses within melorheostotic bone lesion remains to be investigated.

Our study has some limitations and generates intriguing questions about the distinctions and similarities of melorheostosis caused by recurrent mutations in specific residues in MAP2K1 and SMAD3. Our previous report showed that activating MAP2K1 mutations increased bone remodeling within melorheostotic bone lesion. This was attributed to highly proliferative osteoblasts stimulating RANKL expression and osteoclastogenesis (Kang et al., 2018). In melorheostosis caused by SMAD3 mutations, up-regulation of SP7 may increase the increased RANKL expression in cultured affected cells. Although histology of melorheostosis bone tissues showed a minor increase in osteoclasts, the sample size (n = 4) was too small to perform statistical analysis on osteoclastogenesis-associated phenotype. Likewise, only one of four melorheostosis patients with SMAD3 mutations had a high VAF, limiting our RNA-Seq to technical replicates and a compensatory stringent threshold of FDR < 1e-20.

Finally, the clinical, bone tissue, and cellular phenotypes of SMAD3 mutation–positive melorheostosis may be compared with MAP2K1 mutation–positive melorheostosis. Patients with SMAD3 mutation–positive endosteal pattern of melorheostosis on radiographs do not present with a visible bone overgrowth. While patients with MAP2K1 mutation–positive melorheostosis were likely to have cutaneous vascular changes in skin overlying affected bone, no cutaneous changes were noted in patients with SMAD3 mutation–positive melorheostosis. The growth-promoting effects of activating MAP2K1 mutations, which we observed in osteoblasts (Kang et al., 2018), may account for the bone overgrowth outside the cortical margins in these patients, and the vascular skin phenotype may be related to the known ability of MAP2K1 mutations to promote vascular malformations (Couto et al., 2017).

Both mutations in SMAD3 led to bone growth by periosteal apposition of multilayered compact parallel lamellar bone that was subsequently remodeled, although apparently at a slower rate than in the MAP2K1 cases. This provides evidence that the so-called “periosteal reaction” is a common response in both forms of melorheostosis, as well as in other benign or malignant lesions (Bisseret et al., 2015). Although the SMAD3 mutation decreases cell proliferation, there is bone overgrowth. We speculate that this is due to focal apposition of new bone formed by cells of the periosteum, which may not carry a SMAD3 mutation. Also, several cellular distinctions are noted between melorheostosis caused by MAP2K1 mutations and by SMAD3 mutations: (a) mineralization was inhibited in MAP2K1 mutation–positive melorheostosis, with osteoid accumulation and decreased mineral content of the bone matrix, while osteogenesis was increased in SMAD3 mutation–positive melorheostosis, with a shift of the BMDD peak toward higher mineral content; (b) cell growth was stimulated by MAP2K1 mutations but inhibited by SMAD3 mutations; and (c) effects of the activating SMAD3 mutation on osteogenic differentiation were in opposite directions depending on whether BMP2 stimulation was applied or not, but the inhibitory effect of the activating MAP2K1 mutations on osteogenic differentiation was independent of BMP2.

The study presented here demonstrated the critical role of dysregulated TGF-β/SMAD3 pathway signaling in the pathogenesis of melorheostosis with an endosteal radiographical pattern (graphical summaries in Fig. 7). The potential connections between the genes and pathways currently identified in melorheostosis are intriguing. The nature of the crosstalk between the TGF-β and BMP pathways in melorheostosis progression remains to be further elucidated. There was no evidence of alteration in LEMD3 expression by the SMAD3 or MAP2K1 mutations, and LEMD3 transcript level remained similar during osteoblast differentiation. The SMAD3 mutations (p.S264Y or p.S264F) would not be expected to impact physical interaction of SMAD3 and MAN1 (LEMD3 gene product), given that MAN1 recognizes the L3 loop of the MH2 domains of R-SMAD proteins within the helix-bundle region (Miyazono et al., 2018). However, it still remains to be determined whether the SMAD3 mutations have an influence on activities of MAN1 and PPM1A at the post-transcriptional level.

Together with our previous report (Kang et al., 2018), this study demonstrated that melorheostosis is a genetically heterogeneous disorder, with somatic mutations in different genes responsible for distinct radiographical patterns of melorheostosis through different molecular and cellular mechanisms. Our identification of MAP2K1 and SMAD3 as melorheostosis-causing genes accounts for 12 of 15 biopsied patients in our cohort. These data further underscore dysregulation of signals in the RAS/MEK/MAPK pathway and the TGF-β/SMAD pathway as the culprit in this group of skeletal disorders.

All studies were approved by NIH Intramural Institutional Review Board (protocol no. NCT02504879). Unrelated adults with a diagnosis of melorheostosis were recruited. Study subjects were not compensated for their participation. Informed consent was obtained from all patients/parents. The patients underwent open surgical biopsy of the melorheostotic lesion and of a contralateral, generally smaller, sample of unaffected bone, preferentially from the iliac crest (Jha et al., 2019; Kang et al., 2018). Two patients requested a different skeletal site for this elective procedure, and biopsy samples were obtained from the tibia and ulna (Melo-8 and Melo-17, respectively). Affected and unaffected bone samples were divided and immediately immersed either in culture media for genomic DNA extraction and cell culture or fixed in 70% ethanol before being embedded in polymethylmethacrylate for bone histology using standard procedures.

High-coverage (100×) whole exome sequencing was performed by a sequencing service provider (Otogenetics). Agilent SureSelect Human V5 (51 Mbp) kit was used for exome capture. After paired-end sequencing on Illumina HiSeq sequencers, FASTQ files were aligned to Human Genome build 37 (hg19) using Burrows-Wheeler Aligner. 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. Two somatic mutation callers (muTect v1.1.7 and Strelka v1.0.14) were used, with default parameters for whole exome sequencing 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 sequence changes or splicing changes, had <1% frequency in ExAC and 1000 Genomes Project databases, were not in a duplicated genomic region, and were called by both muTect and Strelka. The whole exome sequencing data have been deposited to the database of Genotypes and Phenotypes (dbGaP) under the accession code phs001976.v1.p1. To access these data, users may apply for access to the dbGaP data repository (https://dbgap.ncbi.nlm.nih.gov/aa/wga.cgi?login=&page=login).

In collaboration with Swift Biosciences and IDT, we designed a custom amplicon-based sequencing panel to screen genes of interest. The sensitivity of this platform can detect somatic variants with frequency as low as 0.5% VAF. Sequencing was performed to an average depth of ∼1,000× on the Illumina MiSeq instrument, with 2 × 150 bp paired-end reads. Alignment to the hg19 reference genome was done with Burrows-Wheeler Aligner mem aligner, reads were processed according to GATK best practices, and variants were called using LoFreq somatic variant caller. Subsequent filtering of variants excluded those with poor quality and common population variants.

ddPCR assays were designed for the SMAD3 p.S264Y and p.S264F mutations for the purposes of validation and variant allele quantification. The FAM- and HEX-labeled probes were specific to wild-type and mutant versions of the DNA, respectively. Droplet generation and signal quantification were performed by the Bio-Rad Q×200 ddPCR System. Calculation of VAF was performed by QuantaSoft software (Bio-Rad).

For histological evaluations, 4-µm sections were cut from the tissue block and stained with modified Goldner’s Trichrome using standard procedures (Glorieux et al., 2000). A light microscope equipped with a Zeiss AxioCam video camera was used to obtain digital images of the sections that were analyzed using NIH ImageJ software (Ver. 1.63). We evaluated osteoid thickness, osteoid surface per bone surface, osteoblast surface per bone surface, osteoclast surface per bone surface, and eroded surface per bone surface on four randomly chosen areas throughout each bone section. Bone lamellar organization was observed under polarized light. Subsequently, the residual blocks were prepared for qBEI by grinding and polishing to obtain plane parallel surfaces, then carbon coated (Roschger et al., 1998). 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 (Zeiss Supra 40; Oberkochen) equipped with a four-quadrant semiconductor backscatter electron detector. The field emission scanning electron microscope was operated with an electron energy of 20 keV. The gray levels reflecting the mineral/calcium content were calibrated by the material contrast of pure carbon and aluminum. Gray-level histograms were deduced and transformed into calcium weight percent, resulting in a corresponding BMDD curve described by five variables: CaMean, the average calcium concentration (weighted mean); CaPeak, the most frequently occurring calcium concentration (the peak position of the BMDD) in the sample; CaWidth, the width of the BMDD distribution (full width at half maximum) reflecting the heterogeneity in matrix mineralization; CaLow, the percentage of low mineralized bone area, which is mineralized below 17.68 weight percent calcium, reflecting bone areas undergoing primary mineralization; and CaHigh, the percentage of highly mineralized bone matrix, having calcium content above 25.30 weight percent calcium (Roschger et al., 2008). Statistical evaluation was performed with GraphPad Prism 6.0h. Comparison of histomorphometric and BMDD indices was based either on paired t test for comparison of affected and unaffected bone from the same patient or on unpaired t test comparing affected tissue from the present study with affected tissue from patients with somatic activating in MAP2K1 reported previously (Fratzl-Zelman et al., 2019; Kang et al., 2018). Statistical significance was considered as P < 0.05.

Primary osteoblasts were grown from freshly minced bone chips as described previously (Kang et al., 2018). Osteoblasts were plated (3,000 cells/well) in 96-well tissue culture plates and imaged at 2-h intervals for 120 h at 10× magnification using the IncuCyte ZOOM Kinetic Imaging System (Essen Bioscience). During the live-cell imaging, cells were incubated in 37°C and 8% CO₂. Percentage of confluence was calculated from the percentage of the well area occupied by cells. Doubling time was calculated using GraphPad Prism 6.0 using nonlinear regression. To analyze comparable populations of cells between groups, only wells with confluency at the first time point of image acquisition within one-half SD of the mean were included in the analysis. Statistical comparison was done by two-way ANOVA.

Constructs encoding the human SMAD3 (p.S264Y) and SMAD3 (p.S264F) mutants were generated by site-directed mutagenesis using the In-Fusion kit (Clontech). The primer pairs to introduce each mutation were as follows. SMAD3 (p.S264Y) forward: 5′-GGC​TTC​ACC​GAC​CCC​TAC​AAT​TCG​GAG​CGC​TTC​TGC​CTA​G-3′, SMAD3 (p.S264Y) reverse: 5′-GAA​GCG​CTC​CGA​ATT​GTA​GGG​GTC​GGT​GAA​GCC​ATC​CAC​AGT​C-3′; and SMAD3 (p.S264F) forward: 5′-GGC​TTC​ACC​GAC​CCC​TTC​AAT​TCG​GAG​CGC​TTC​TGC​CTA​G-3′, SMAD3 (p.S264F) reverse: 5′-GAA​GCG​CTC​CGA​ATT​GAA​GGG​GTC​GGT​GAA​GCC​ATC​CAC​AGT​C-3′. Plasmid CS2-Flag-SMAD3 (WT) was used as a template for PCR cloning to generate plasmid vectors expressing the SMAD3 mutants described in the study. The plasmid constructs were verified by Sanger sequencing. CS2-Flag-SMAD3 was a gift from Joan Massagué (Cancer Biology and Genetics Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY; Addgene plasmid #14052; Kretzschmar et al., 1997).

MC3T3-E1 cells (American Type Culture Collection) were maintained in α-MEM base media supplemented with 10% FBS (GemCell, Gemini), 100 µg/ml streptomycin, and 100 U/ml penicillin (Gibco). Cells were cultured at 37°C and 5% CO₂ and were sub-cultured upon reaching 80% confluence. For transfection, 8 × 104 cells per well were plated overnight in 12-well plates and transfected the next day with plasmid DNA constructs (1–2 µg/well) using X-tremeGENE HP DNA transfection reagent (Sigma). In some experiments, cells were serum starved 3 h with or without TGF-β receptor inhibitor SB431542 (10 µM; Sigma) and then incubated with recombinant TGF-β1 (R&D Systems) for various times as indicated.

MC3T3-E1 cells were transfected with plasmid DNA constructs expressing SMAD3 (wild-type and p.S264Y) using X-tremeGENE HP DNA Transfection Reagent. 24 h after transfection, cells were treated with or without TGF-β1 (10 ng/ml) for 30 min and fixed with 2% paraformaldehyde for 10 min at 37°C and washed five times with 1× PBS. Cells were permeabilized with 1% BSA/0.1% Triton X-100/0.2% Saponin in 1× PBS for 40 min at room temperature and washed five times with 1× PBS before incubating at 4°C overnight with primary antibodies (anti-FLAG; clone M2, #F1804; Sigma) or anti-SMAD3 (C67H9, rabbit, #9523; Cell Signaling Technology) diluted in 2% BSA/1× PBS. Cells were washed five times with 1× PBS before blocking with 5% goat serum diluted in 1× PBS for 30 min at room temperature. Cells were then stained with secondary antibodies and nuclear dye at 4°C for 1 h protecting from light. The following secondary antibodies and fluorescent dyes were used for imaging using the IncuCyte system: Alexa Fluor 568 goat anti-mouse IgG (H+L; #A11004; Invitrogen), Alexa Fluor 488 donkey anti-rabbit IgG (H+L; #A11008; Invitrogen), and Syto61 (Invitrogen).

Cultured and treated cells were lysed in radioimmunoprecipitation assay buffer and subjected to Western blot analysis, as described in the previous report (Kang et al., 2018). The following primary and secondary antibodies were used for Western blotting: phospho-Smad3 (Ser423/425; C25A9) rabbit mAb (#9520; Cell Signaling Technology), phospho-Smad3 (Ser423/425) rabbit polyclonal Ab (#600–401-919; Rockland Inc.), Smad2/3 (D7G7) XP rabbit mAb (#8685; Cell Signaling Technology), SMAD3 (C67H9) rabbit mAb (#9523; Cell Signaling Technology), β-Actin (AC-15) mouse mAb (#A5441; Sigma), GAPDH (D16H11) XP rabbit mAb (#5174; Cell Signaling Technology), histone H3 (D1H2) XP rabbit mAb (#4499; Cell Signaling Technology), and FLAG (M2) mouse mAb (#F1804; Sigma). Secondary antibodies (LI-COR Biosciences) were IRDye 680RD anti-rabbit IgG (#926-68071) or IRDye 800CW anti-mouse IgG (#926-32210).

Total RNA was extracted from osteoblasts using the RNeasy Mini kit (Qiagen) following the manufacturer’s instruction. Genomic DNA was removed from RNA preps by on-column DNase digestion. Synthesis of cDNA was performed using 500 ng of total RNA and iScript reverse transcription Supermix (Bio-Rad), following the manufacturer’s instructions. Comparative real-time qPCR was performed in TaqMan Universal PCR Master Mix in triplicate using Applied Biosystems Prism 7500 Fast Sequence Detection System, according to the manufacturer’s protocol. TaqMan probes and primers were purchased from Applied Biosystems: Human COL1A1 (Hs00164004_m1), COL1A2 (Hs01028956_m1), CDKN2B (Hs00793225_m1), SERPINE1 (Hs00167155_m1), ALPL (Hs01029144_m1), SP7 (Hs01866874_s1), CDH11 (Hs00901479_m1), MMP13 (Hs00942584_m1), and SPP1 (Hs00959010_m1). Mouse Ctgf (Mm01192933_g1), Serpine1 (Mm00435858_m1), and Vegfa (Mm00437306_m1). Human GAPDH (Hs02786624_g1), human TBP (Hs00427620_m1), mouse Gapdh (Mm99999915_g1), and mouse Tbp (Mm00446971_m1) were used as endogenous controls for normalization of real-time qPCR. Relative expression was calculated using the comparative ΔΔCt method.

Primary patient osteoblasts were plated at passage 2 with 50,000 cells/well in 12-well tissue culture plates and cultured to confluency in α-MEM supplemented with 10% FBS and antibiotics at 37°C and 8% CO₂. After confluence, cells were further cultured in osteogenic media (50 µg/ml L-ascorbic acid, 2.5 mM β-glycerophosphate, and 10 nM dexamethasone, all from Sigma), with or without recombinant BMP2 (100 ng/ml; R&D Systems), refreshing osteogenic media every 3 d. To visualize ECM mineralization, cells were washed once with 1× PBS and fixed with 4% paraformaldehyde for 10 min at room temperature. Fixed cells were washed three times with 1× PBS and stained with 2% Alizarin Red S solution (pH 4.2) for 10 min at room temperature. Excess Alizarin Red S stain was removed by washing five times with distilled water. Cells were air-dried in the dark before imaging.

Total RNA samples (∼1–4 µg) were purified with Poly-A extraction, and then purified mRNAs were used to construct RNA-Seq libraries with specific bar codes using the Illumina TruSeq Stranded mRNA Library Prep Kit. All the RNA-Seq libraries were pooled together and sequenced using Illumina HiSeq 2500 to generate ∼40 million 2 × 100 paired-end reads for each sample. Demultiplexed 100-bp paired-end fastq reads were aligned to GRCh38 human reference using GENCODE release 27 gene definitions using RNA STAR v2.5.4a (Dobin et al., 2013). Following alignment, reads were quantitated using the subread package featureCounts v1.5.2 (Liao et al., 2014) against GENCODE release 27 gene definitions, counting multi-mapping and multi-overlapping reads at the level of gene names. All downstream analyses were performed in R v3.5.1, Ingenuity Pathway Analysis (Qiagen Bioinformatics), and Metascape. For each of these conditions, we loaded the counts data into DESeq2 v1.20 (Love et al., 2014) and performed differential expression analysis between affected and unaffected samples at each time point. Since the samples are technical replicates of cells from one melorheostosis patient, this resulted in lower reported P values than biological replicates. To compensate for this, we chose a more stringent threshold of FDR < 1e-20 to select DEGs. Other thresholds return largely the same results.

To identify patterns of expression throughout the data, we used the degPatterns function from the DEGreport v1.16 Bioconductor package (Steinbaugh et al., 2018). Briefly, this analysis takes normalized counts from genes of interest, generates a correlation matrix of genes across condition and time, clusters that matrix, and cuts the resulting tree into modules of coexpressed genes. We used counts normalized with the DESeq2 variance stabilizing transformation and provided counts for the union of genes differentially expressed between affected and unaffected at any time point. Resulting clusters were merged together if the pairwise correlation coefficient was >0.5, and clusters with <15 genes after merging were dropped. The plots show normalized counts for the genes in each cluster. MA plots show the log of mean normalized counts for each gene on the x axis and the log₂ fold change on the y axis. These plots reflect the default behavior of DESeq2 to perform fold-change shrinkage, shifting the reported log fold change of genes with low information toward zero. Heatmaps show normalized counts after applying the variance-stabilizing transform, with genes subsetted from the full normalized counts depending on the gene list of interest. Plots were generated with the heatmaply v0.15.2 package, scaling by row, disabling column clustering, and otherwise using defaults. Gene set enrichment analysis (GSEA) analysis was performed using the clusterProfiler v3.8.1 package’s GSEA function (Yu et al., 2012) after providing a custom annotation containing only the set of genes of interest. The RNA-Seq data have been deposited to dbGaP under the accession code phs001976.v1.p1. To access these data, users may apply for access to the dbGaP data repository (https://dbgap.ncbi.nlm.nih.gov/aa/wga.cgi?login=&page=login).

For bone histology, statistical evaluation was performed with GraphPad Prism 6.0h. Comparison of histomorphometric and BMDD indices was based either on paired t test for comparison of affected and unaffected bone from the same patient or on unpaired t test comparing affected tissue from the present study with affected tissue from MAP2K1 mutation–positive patients. Statistical significance was considered as P value <0.05. For cell proliferation assay and SMAD3 nuclear localization, statistical comparison was done by two-way ANOVA. For real-time qPCR data analysis, differences between groups were evaluated with the Student’s t test using a two-tailed distribution. Results were considered significant when P value was <0.05. Data are expressed as mean ± SD.

Fig. S1 shows radiographs of three patients with SMAD3 mutation–positive melorheostosis demonstrating the endosteal appearance. Fig. S2 shows the results of a genetic study identifying SMAD3 somatic mutations in melorheostosis. Fig. S2 also illustrates human SMAD3 missense mutations found in melorheostosis, CRC, and LDS. Fig. S3 shows effects of SMAD3 mutations on TGF-β signaling. Fig. S4 shows the results of transcriptome profiling by RNA-Seq. Fig. S5 shows BMP2-independent inhibition of ECM mineralization by mosaic-activating MAP2K1 mutation. Table S1 summarizes the clinical findings on SMAD3 mutation–positive melorheostosis patients. Table S2 summarizes the bone histomorphometry and qBEI analyses in melorheostosis bone. Table S3 shows β-catenin (CTNNB1) up-regulation and changes of downstream target genes in affected cells of Melo-11. Table S4 shows up-regulation of interferon-inducible signature genes (21 genes) in affected cells of melorheostosis (Melo-11, SMAD3, p.S264Y). Data S1 lists the genes in each cluster generated by degPatterns analysis of RNA-Seq data of BMP2(−) condition. Data S2 lists the genes in each cluster generated by degPatterns analysis of RNA-Seq data of BMP2(+) condition. Data S3 lists enriched pathway terms depicted as an interaction network for each cluster generated through the Metascape analysis (BMP2[−] condition). Data S4 lists enriched pathway terms depicted as an interaction network for each cluster generated through the Metascape analysis (BMP2[+] condition).

The authors thank the Center for Cancer Research Genomics Core facility (National Cancer Institute) for conducting ddPCR assays and Daniela Gabriel, Petra Keplinger, Sonja Lueger, and Phaedra Messmer (Ludwig Boltzmann Institute of Osteology at the Hanusch Hospital of Wiener Gebietskrankenkasse and Allgemeine Unfallversicherungsanstalt Trauma Center Meidling) for careful sample preparations and qBEI measurements. We thank Tianwei Li and James Iben at the Molecular Genomics Core (National Institute of Child Health and Human Development) for RNA-Seq and primary data processing, and Nicole Swan (National Institute of Child Health and Human Development) for assistance with graphics. We are especially grateful to Kathleen Harper, the Melorheostosis Association, and the patients for their support of this study.

This work was supported by National Institutes of Health intramural funds to J.C. Marini (National Institute of Child Health and Human Development, 1ZIAHD008973-02), R.M. Siegel (National Institute of Arthritis and Musculoskeletal and Skin Diseases, AR041173-11, AR041207-03, and AR041180-11), and T. Bhattacharyya (National Institute of Arthritis and Musculoskeletal and Skin Diseases, 1ZIDAR041180-09), and the Allgemeine Unfallversicherungsanstalt (research funds of the Austrian Workers’ Compensation Board) and Wiener Gebietskrankenkasse (Vienna Regional Health Insurance Fund) to N. Fratzl-Zelman. A. Ivovic was supported by the National Institutes of Health-Oxford-Cambridge PhD Scholars Program.

Author contributions: T. Bhattacharyya, S. Jha, E. Lange, E.W. Cowen, and J. Katz wrote the clinical protocol and provided patient care. H. Kang, A. Ivovic, S. Jha, Z. Deng, N. Fratzl-Zelman, A. Mitra, W.A. Cabral, E.P. Hanson, P. Roschger, K. Klaushofer, R.K. Dale, and T. Bhattacharyya designed and conducted experiments and analyzed data. The manuscript was drafted by H. Kang and edited by J.C. Marini with significant input from A. Ivovic, S. Jha, Z. Deng, N. Fratzl-Zelman, R.M. Siegel, and T. Bhattacharyya.