Precocious chondrocyte differentiation disrupts skeletal growth in Kabuki syndrome mice

Mendelian disorders of the epigenetic machinery disrupt the fundamental processes of neurological development and growth (1). This rapidly growing group of inherited conditions resulting from germ-line mutations in components of the epigenetic machinery is expected to have broad epigenomic consequences. Despite growth abnormalities being the second most common disease manifestation, molecular underpinnings have not been examined in detail but could provide insight into disease mechanisms that may be broadly applicable to other more common growth disorders, like idiopathic short stature.

Kabuki syndrome (KS; MIM 147920), an autosomal dominant Mendelian disorder of the epigenetic machinery (1), results from heterozygous, typically de novo inactivating mutations in KMT2D (2), which encodes an epigenetic writer that normally catalyzes histone methylation on H3 at lysine 4 (H3K4me). Individuals with Kabuki syndrome 1 (KS1) exhibit the cardinal phenotypic features of postnatal growth retardation, intellectual disability, and craniofacial abnormalities (3, 4). The distinct craniofacial features that are characteristic of KS1 include flattening of the facial profile, elongated palpebral fissures with eversion of the lower eyelids, highly arched and interrupted eyebrows, short columella with a depressed nasal tip, prominent ears, and palate abnormalities; these often provide the best clue for clinical diagnosis in affected individuals (5–7). Schott et al. systematically evaluated growth patterns in individuals with molecularly confirmed KS1, revealing postnatal growth retardation in the vast majority (8) and confirming previous reports from clinically diagnosed patients (9, 10). Most of the molecularly confirmed individuals with KS1 did not meet criteria for growth hormone deficiency (11). Despite this, treatment with recombinant growth hormone therapy improved linear growth in some but not all individuals with KS1 (12). Such observations demonstrate the lack of knowledge about the underlying mechanism of growth retardation in KS1 and suggest that it is more complex than isolated growth hormone deficiency.

Using morpholinos to knock down Kmt2d in zebrafish, Van Laarhoven et al. postulated a defect in neural crest–derived cell function in the development of facial flattening in fish deficient in Kmt2d because multiple cartilaginous structures were underdeveloped (13). In our previous work, we characterized a mouse model of KS1, which exhibits many features seen in patients with the disorder (14). Kmt2d^+/βGeo mice have disrupted H3K4 trimethylation (H3K4me3) in the dentate gyrus of the hippocampus and associated neurogenesis defects and memory deficits. All 3 features could be reversed by postnatal administration of agents that favor chromatin opening such as HDAC inhibitors or the HDAC inhibitor–like ketone body, β-hydroxybutyrate (14, 15). Kmt2d^+/βGeo mice also weigh less than Kmt2d^+/+ littermates and exhibit a flattened facial profile like individuals with KS1 (14). The cellular and molecular basis for this phenotype has not been examined.

Here, using high-resolution micro-CT, we elucidate a robust and quantitative skeletal growth retardation phenotype impacting the long bones and cranial base of the skull in Kmt2d^+/βGeo mice. Histological data from growth plates within both sites suggest a unifying mechanism of disrupted endochondral ossification. Our in vitro studies of stable chondrocyte cell lines harboring loss-of-function Kmt2d mutations provide further evidence that precocious differentiation of chondrocytes — and increased proliferation, to a lesser extent — play a key role in KS1 pathogenesis. Targeted and genome-wide transcriptome profiling supports precocious differentiation of chondrocytes and, together with targeted chromatin immunoprecipitation studies, suggests the basis for Sox9 dysregulation. This involves loss of KMT2D-mediated H3K4me3 at an unsuspected target, Shox2, causing release of Sox9 inhibition and, thereby, precocious chondrocyte differentiation. Our findings provide mechanistic insight into the pathogenesis of growth retardation in KS1 and related disorders and suggest potentially novel therapeutic targets.

Kmt2d^+/βGeo mice exhibit a specific skeletal growth retardation phenotype. Kmt2d^+/βGeo mice are smaller than Kmt2d^+/+ littermates, both grossly (14) and on lateral radiographs (Figure 1A). Quantification of body weight (Figure 1B) and length (Figure 1C) with repeated observations using mice at multiple ages (Supplemental Figure 1A; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.129380DS1) confirmed growth retardation in Kmt2d^+/βGeo mice. Similar to children with KS (8), decreases in weight and length were proportional, as comparison of weights normalized to lengths in Kmt2d^+/βGeo and Kmt2d^+/+ mice showed similar trends but were mostly not significant (data not shown). Because detailed quantification on radiographs revealed a decrease in upper jaw length in Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates (data not shown), we first performed high-resolution micro-CT analysis of the skulls of Kmt2d^+/βGeo mice and Kmt2d^+/+ littermates to examine the craniofacial phenotype. Reconstructed micro-CT images unequivocally confirmed a striking flattening of the facial profile in Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates (Figure 1D), which resembles the facial flattening seen in individuals with KS (16).

Figure 1

Kmt2d^+/βGeo mice exhibit generalized growth retardation and a specific craniofacial phenotype reminiscent of individuals with KS1. (A) Representative radiographs of Kmt2d^+/βGeo mice and Kmt2d^+/+ littermates illustrating growth retardation and flattening of the facial profile. Quantification of (B) body weight and (C) body length in 6-week-old Kmt2d^+/βGeo male (n = 6) and female (n = 3) mice and Kmt2d^+/+ male (n = 3) and female (n = 7) littermates. Data represent mean ± SD, and similar results were obtained with multiple cohorts of mice. One-sided unpaired Student’s t test *P < 0.05; **P < 0.01; ***P < 0.001. (D) Representative reconstructions of high-resolution craniofacial micro-CTs in the left lateral view from Kmt2d^+/βGeo mice and Kmt2d^+/+ littermates illustrating craniofacial phenotype and (E) showing, in green, 4 pairs of bilateral landmarks (top) and 10 midline landmarks (bottom) used for morphometric analysis. Scale bar: 3 mm. (F) PCA of shape revealing separation of 2 distinct groups along PC1 and PC3, with Kmt2d^+/+ mice (n = 21) toward the lower end and Kmt2d^+/βGeo mice (n = 13) toward the upper end of PC1. Two-way Procrustes ANOVA confirmed a statistically significant effect of genotype on overall cranial shape (Pillai’s trace = 0.98; P = 0.0086). (G) Overlay of wire frames in left lateral view illustrating relative differences in shape of Kmt2d^+/+ mice (blue) and Kmt2d^+/βGeo mice (red). Black vectors show displacement of landmarks associated with the range of shape variation observed on PC1 and indicate ventral bowing, dorsal expansion, and brachycephaly in KS1 (thick black arrows).

To further characterize the craniofacial phenotype, we performed morphometric analyses of 3-dimensional (3-D) reconstructions of the micro-CT data. Principal component analysis (PCA) using 18 landmarks (Figure 1E) revealed that Kmt2d^+/βGeo and Kmt2d^+/+ mice fall into 2 distinct groups, with Kmt2d^+/βGeo mice toward the upper end of principal component 1 (PC1) and approximately half of the total variance across the combined sample being attributable to ventral bowing and brachycephaly (PC1) in Kmt2d^+/βGeo mice (Figure 1, F and G, and Supplemental Figure 2). This analysis also indicates that the Kmt2d^+/βGeo mice exhibit a dorsal expansion of the cranial vault compared with Kmt2d^+/+ littermates (Figure 1G and Supplemental Figure 2).

Because long bones also appeared to be shortened on lateral radiographs (Figure 1A) and generalized postnatal growth retardation is a key aspect of the KS skeletal phenotype in patients, we performed analyses of femur and tibial length and high-resolution micro-CT evaluation of bone architecture. Both femurs (Figure 2, A and B) and tibias (Figure 2, C and D) were significantly shorter in Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates at 6 and 18 weeks of age (Figure 2, A–D, and Supplemental Figure 1, B and C). High-resolution micro-CT revealed differences in cortical and trabecular bone structure in Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates. In the femur, cortical cross-sectional tissue area did not differ between groups (Figure 2, E and F), and mineralized bone area was reduced in Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates (Figure 2, E and G), which led to a significant decrease in the percent bone area per tissue area in Kmt2d^+/βGeo males (Figure 2, E and H). Trabecular bone volume examined in the distal femur was reduced in Kmt2d^+/βGeo male mice (Figure 2, I and J) secondary to decreases in trabecular number (Figure 2, I and K) and trabecular thickness (Figure 2, I and L); however, in female Kmt2d^+/βGeo mice, only trabecular thickness was significantly reduced (Figure 2, I–L). Upon adjusting the region of interest for analysis of trabecular bone parameters to account for reduced femur length in Kmt2d^+/βGeo male mice, similar results were obtained (data not shown). Therefore, Kmt2d^+/βGeo mice exhibit a distinct defect in skeletal growth that is evident at multiple sites and resembles that observed clinically in KS1.

Figure 2

High-resolution micro-CT analysis of long bones in KS1 reveals shortening, thinning, and altered trabecular bone formation in Kmt2d^+/βGeo mice. (A and B) Femurs and (C and D) tibias are shorter in Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates. (E and F) Overall cross-sectional area does not differ between Kmt2d^+/βGeo and Kmt2d^+/+ femurs. (E and G) Cross-sectional area of mineralized bone is reduced in Kmt2d^+/βGeo femurs compared with Kmt2d^+/+ femurs. (E and H) The percent of cross-sectional area made up of mineralized bone is reduced in male Kmt2d^+/βGeo femurs compared with Kmt2d^+/+ femurs, whereas the difference is not significant in females. (I) Kmt2d^+/βGeo femurs appear to have decreased trabecular bone near the growth plate compared with Kmt2d^+/+ femurs. (I and J) Specifically, percent of tissue volume made up of bone is decreased in male Kmt2d^+/βGeo femurs, and (I and K) trabecular number is decreased in male Kmt2d^+/βGeo femurs. (I and L) Trabecular thickness is decreased in male and female Kmt2d^+/βGeo femurs. For Kmt2d^+/+ femurs, n = 18 (8 male and 10 female); for Kmt2d^+/+ tibias n = 15 (7 male and 8 female); for Kmt2d^+/βGeo femurs, n = 13 (7 male and 6 female, except for femur length, where only 6 male mutants could be measured); for Kmt2d^+/βGeo tibias, n = 11 (5 male, 6 female). Data represent mean ± SD. One-tailed (B and D) or 2-tailed (F–H and J–L) unpaired Student’s t tests were used. *P < 0.05; **P < 0.01; ***P < 0.001.

In vivo and in vitro studies suggest disrupted endochondral ossification in KS1. Long bones, as well as a few bones in the skull, increase in length by means of endochondral ossification at growth plates, which involves proliferation and hypertrophy of chondrocytes and production of the cartilaginous matrix that is ultimately replaced by bone laid down by osteoblasts. The pattern of growth abnormalities in Kmt2d^+/βGeo mice suggests alterations in growth plate dynamics as a potential mechanism for the defects observed in long bones and the cranial base. Therefore, we stained longitudinal sections from Kmt2d^+/βGeo and Kmt2d^+/+ proximal tibia growth plates with H&E. While overall morphology was intact, we found that the growth plates were expanded in Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates (Figure 3A). Growth plate height (Figure 3, A and B), proliferative zone height (Figure 3, A and C), and hypertrophic zone height (Figure 3, A and D) were increased in both male and female Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates at comparable growth plate regions (Supplemental Figure 3A). Counting of the number of chondrocytes per column in the proliferative zone (Figure 3, A and E) and in the hypertrophic zone (Figure 3, A and F) revealed increased cell numbers in Kmt2d^+/βGeo mice compared with Kmt2d^+/+ mice at both sites. Whereas no significant difference in cell size was evident in males, female Kmt2d^+/βGeo mice showed increased hypertrophic cell size compared with Kmt2d^+/+ mice (Supplemental Figure 3B). There was no difference in cell death between genotypes (Supplemental Figure 3, C and D).

Figure 3

Growth plates from long bones and within the cranial base are expanded in Kmt2d^+/βGeo mice. (A and B) Proximal tibia growth plates and their (A and C) proliferative and (A and D) hypertrophic zones (PZ and HZ, respectively) are expanded in Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates. The mechanism involves increased cell numbers per column in both the (E) PZ and (F) HZ. (G and H) Cranial base intrasphenoidal synchondroses (ISS) from Kmt2d^+/βGeo mice are expanded compared with Kmt2d^+/+ littermates. GP, growth plate. For Kmt2d^+/+ proximal tibia growth plates, n = 16; 7 male and 9 female. For Kmt2d^+/+ ISS growth plates, n = 20; 11 male, 9 female. For Kmt2d^+/βGeo proximal tibia growth plates, n = 11; 5 male and 6 female. For Kmt2d^+/βGeo ISS growth plates, n = 13; 7 male and 6 female. Data represent mean ± SD. One-sided unpaired Student’s t test; *P < 0.05; **P < 0.01; ***P < 0.001. Scale bar: 100 μM.

We next examined the cranial base, a key area of the skull that grows by endochondral ossification, and another site that also appears to be shortened in individuals and mice with KS1 (Figure 1). Examination of intrasphenoidal synchondroses, which are partly responsible for growth of the cranial base in the anterior-posterior direction (17), revealed a similar expansion in Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates (Figure 3, G and H). While counterintuitive that growth plate height could be expanded while the cranial base and long bones themselves were shortened in length, these findings suggest that abnormal endochondral ossification of long bones and of the cranial base are a unifying mechanism for the skeletal phenotype of KS1 mice.

To understand the cellular basis for this phenotype, we employed a more manipulatable system of chondrocyte development, the well-established ATDC5 cell system; upon induction, these mesenchymal progenitors progress through the stages of chondrocyte development, including proliferation and hypertrophy (18). Using CRISPR-Cas9 genome editing technology, we created stable cell lines with compound heterozygous mutations of varying severity — Kmt2d^ΔR5551/–, which has a deletion of a single amino acid on 1 allele and deletion of the SET domain on the other allele, and Kmt2d^–/–, which has biallelic deletions of the SET domain — and compared them with vector-only parental ATDC5 cells (Kmt2d^+/+; Figure 4, Supplemental Figure 4, and Supplemental Figure 5) before and at distinct independent time points after induction of chondrocyte differentiation. Assessment of cellular proliferation by the MTT assay (Supplemental Figure 5A) and direct cell counting (Supplemental Figure 5B) suggested that loss of Kmt2d function may lead to increased chondrocyte replication in accordance with the severity of the mutation, supporting our previous observation of proliferative zone expansion in proximal tibia growth plates from Kmt2d^+/βGeo mice (Figure 3, A, C, and E); however, differences between genotypes were not significant. Interestingly, cultures of Kmt2d^–/– cells exhibited more dramatic increases in matrix deposition, as indicated by Alcian blue staining, when compared with Kmt2d^+/+ controls, with Kmt2d^ΔR5551/– cells having an intermediate phenotype (Figure 4, A and B). By 14 days after induction of differentiation, we observed a significant increase in Alcian blue staining of Kmt2d^–/– chondrocytes compared with Kmt2d^ΔR5551/– and Kmt2d^+/+ chondrocytes (Figure 4, A and B). This increase was progressive and exceedingly apparent at day 21. These data suggest precocious differentiation of Kmt2d^–/– cells and fit with the expansion of hypertrophic zones observed in proximal tibia growth plates from Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates (Figure 3, A, D, and F). Together, these findings show that all cell lines (Kmt2d^+/+, Kmt2d^ΔR5551/–, and Kmt2d^–/–) proceed through the phases of chondrocyte development (i.e., proliferation [Supplemental Figure 5, A and B] and hypertrophy [Figure 4, A and B]); however, Kmt2d^–/– chondrocytes trend toward increased proliferation and exhibit striking precocious differentiation compared with Kmt2d^+/+ chondrocytes, with Kmt2d^ΔR5551/– cells having an intermediate phenotype (Figure 4, A and B, and Supplemental Figure 5, A and B). These observations support our in vivo findings of growth plate expansion in both tibias and cranial base intrasphenoidal synchondroses of Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates (Figure 3), and together, these data support precocious differentiation (and increased proliferation, to a lesser extent) as a unifying feature of KMT2D-deficient chondrocytes both in vitro and in vivo.

Figure 4

Precocious differentiation of Kmt2d^–/– and Kmt2d^ΔR5551/– chondrocytes. Stable ATDC5 cell lines with Kmt2d mutations of varying severity were created using CRISPR/Cas9 genome editing technology. Chondrocyte differentiation was induced at day 0. Alcian blue staining was used to (A) visualize and (B) quantify chondrocyte differentiation over time. Scale bar: 1 mm. Black asterisks represent differences between Kmt2d^+/+ cells and Kmt2d^–/– or Kmt2d^ΔR5551/– cells; gray asterisks represent differences between Kmt2d^–/– and Kmt2d^ΔR5551/– cells. qPCR with primers specific for (C) Col2a1, (D) Col10a1, and (E) Sox9 was performed on RNA isolated from undifferentiated (day 0) and differentiated Kmt2d^+/+, Kmt2d^ΔR5551/–, and Kmt2d^–/– stable chondrocyte cell lines at 4, 7, 14, and 21 days after induction of differentiation. Fold change was calculated relative to undifferentiated Kmt2d^+/+ cells (day 0). (B–E) Data represent mean ± SEM; the experiments were performed 2 times with 3 independent stable cell lines per time point (n = 6) for Kmt2d^+/+ and Kmt2d^ΔR5551/– and with 2 independent stable cell lines per time point (n = 4) for Kmt2d^–/–. Mixed model ANOVA with Tukey’s adjustment method within each time point was used; *P < 0.05; **P < 0.01. (F) Preliminary model for precocious chondrocyte differentiation in KS1.

Precocious chondrocyte expression patterns in a cellular model of KS1. Using quantitative PCR (qPCR), we measured expression of chondrocyte differentiation markers Col2a1 and Col10a1 in undifferentiated mesenchymal cells and in differentiated chondrocytes with independent experiments performed at each time point (Figure 4, C and D). Col2a1 is normally expressed early in the chondrocyte developmental program, initially within the proliferative phase, and Col10a1 is expressed later in more advanced stages of chondrocyte development beginning in the early hypertrophic phase (19). Upon differentiation induction, we observed a striking dose-dependent increase in Col2a1 expression in Kmt2d^ΔR5551/– and Kmt2d^–/– chondrocytes compared with Kmt2d^+/+ chondrocytes; these results were significant between Kmt2d^–/– and Kmt2d^+/+ chondrocytes at 4–14 days after induction of differentiation (Figure 4C). While a gradual increase in Col2a1 expression is expected and this trend was observed over time (days 4–14) in Kmt2d^+/+ chondrocytes, in Kmt2d^ΔR5551/– and Kmt2d^–/– chondrocytes, Col2a1 expression increased precociously, attaining higher maximal levels earlier (Figure 4C). Similarly, upon induction of chondrocyte differentiation in Kmt2d^+/+ cells, expression of the normally later Col10a1 chondrocyte differentiation marker initially dropped to almost undetectable levels (day 4) but then increased to maximal levels by days 14–21 (Figure 4D), consistent with Col10a1 being a late marker of hypertrophic chondrocyte differentiation (20). Over time, both Kmt2d^ΔR5551/– and Kmt2d^–/– chondrocyte cell lines followed a pattern similar to Kmt2d^+/+ (Figure 4D); however, Col10a1 transcript levels were significantly higher in Kmt2d^–/– cells compared with Kmt2d^+/+ cells by day 4, after differentiation induction (Figure 4D), suggesting precocious expression of this normally late chondrogenic marker. Kmt2d^ΔR5551/– cells showed intermediate levels of Col10a1 expression overall, more closely resembling Kmt2d^–/– cells until day 4 and Kmt2d^+/+ cells thereafter (Figure 4D).

Overexpression of Sox9 mediates the precocious differentiation phenotype. We hypothesized that Sox9, which encodes a key transcription factor that promotes chondrocyte differentiation and is a known direct activator of Col2a1 and Col10a1 (19), could be a potential candidate gene to mediate our cellular phenotype. SOX9 facilitates chondrocyte differentiation through proliferation and hypertrophy and prevents further progression toward bone-forming osteoblasts (19). Upon induction of chondrocyte differentiation, Sox9 expression in Kmt2d^+/+ cells gradually increased by 3- to 4-fold over 21 days, and by 4 days after differentiation induction, we observed a clear trend toward increased Sox9 expression in Kmt2d^ΔR5551/– and Kmt2d^–/– cells dependent upon the severity of the KMT2D mutation (Figure 4E). Like Col2a1 and Col10a1, which are induced by SOX9, Sox9 is expressed at higher levels in Kmt2d^–/– compared with Kmt2d^+/+ chondrocytes (and in Kmt2d^–/– compared with Kmt2d^ΔR5551/– chondrocytes, to a lesser extent) up to 14 days after induction of differentiation, and these results are significant at day 14 (Figure 4E). For Sox9, this culminates in a striking 4-fold increase in expression in Kmt2d^–/– chondrocytes compared with Kmt2d^+/+ chondrocytes at day 14 (Figure 4E). The observation that genetic ablation of KMT2D, a histone methyltransferase writer of an activating mark, increases Sox9 suggests a model by which loss of KMT2D may prevent expression of an inhibitor of Sox9, thus allowing subsequent chondrocyte differentiation to proceed precociously (Figure 4F).

Loss of Shox2 expression and release of Sox9 inhibition in a KS1 chondrocyte model. We performed genome-wide transcriptome profiling by RNA sequencing (RNA-seq) using RNA from Kmt2d^–/– and Kmt2d^+/+ chondrocytes 7 days after induction of differentiation to identify KMT2D targets that can inhibit Sox9 (Figure 5A, Supplemental Figure 6A, and Supplemental Table 1). Day 7 appears to be a key time point for gene expression changes that lead to precocious differentiation because, after day 7, Alcian blue staining (Figure 4, A and B) and expression of Sox9 (Figure 4E) increase significantly in Kmt2d^–/– chondrocytes compared with Kmt2d^+/+ chondrocytes. For comparison, we also performed RNA-seq on undifferentiated Kmt2d^–/– and Kmt2d^+/+ mesenchymal cells (Figure 5B, Supplemental Figure 6B, and Supplemental Table 2).

Figure 5

Genome-wide transcriptome profiling identifies Shox2 as a target of KMT2D in vitro and in vivo. Volcano plots showing genome-wide differential gene expression in Kmt2d^–/– and Kmt2d^+/+ (A) chondrocytes differentiated for 7 days and (B) undifferentiated cells. (C) PCA reveals tight clustering within and distinct separation between genotypes and differentiation states, separating the cells into 4 distinct groups. Kmt2d^+/+ cells cluster toward the upper end of PC2 and Kmt2d^–/– cells cluster toward the lower end of PC2. Undifferentiated cells cluster toward the lower end of PC1 and chondrocytes cluster toward the upper end of PC1. Each color represents a biological replicate; each biological replicate is a distinct clonal cell line. Three (Kmt2d^+/+) or 2 (Kmt2d^–/–) biological replicates were used for each differentiation state (chondrocytes vs. undifferentiated cells). Each point represents a technical replicate; 2 technical replicates were performed for each cell line. Sox9 and Col2a1 are upregulated and Shox2 is downregulated in (A) Kmt2d^–/– chondrocytes and (B) Kmt2d^–/– undifferentiated cells. (D) Venn diagrams summarizing differential gene expression from RNA-seq analyses in Kmt2d^+/+ and Kmt2d^–/– chondrocytes and undifferentiated cells. (E) Validation of downregulation of Shox2 in 3 independent experiments in vitro; mean ± SEM; n = 9 (Kmt2d^+/+) or n = 6 (Kmt2d^–/–) per group (undifferentiated cells or chondrocytes); 1-sided unpaired Student’s t test; ***P < 0.001; ****P < 0.0001. (F) Immunofluorescence and (G) quantification of SHOX2 protein levels in vivo in Kmt2d^+/βGeo and Kmt2d^+/+ proximal tibia growth plates (outlined in yellow dashed lines). SHOX2 intensity was normalized to number of DAPI-stained cells within the growth plate. Data represent mean ± SD; n = 5 male Kmt2d^+/+, n = 5 male Kmt2d^+/βGeo, n = 3 female Kmt2d^+/+, n = 3 female Kmt2d^+/βGeo; 1-sided unpaired Student’s t test; **P < 0.01. Scale bar: 400 μM. Chondro, chondrocytes differentiated for 7 days; Undiff, undifferentiated cells.

Multiple observations emerged from our transcriptome profiling. First, PCA revealed highly correlative data with tight clustering within and distinct separation between genotypes and differentiation states, separating the cells into 4 distinct groups (Kmt2d^+/+ undifferentiated cells, Kmt2d^+/+ chondrocytes, Kmt2d^–/– undifferentiated cells, and Kmt2d^–/– chondrocytes) with differentiation state (chondrocytes vs. undifferentiated cells) explaining a greater percentage of the variance than genotype (Kmt2d^–/– vs. Kmt2d^+/+; Figure 5C). Second, differential expression analysis revealed multiple expected gene expression changes, including greater than 50-fold upregulation of Col2a1 in Kmt2d^–/– cells compared with Kmt2d^+/+ cells, both in chondrocytes (day 7; Figure 5A, Supplemental Figure 6A, and Supplemental Table 1) and in the undifferentiated state (day 0; Figure 5B, Supplemental Figure 6B, and Supplemental Table 2), similar to our previous observations (Figure 4C). In addition, we observed consistent 3- to 4-fold upregulation of the master transcriptional regulator of chondrogenesis, Sox9, in Kmt2d^–/– cells compared with Kmt2d^+/+ cells in both differentiation states (Figure 5, A and B, and Supplemental Tables 1–3), again supporting our previous findings (Figure 4E). Additional well-established markers of chondrocyte differentiation (21) were upregulated, as well, in Kmt2d^–/– cells compared with Kmt2d^+/+ cells, including Col11a1 (21) (over 10-fold in both chondrocytes and undifferentiated cells; Supplemental Tables 1–4), Col9a1 (21) (over 10-fold in chondrocytes; Supplemental Table 1), and Col10a1 (4-fold in undifferentiated cells; Supplemental Table 2; compare with Figure 4D), further supporting a mechanism involving precocious chondrocyte differentiation. Lastly, we observed more differential gene expression within a single genotype, Kmt2d^+/+ (7,201 genes; Supplemental Figure 6, C and E, and Supplemental Table 5) or Kmt2d^–/– (4,764 genes; Supplemental Figure 6, D and F, and Supplemental Table 6), over the course of differentiation from mesenchymal cells to chondrocytes compared with the differential expression observed between genotypes within a particular differentiated state (i.e., chondrocytes Kmt2d^–/– vs. Kmt2d^+/+ [932 genes; Figure 5D and Supplemental Table 1] or undifferentiated mesenchymal cells Kmt2d^–/– vs. Kmt2d^+/+ [1,009 genes; Figure 5D and Supplemental Table 2]). This is expected and suggests that the number of gene expression changes required for a cell to undergo a differentiation program is greater than the number that result from loss of function of a single gene, even if that gene is a component of the epigenetic machinery with many downstream direct targets and indirect consequences.

Comparing Kmt2d^–/– chondrocytes with Kmt2d^+/+ chondrocytes, we found 932 differentially expressed genes (Figure 5D) with 486 upregulated and 446 downregulated (Figure 5A, Supplemental Figure 6A, and Supplemental Table 1). There were 1,009 differentially expressed genes when Kmt2d^–/– and Kmt2d^+/+ undifferentiated cells were compared (Figure 5D), with 474 being upregulated and 535 being downregulated (Figure 5B, Supplemental Figure 6B, and Supplemental Table 2). Of those differentially expressed genes, 382 are common across differentiation state (i.e., found in both chondrocytes and undifferentiated cells; Figure 5D and Supplemental Table 3), and this represents a highly significant overrepresentation (Fisher’s test, P < 2.2 × 10^–16, OR = 14.2). Of the 932 differentially expressed genes in chondrocytes, 265 were differentially expressed at least 4-fold (Figure 5, A and D, and Supplemental Figure 6A), and of the 1,009 differentially expressed genes in undifferentiated cells, 180 were differentially expressed least 4-fold (Figure 5, B and D, and Supplemental Figure 6B). Seventy-eight are common across chondrocytes and undifferentiated cells (Figure 5D and Supplemental Table 4), and this overlap is again highly significant (Fisher’s test, P < 2.2 × 10^–16, OR = 78.1). Of those 78 differentially expressed transcripts (≥4-fold) common to both chondrocytes and undifferentiated cells, 74 had established gene names (Figure 5D). Of those 74 genes, 30 were upregulated and 44 were downregulated in Kmt2d^–/– cells (Figure 5D and Supplemental Table 4). We hypothesized that direct KMT2D targets would have more dramatic gene expression abnormalities and, thus, focused on the subset that were downregulated at least 4-fold in Kmt2d^–/– cells (Supplemental Table 4). We searched the literature to identify whether any of the 44 genes had previously been reported as targets of KMT2D, as regulators (particularly inhibitors) of Sox9, and/or as involved in chondrocyte development. No known KMT2D targets were among the 44 genes identified by our transcriptome analysis to be downregulated at least 4-fold in Kmt2d^–/– cells. Three genes were identified as being involved in chondrocyte development and/or endochondral ossification and potentially connected with Sox9 in the literature. Loxl1 encodes lysyl oxidase-like protein 1 (LOXL1) involved in extracellular matrix metabolism and endochondral ossification; however, Loxl1^–/– mice exhibit decreased rather than increased cellularity at growth plates and do not exhibit long bone shortening, as in Kmt2d^+/βGeo mice (22). Ccnd1 encodes cyclin D1, a cell cycle regulator expressed in and required for chondrocyte proliferation (23); expression of this gene has been shown to directly correlate with that of Sox9 (and Col2a1) in an osteoarthritis model, disfavoring a role for Ccnd1 in Sox9 inhibition (24). The product of only 1 gene identified in our analysis was known to play a role in endochondral ossification and had the potential to negatively regulate Sox9 based on previous studies — the transcription factor SHOX2.

Shox2 is the most intriguing candidate gene downregulated over 4-fold (Supplemental Table 4) in our genome-wide transcriptome analysis comparing Kmt2d^–/– cells with Kmt2d^+/+ cells in differentiated chondrocytes (Figure 5A and Supplemental Figure 6A) and undifferentiated cells (Figure 5B and Supplemental Figure 6B). Shox2 is the mouse paralogue to the human SHOX gene, which, when disrupted, causes short stature (25, 26). Shox2 is important in mice for long bone growth and chondrocyte development (27, 28), and importantly, conditional deletion in chondrocytes leads to precocious differentiation with concomitant increase in Col2a1 and Sox9 expression (27), similar to what we observed in Kmt2d^–/– cells (Figure 4, C and E). We confirmed that Shox2 is downregulated over 4-fold in Kmt2d^–/– cells compared with Kmt2d^+/+ cells in the undifferentiated state and in chondrocytes by qPCR in 3 independent experiments (Figure 5E). We subsequently went back to our mouse model to verify these results in vivo; we observed decreased SHOX2 within proximal tibia growth plates in Kmt2d^+/βGeo mice compared with Kmt2d^+/+ littermates (Figure 5, F and G).

Functional testing supports the idea that Shox2 plays a mechanistic role in the cellular differentiation phenotype. Based on these observations, we hypothesized that Shox2 may be a target of KMT2D and a mediator of the precocious differentiation observed in Kmt2d^–/– cells and in Kmt2d^+/βGeo mice. KMT2D places the H3K4me3 activating mark at gene promoters (29), and we previously showed depletion of this mark in the hippocampi of Kmt2d^+/βGeo mice (14). Moreover, examination of ENCODE data from mouse limb bud revealed a strong peak of H3K4me3 at the 5′ end of the Shox2 gene (30, 31) (Supplemental Figure 7A). We therefore performed ChIP-qPCR using an antibody that specifically recognizes H3K4me3 in Kmt2d^–/– and Kmt2d^+/+ cells. We observed a significant depletion of H3K4me3 at multiple sites within this region in Kmt2d^–/– cells compared with Kmt2d^+/+ cells (Figure 6A and Supplemental Figure 7, B–F). In contrast, H3 levels were no different between any of the cell lines (Supplemental Figure 7, B–F). This suggests that, when present, KMT2D may bind to the Shox2 promoter and place the H3K4me3 activating mark in these cells. However, when KMT2D is absent, less of this mark is placed. Reduction of this mark (as opposed to complete loss) makes sense, given the redundancy of the H3K4 histone methyltransferase system. Specifically, reduction of the mark at the Shox2 promoter in Kmt2d^–/– cells, along with decreased Shox2 expression, implicates KMT2D and H3K4me3 in the activation of Shox2 transcription, providing mechanistic insight into precocious chondrocyte differentiation and the epigenetic regulation thereof.

Figure 6

Functional testing supports Shox2 as a target of KMT2D-mediated H3K4me3 and partial mediator of precocious chondrocyte differentiation in KS1. (A) ChIP-qPCR revealed decreased H3K4me3 along the Shox2 promoter in Kmt2d^–/– cells compared with Kmt2d^+/+ cells (mean ± SEM; n = 3 per group per time point; 1-sided unpaired Student’s t test; *P < 0.05). Overexpression of Shox2 in Kmt2d^–/– cells led to partial recovery (downregulation) of (B) Sox9 expression but did not lead to significant recovery of (C) Col2a1 expression or (D and E) Alcian blue staining. The experiment was performed in duplicate (n = 2 per group). Scale bars: 1 mm (in top panels) and 200 μM (in bottom panels). (F) Model for molecular pathogenesis of KS1 growth retardation. TSS, transcription start site; Ctl, control Lentiviral-ORF particles.

Next, we wanted to determine whether overexpression of Shox2 in Kmt2d^–/– cells was sufficient to inhibit precocious chondrocyte differentiation. Transduction of Kmt2d^+/+ and Kmt2d^–/– cells with lentivirus overexpressing Myc-tagged Shox2 led to increased Shox2 transcript levels and increased Myc-tagged SHOX2 protein in 1 Kmt2d^+/+ and 1 Kmt2d^–/– cell line (Supplemental Figure 8, A and B). We first examined Sox9 transcript levels to determine whether overexpression of Shox2 could restore proper Sox9 expression in Kmt2d^–/– cells. Similar to previous observations in nontransduced cells, Kmt2d^–/– cells transduced with control virus showed higher Sox9 transcript levels compared with Kmt2d^+/+ cells transduced with control virus. Notably, upon transduction of Kmt2d^–/– cells with Myc-tagged Shox2, we observed reduced Sox9 transcript levels (Figure 6B), indicating that overexpression of Shox2 is sufficient to inhibit and partially restore Sox9 expression to near WT levels and supports a key role for SOX9 in precocious chondrocyte differentiation in KS1.

To determine if overexpression of Shox2 can fully restore the cellular differentiation phenotype in Kmt2d^–/– cells, we looked at Col2a1 expression and Alcian blue staining. Similar to previous observations in nontransduced cells, Kmt2d^–/– cells transduced with control virus showed higher Col2a1 transcript levels and Alcian blue staining compared with Kmt2d^+/+ cells transduced with control virus (Figure 6, C–E). Upon transduction of Kmt2d^–/– cells with Myc-tagged Shox2, we did not observe a significant difference in Col2a1 expression or Alcian blue staining between Kmt2d^–/– cells overexpressing Shox2 and Kmt2d^–/– cells transduced with control virus, nor did we observe a significant difference in these parameters when we compared Kmt2d^–/– cells overexpressing Shox2 and Kmt2d^+/+ cells transduced with control virus (Figure 6, C–E). Rather, we observed intermediate Col2a1 expression and Alcian blue staining upon overexpression of Shox2 in Kmt2d^–/– cells (Figure 6, C–E). These findings suggest an intermediate differentiation state between Kmt2d^+/+ and Kmt2d^–/– cells and, moreover, that additional cellular targets of KMT2D and of SOX9 contribute to the precocious chondrocyte differentiation observed in KS1 cells. Based on these data, our model for precocious chondrocyte differentiation–mediated growth retardation in KS shows that deficiency of Kmt2d-mediated H3K4me3 leads to impaired activation of Shox2, among other targets (Figure 6F). This allows for Sox9 overexpression and subsequent precocious chondrocyte differentiation via activation of Col2a1, Col10a1, and other chondrocyte-specific genes (Figure 6F), leading to the observed skeletal phenotype associated with KS.

We show for the first time to our knowledge that Kmt2d^+/βGeo mice exhibit a highly specific and quantitative skeletal growth retardation phenotype consisting of decreased body length and weight, shortened long bones, and ventral bowing, dorsal expansion, and brachycephaly of skulls. Our findings resemble those observed classically in individuals with KS, including postnatal short stature and flattened facial profile, and are further supported by more recent detailed observations that patients with KMT2D mutations have relatively increased cranial heights and disproportionately shortened tibias (32). This suggests that our Kmt2d^+/βGeo mice create a good model of the human KS1 skeletal growth retardation phenotype and that our findings provide relevant outcome measures for preclinical therapeutic trials in mice and future clinical studies in humans.

This pattern of skeletal features can be explained, at least in part, by a common mechanism of disrupted endochondral ossification at distinct sites in Kmt2d^+/βGeo mice, namely at long bone growth plates and intrasphenoidal synchondroses within the cranial base. Others have observed similar and additional craniofacial features in a mouse model of the closely related disorder KS2 due to constitutive and neural crest–specific ablation of the H3K27 demethylase Kdm6a/Utx (33). Our findings and theirs could indicate that KS1 and KS2 represent distinct disorders with distinct mechanisms at play. However, we favor the idea that common mechanisms may be present in KS1 and KS2, and the seemingly disparate findings are simply due to timing and experimental design. Shpargel et al. examined the explicitly neural crest–derived viscerocranium/facial skeleton during embryonic development (33), while our studies focus on the ongoing development of the cranial base in the adult animal. However, both likely influence one another, and the impacts of each will be difficult to parse out. Perhaps, simultaneous ongoing decreased cartilaginous growth of the cranial base as suggested here in KS1 — which leads to ventral bowing, dorsal expansion, and brachycephaly of the skull — might contribute to or exacerbate a neural crest–derived craniofacial phenotype previously observed in KS2 (33) but also likely present in KS1. Additional studies will be required to fully understand the pathogenesis of the craniofacial phenotype in KS1 and KS2.

We observed sex-based phenotypic differences of unclear etiology in Kmt2d^+/βGeo mice. KMT2D (formerly MLL, MLL2) promotes estrogen receptor α–mediated (ERα-mediated) gene expression and increased H3K4me3 in cancer cells (34), and knockdown of KMT2D leads to decreased activation of estrogen-responsive genes (35, 36), with less known about interactions between KMT2D and ERα in skeletal tissues. Based on these observations, if an estrogen effect is responsible for the sex-based differences in our mice, deficiency of Kmt2d would be expected to lead to reduced ERα-mediated gene expression. In addition, we would expect Kmt2d^+/βGeo mice to phenocopy ERα-deficient mice and, thus, for only females (not males) to have shortened long bones (37), which was not the case. In addition, we observed opposite effects on trabecular bone parameters compared with those observed in both male (38) and female (37) ERα^–/– mice. These observations suggest that sex-based differences in KS1 mice are not solely estrogen dependent; the mechanism is likely more complex and will require additional studies, beyond the scope of this manuscript, to fully elucidate.

Here, we implicate SOX9 in the pathogenesis of the KS1 skeletal phenotype. Our initially unexpected observations of growth plate expansion in the setting of decreased growth at distinct sites in KS1 mice suggested disrupted endochondral ossification as a mechanism and the chondrocyte as a relevant cell type. A more dynamic system of chondrocyte development, stable ATDC5 cell lines with Kmt2d loss-of-function mutations, confirmed precocious chondrocyte differentiation at the cellular and gene expression levels, thereby suggesting a cell-autonomous phenotype reconstituted in an isogenic model. Our findings point to a mechanism of loss of KMT2D-mediated release of inhibition of Sox9, leading to unchecked and precocious chondrocyte differentiation. Supporting our findings, Sox9 is a well-known mediator of chondrocyte differentiation required for cartilage formation (39–41). Monoallelic loss-of-function mutations in humans lead to campomelic dysplasia, a severe chondrodysplasia (42, 43), and heterozygous ablation in mice resembles the human phenotype (40). Conditional ablation in differentiated chondrocytes disrupts endochondral ossification, resulting in severely shortened growth plate heights, and reveals that Sox9 expression is necessary for chondrocyte proliferation and late stages of hypertrophy (19). This is consistent with our findings of precocious chondrocyte differentiation (hypertrophy) and subtly increased proliferation in Kmt2d^–/– stable cell lines in the setting of increased expression of Sox9. SOX9 is known to directly activate multiple chondrocyte-specific factors, including early and late markers of chondrogenesis — such as Col2a1 and Col10a1, respectively (19, 39) — which are also dysregulated and may be a measure of SOX9’s downstream effects, although it is unlikely that these are the only genes dysregulated secondary to increased Sox9 expression. Moreover, loss of KMT2D leading to release of Sox9 inhibition and subsequent activation of early and late effectors of chondrogenic differentiation would explain our observed increases in proliferative zone, hypertrophic zone, and overall growth plate heights in Kmt2d^+/βGeo mice.

In addition to facilitating chondrocyte differentiation, SOX9 also prevents osteoblast differentiation (19). Specifically, abnormal persistent Sox9 expression prevents transdifferentiation of terminal hypertrophic chondrocytes into osteoblasts (44); however, this does not directly address the contribution of BM mesenchymal stem cell–derived osteoblasts. The net effect of this imbalance between precocious chondrocyte differentiation and reduced osteoblast differentiation may lead to growth retardation in KS1, specifically to shortening of the long bones and decreased bone formation in Kmt2d^+/βGeo mice. However, future conditional studies will be required to reveal the exact impact in individual lineages.

We used genome-wide transcriptome profiling by RNA-seq to identify KMT2D target genes and associated molecular pathways that may play a role in KS1 pathogenesis. We identified Shox2 as a potentially novel target of KMT2D. Shox2 was downregulated over 4-fold in Kmt2d^–/– chondrocytes and undifferentiated mesenchymal cells. Depletion of the H3K4me3 activating mark placed by KMT2D at the Shox2 promoter in Kmt2d^–/– cells supports this new finding as a direct interaction between KMT2D and the Shox2 promoter. We propose that loss of KMT2D (and thus H3K4me3) in chondrocytes and their precursors leads to decreased Shox2 expression and subsequent release of Sox9 inhibition; this leads to increased SOX9, which activates the chondrocyte differentiation program prematurely (including overexpression of Col2a1 and Col10a1) and leads to precocious differentiation of chondrocytes, while inhibiting chondrocyte transdifferentiation to osteoblasts.

Our finding of Shox2 as a potential direct target of KMT2D is intriguing and quite relevant to the KS1 phenotype. A transcription factor highly conserved between mouse and humans (45, 46), Shox2 is expressed in developing limbs, the palate, the CNS, and the heart in both species (45), as well as in additional pharyngeal arch–derived craniofacial structures and the nasal process in humans. Notably, all of these tissues are affected in KS1 with its associated findings of developmental delay/intellectual disability, short stature/lower limb shortening, congenital heart disease, and well-documented craniofacial phenotype, including flattening of the facial profile due to short columella and depressed nasal tip, as well as ear and palate abnormalities and hearing loss (5, 47). No congenital disorders in humans have yet been attributed to SHOX2 intragenic mutations to determine whether this might lead to a phenocopy of KS. Contiguous gene deletions and duplications involving human chromosome 3q25–3q26 including Shox2 and/or its known regulatory elements have been described, however, and these individuals have intellectual disability, skeletal and growth abnormalities, and dysmorphic facial features (48–50), resembling aspects of the KS phenotype. Neural-specific conditional ablation of Shox2 leads to poor cerebellar development associated with precocious differentiation of neural progenitors and motor coordination deficits (51), supporting a functional role for SHOX2 in neurological development and potentially implicating it in motor delays, a key aspect of the KS1 phenotype. Yu et al. (52) showed that Shox2^–/– mouse embryos die during midgestation and exhibit cleft palate and heart defects, features often seen in KS, though it is unclear whether additional skeletal growth phenotypes were examined. Conditional ablation of Shox2 in mesenchymal cells and early chondrocytes leads to shortening of femurs and humeri (27), resembling the femur shortening observed in Kmt2d^+/βGeo mice. This was associated with precocious chondrocyte differentiation and hypertrophy, similar to what we observed in Kmt2d^–/– stable chondrocyte cell lines and at Kmt2d^+/βGeo growth plates (27). Moreover, the conditional ablation in mesenchymal cells at an early time point is associated with increased expression of Col2a1, an early marker of chondrocyte differentiation, and Sox9, the master regulator of chondrogenesis (27), similar to our findings in Kmt2d^–/– and Kmt2d^ΔR5551/– stable chondrocyte cell lines. Finally, later conditional ablation in developing chondrocytes results in increased expression of late markers of chondrocyte differentiation like Col10a1 (27), similar to our findings in Kmt2d^–/– cells. These results support our observations, and together, the findings suggest a mechanism whereby loss of Shox2 expression in chondrocytes or their precursors, whether genetic or epigenetic in nature, leads to Sox9-mediated precocious chondrocyte differentiation and, ultimately, to shortening of long bones.

Transient overexpression of Shox2 in Kmt2d^–/– cells reduced Sox9 transcript levels back toward those observed in Kmt2d^+/+ cells. The reduction in Sox9 appears insufficient to restore proper expression levels of Col2a1 (and Col10a1) and rescue the cellular phenotype, as evidenced by the lack of reduction in Alcian blue staining upon overexpression of Shox2 in Kmt2d^–/– cells. While we observed decreased Alcian blue staining in some images, this was not consistent. This discrepancy may be due to the transient nature of Shox2 overexpression in our system and the extended length of time required to observe changes in Alcian blue staining over the course of differentiation (14 days). Therefore, even a slight decrease in Alcian blue staining may be suggestive of partial restoration of the cellular differentiation phenotype. Alternatively, there are many other target genes and pathways involved in KS-associated precocious chondrocyte differentiation. Indeed, we observed thousands of differentially expressed genes upon differentiation of mesenchymal cells to chondrocytes and roughly a thousand differentially expressed genes when we compared Kmt2d^–/– and Kmt2d^+/+ cells in both the differentiated and undifferentiated states. Restoration of 1 KMT2D target gene, even if critically important, simply may not be sufficient to rescue the cellular differentiation phenotype. Our findings here support the idea that multiple key cell type–specific target genes are disrupted and lead to relevant phenotypes in KS and other Mendelian disorders of the epigenetic machinery (1).

It is worth mentioning that murine Shox2 (mShox2), our potentially novel KMT2D target, is the most closely related gene to the better-known human SHOX (hSHOX) gene, which is known to cause multiple short stature syndromes with associated mesomelia (25, 26, 53) and to contribute to the short stature seen in Turner syndrome (25, 54). Moreover, a subset of individuals with Turner syndrome have a KS phenotype (55–58). SHOX and SHOX2 are thought to exhibit functional redundancy in humans with no known mouse SHOX homologue (45, 46), suggesting that SHOX2 may perform the overlapping functions of both genes in mice (46). In humans, SHOX is thought to mainly pattern the distal limb segment (tibia/radius), while SHOX2 may be more important in patterning the proximal limb segment (femur/humerus; ref. 46). Because both proximal femurs and distal tibias are shortened in Kmt2d^+/βGeo mice and SHOX2 staining is reduced in proximal tibia growth plates in Kmt2d^+/βGeo mice, our studies fit with prior work and suggest that SHOX2 may indeed serve dual functions. Based on our molecular studies in mice and data from humans with KS1, we hypothesize that SHOX2 may be a target of KMT2D in humans, as well, and we argue that SHOX2 is a more likely target than SHOX based on 2 observations. First, a hallmark of every SHOX-associated disorder is mesomelic upper limb shortening, and humans with KS1 were shown to have disproportionately long arms (32), disfavoring SHOX as a key target of KMT2D and mediator of skeletal growth retardation in human KS1. Second, the expression patterns of mShox2 and hSHOX2 are more similar than mShox2 and SHOX, and the tissue distribution of expression of hSHOX2 overlaps almost exactly with the key tissues affected in KS1, including skeletal and others. Because SHOX is expressed in a subset of those tissues, we cannot rule out its involvement; however, our findings favor SHOX2 as a mediator of key aspects of KS1.

Our findings reveal the first mechanistic insights to our knowledge into the molecular basis of the skeletal growth retardation phenotype and the highly characteristic facial appearance of KS1. We have characterized and quantified a robust skeletal growth retardation phenotype in KS1 mice, which resembles the human condition and involves disrupted endochondral ossification at distinct sites. Furthermore, we have elucidated a mechanism involving loss of KMT2D-mediated H3K4me3 at a target (Shox2) and release of SOX9 inhibition, which allows precocious chondrocyte differentiation to proceed unchecked. Partial rescue of the KS1 chondrocyte gene expression profile through modulation of a single disrupted target gene, Shox2, fits with our model that KMT2D acts on multiple targets, some of which have key functional consequences. This is encouraging for future development of therapeutics directed against key epigenetic machinery targets implicated in KS1 and related disorders but suggests that a broader therapeutic approach may be required in some cases. The previously identified roles of SHOX2 in skeletal tissues, including long bones and the palate, as well as in other tissues relevant to KS, like cerebellum, point to additional potential roles for this gene in the pathogenesis of KS1 and related disorders and suggest additional therapeutic benefits to targeting SHOX2 for the treatment of diverse manifestations of KS, other Mendelian disorders of the epigenetic machinery, and more common disease states that disrupt normal growth and development.

Supplemental Methods are available online with this article.

Mice. Kmt2d^+/βGeo mice were previously described (14, 15). Mice were housed in a clean, specific pathogen-free state-of-the-art animal facility in ventilated racks and provided ad libitum access to a standard rodent diet and to filtered water via an automatic watering system. All mice used in these studies were approximately 6 weeks of age unless otherwise noted, and littermates were used as controls in all experiments. Analyses were performed in a blinded manner unless otherwise noted. A single Kmt2d^+/+ mouse was excluded from analyses due to being obviously runted and for having gross midfacial asymmetry related to incisor malocclusion, making this animal an outlier with respect to size, craniofacial shape/asymmetry, and micro-CT parameters. Otherwise, no outliers were excluded from any comparisons or statistical calculation.

High-resolution micro-CT. Femurs, tibias, and skulls were fixed in 4% paraformaldehyde, washed, and transferred to 70% ethanol. High-resolution images were acquired using a desktop micro-tomographic imaging system (Skyscan 1172, Bruker) in accordance with the recommendations of the American Society for Bone and Mineral Research (ASBMR; ref. 59). Full-length scans were reconstructed with NRecon software (Bruker). Femurs and tibias were scanned at 65 keV and 153 µA using a 1.0-mm aluminum filter with an isotropic voxel size of 10 μm. Femur trabecular bone parameters were assessed in a region of interest 500 μm proximal to the growth plate and extending for 2 mm (200 CT slices) using CtAn software (Bruker). Femur cortical bone structure was assessed using a 500 μm region of interest centered on the middiaphysis. Skulls were scanned at an isotropic voxel size of 10–28 μm at 80 keV and 120 μA.

Craniofacial morphometric analysis. We used geometric morphometrics to test the effect of the Kmt2d^+/βGeo mutation on craniofacial structure. Micro-CT image volumes of the head were obtained for a comparative adult sample of Kmt2d^+/βGeo mice (n = 13) and Kmt2d^+/+ littermates (n = 21), as described above (scan parameters, 80 keV and 120 μA; reconstruction, 10–28 μm cubic voxels). Images were reconstructed using Amira postprocessing software (v. 6.1.1, FEI), and 3-D models of the skull were extracted based on density thresholds. The shape of the cranium was estimated by collecting 3-D coordinate data in Amira for biologically relevant, homologous landmarks (K=18; Figure 1E). We used MorphoJ software to produce a Procrustes superimposition of all landmark configurations (60) and, based on the assumption of object symmetry for the skull, analyzed only the symmetric portion of variance to minimize the effects of subtle asymmetries (61). PCA was applied to visualize the most influential patterns of shape variance within the combined sample, and 2-way Procrustes ANOVA was used to test the effect of the mutation on overall cranial shape.

Histology and histomorphometry. Tibias and skulls were fixed in 4% paraformaldehyde and transferred to 10%–14% EDTA for decalcification. Decalcified tibias and skulls were then processed, embedded, sectioned, and stained with H&E. Images were taken using a Nikon 80i microscope and analyzed using NIS elements software. For tibias, longitudinal growth plate sections were used, and proliferative zone, hypertrophic zone, and total growth plate heights were measured in at least 3 sites per section in 4 sections per mouse within the central part of the growth plate. The number of cells per column within the proliferative zone and within the hypertrophic zone were also counted. For skulls, parasagittal sections were used, and intrasphenoidal synchondrosis growth plate heights were measured at 5 sites per section and in 3 sections per mouse within the central two-thirds of the growth plate. For all, an average measurement or cell count per mouse was generated, and then the mean within each experimental group was calculated.

Immunofluorescence. Paraffin-embedded proximal tibial growth plate sections were deparaffinized and rehydrated, followed by antigen retrieval with Pepsin Antigen Retriever (Sigma-Aldrich, R2238) and 0.2 N HCl sequentially at 37°C for 10 minutes. The sections were then permeabilized with 0.3% Triton X-100 (Fisher BioReagents) detergent in PBS for 30 minutes and blocked using 5% donkey serum in PBS for 1 hour before applying primary antibody (SHOX2 1:100, Origene, TA333910) at 4°C overnight. Secondary antibody (Alexa Fluor 594–conjugated donkey anti-rabbit; Thermo Fisher Scientific; A32754) was applied for 1 hour, followed by DAPI staining. Quantification was performed using fluorescence intensity (total fluorescence intensity with background intensity subtraction) per DAPI-stained cell in the growth plate using 2 sections per mouse. Total fluorescence intensity and background intensity were employed in ImageJ. An image-based tool for counting nuclei in ImageJ was used to count the total DAPI-stained cell number in the growth plate.

Cell culture. ATDC5 cell line (18) (Sigma-Aldrich; 99072806) was obtained from European Collection of Authenticated Cell Cultures. ATDC5 cells were maintained in standard medium: DMEM/Ham’s F-12 medium (Thermo Fisher Scientific) containing 5% FBS (MilliporeSigma), antibiotics (100 units/mL penicillin, 100 μg/mL streptomycin; Corning), and 2 mM L-glutamine (Corning) in 5% CO₂ in a 37°C incubator (Thermo Fisher Heracell CO₂ incubator).

Generation of stable cell lines. ATDC5 cells were seeded at 2 × 10⁵ cells/well in 6-well plates overnight. CRISPR-Cas9 constructs were a gift from J. Robertson and L. Goff (Johns Hopkins University, Baltimore, Maryland, USA). pSpCas9(BB)-2A-Puro (PX459) V2.0 vector containing gRNA inserts targeting Kmt2d exon 51 (5′-TCTGGCTCGTTCG CGTATCC-3′) and exon 53 (5′-TCCTTTGGGGATTCGCCGGC-3′) or empty vector were transfected using Lipofectamine 3000 (Invitrogen) according to the manufacturer’s instructions. Twenty-four hours after transfection, the cells were treated with puromycin at a final concentration of 5 μg/mL for 3 days, and the cells were allowed to recover. For single cell clonal analysis, cells were trypsinized and plated in 96-well plates at average 1 cell/well or plated in a 100-mm cell culture dish at an average of 168 or 336 cells/dish and incubated at 37°C for 2 weeks. Each well or plate was microscopically evaluated, and single cell–derived clones were selected, expanded, and genotyped using colony PCR and Sanger sequencing. PCR conditions were 1 cycle of 10 minutes at 94°C, followed by 35 cycles of 20 seconds at 94°C, 30 seconds at 60°C, 1 minute at 72°C, and finishing with a 5-minute incubation at 72°C (primers: 5′-ACTCCAAGTCATCTCAGTAC-3′ and 5′-ACTGATAGTCATAGGTCAGC-3′). Kmt2d^–/– stable cell lines have biallelic deletions within the catalytic SET domain, and Kmt2d^ΔR5551/– stable cell lines each have a monoallelic deletion within the SET domain on 1 allele and deletion of a single amino acid Arg5551 on the other allele. Arg5551 has not been associated with disease and corresponds to the site in the mRNA transcript expected to be cut by Cas9 based on the targeting strategy. KMT2D protein expression was detected by Western blot using anti-KMT2D antibody (MilliporeSigma; ABE1867).

Chondrocyte differentiation. ATDC5 stable cell lines (Kmt2d^+/+, Kmt2d^ΔR5551/–, Kmt2d^–/–), which were generated in house, were seeded in the above medium at 1 × 10⁵ cells/well in 6-well plates for Alcian blue staining and RNA isolation. The next day and every 2–3 days thereafter, medium was replaced with chondrogenic differentiation medium: DMEM/Ham’s F-12 supplemented with 5% FBS, 100 units/mL penicillin, 100 μg/mL streptomycin, 2 mM L-glutamine, 1× Insulin-Transferrin-Selenium (Thermo Fisher Scientific), 50 μg/mL L-ascorbic acid (MilliporeSigma), and 10 mM β-glycerophosphate (MilliporeSigma). Cells were harvested from independent cultures at 0, 4, 7, 14, and 21 days.

Alcian blue staining. At the indicated time points, cells were fixed in 4% paraformaldehyde, stained with Alcian blue (Sigma-Aldrich), and microscopically evaluated. After permeabilization with 1% SDS, absorbance was measured (605 nm) using a BioTek Synergy 2 Multi-Mode Microplate Reader.

qPCR. Total RNA was isolated using TRIZOL reagent according to the manufacturer’s instructions (Invitrogen) at the indicated time points. Reverse transcription was performed using the SuperScript IV First-Strand Synthesis System (Thermo Fisher Scientific). qPCR analysis using the comparative Ct method was performed on the Applied Biosystem Vii 7 system (Thermo Fisher Scientific) with the Powerup SYBR Master Mix (Thermo Fisher Scientific) according to the manufacturer’s instructions. Gapdh was used as reference gene for normalization.

ChIP. ChIP assay was performed according to Cold Spring Harbor protocol (62). Briefly, cells were cross-linked for 10 minutes with formaldehyde at a 1% final concentration, and the reaction was quenched by adding glycine to a final concentration of 137.5 mM. The cells were then rinsed in ice-cold PBS, scraped, and resuspended in a lysis buffer with the addition of complete protease inhibitor cocktail (Cell Signaling Technology). Chromatin was sheared by sonication 3 × 5 cycles of 30 seconds on/30 seconds off with high intensity (Bioruptor) to generate DNA fragments between 300–500 bp, used for immunoprecipitation with anti-H3K4me3 (MilliporeSigma; 07-473) and then captured by protein G-agarose/salmon sperm DNA (MilliporeSigma). Precipitated DNA was reverse cross-linked and then amplified by qPCR using primers amplifying amplicons with the following midpoints relative to the transcription start site of Shox2: −500, −250, 0 (transcription start site; TSS), +200, and +525 and compared with the amount of input DNA before immunoprecipitation.

RNA-seq. ATDC5 stable cell lines (Kmt2d^+/+, Kmt2d^–/–) were seeded as above and allowed to adhere overnight in standard medium. The next day and every 2–3 days thereafter, medium was replaced with chondrogenic differentiation medium. Total RNA was isolated using TRIZOL reagent and RNA Clean & Concentrator-5 kit (Zymo Research) according to the manufacturers’ instructions. Contaminating genomic DNA was removed by treatment with DNase-I (Zymo Research). RNA was assessed for quantity and quality using Qubit RNA BR Assay Kit (Thermo Fisher Scientific) and RNA 6000 Nano Kit on the Agilent Bioanalyzer 2100 system (Agilent Technologies), respectively, according to the manufacturers’ instructions. mRNA was purified from 1 μg of total RNA (≥100 ng/μL) using NEBNext Poly(A) mRNA Magnetic Isolation Module for Illumina (New England Biolabs). Sequencing libraries were generated using NEBNext Ultra II RNA Library Prep with Sample Purification Beads for Illumina (New England Biolabs) in accordance with the manufacturer’s recommendations and were validated using the Agilent High Sensitivity DNA assay on the Agilent Bioanalyzer 2100 system and quantified by NEBNext Library Quant Kit for Illumina (New England Biolabs). After clustering of the index-coded samples, libraries were sequenced on an Illumina HiSeq 2500 platform, which generated 100-bp single-end reads. RNA-seq was performed on 2 (Kmt2d^–/–) or 3 (Kmt2d^+/+) biological replicates for each differentiation state (chondrocytes vs. undifferentiated cells), and each biological replicate is a distinct clonal cell line. Two technical replicates were performed for each cell line.

RNA-seq bioinformatics analysis. We first pseudoaligned the reads to a fasta file (Mus_musculus.GRCm38.cdna.all.fa.gz) obtained from Ensembl (http://uswest.ensembl.org/Mus_musculus/Info/Index, version 91, downloaded January 2018), which contained all mouse cDNA sequences, and then performed the quantification of transcript abundances using Salmon (63). Then, we utilized the tximport R package (64) in order to obtain normalized gene-level counts from the transcript abundances. To achieve this, we set the “countsFromAbundance” parameter equal to “lengthScaledTPM”. Subsequently, using the edgeR (65) and limma (66) R packages, we applied a log₂ transformation to the gene-level counts and normalized each sample with the “voom” function in limma, using the effective library size (i.e., the product of the library size and the normalization factors, calculated using the “calcNormFactors” function in edgeR). We then estimated the mean-variance relationship and computed weights for each observation. Since the differential expression analysis included technical replicates for each of 3 WT and 2 mutant clones, we fit a mixed linear model using the function “duplicateCorrelation” from the statmod R package (67) by blocking on the clone to account for the correlation among technical replicates. Finally, we performed the differential analysis with the limma R package with an FDR of 0.05 as the threshold for statistical significance. Prior to performing the PCA, we converted transcript abundances to gene-level counts using the tximport R package by setting the “countsFromAbundance” parameter to the “no” designation. Then, we first used the “vst” function from the DESeq2 R package (68); with the parameter “blind” set to “TRUE”) in order to apply a variance stabilizing transformation to the obtained gene-level counts. Then, without standardizing the resulting expression matrix, we used the 1,000 most variable genes to estimate the principal components. RNA-seq data is available from NCBI GEO (accession number GSE129365).

Lentiviral transduction. Lenti-ORF particles, Shox2 (Myc-DDK–tagged), and control (Myc-DDK-tagged) were purchased from Origene. To transduce the ATDC5 cells, 1.2 × 10⁵ undifferentiated cells/well were seeded in 6-well plates overnight and then transduced with the lentivirus at a multiplicity of infection (MOI) of 5 and with a final polybrene concentration of 8 μg/mL. Transduced cells were then differentiated and analyzed by qPCR and Alcian blue staining as described above.

Statistics. One- or 2-tailed, unpaired Student’s t tests were used to calculate P values in most cases; P < 0.05 indicated statistical significance. For Figures 4, B–E, and Supplemental Figures 5, A and B, we fit a mixed-effects model using the lme4 R package (69), separately within each time point. We modeled the genotype as a 3-level factor (i.e., we did not assume a dose-response relationship between intact KMT2D copies and the expression of genes under interrogation). For Figure 1, D–F, 2-way procrustes ANOVA was used to test the effect of the Kmt2d mutation on overall cranial shape. The model contained a fixed effect for the genotype and a random effect to account for the correlation between the technical replicates of each clone. At time points where the variance of the random effect was estimated equal to 0, we fit a fixed-effects model using the genotype only. Finally, within each time point, we performed the post hoc pairwise comparisons among the 3 different genotypes using the emmeans R package, and we corrected for multiple testing using Tukey’s adjustment method.

Study approval. All experiments using laboratory mice were performed in accordance with the NIH Guide for the Care and Use of Laboratory Animals (National Academies Press, 2011) and were approved by the Johns Hopkins University IACUC and performed in accordance with their guidelines.