Bone morphology is regulated modularly by global and regional genetic programs

The vertebrate skeleton is composed of numerous bones, each displaying a unique shape and size. Yet, despite this morphological diversity, most bones are formed by a common developmental process, namely endochondral ossification (Berendsen and Olsen, 2015; Cervantes-Diaz et al., 2017; Long and Ornitz, 2013; Olsen et al., 2000). During this process, mesenchymal cells derived from the lateral plate mesoderm under the regulation of the transcription factor SRY-box 9 (SOX9) condense and differentiate into chondroprogenitors, forming the cartilaginous anlage of the future bones (Kawakami et al., 2006; Zhao et al., 1997). Next, a cascade of chondrocyte differentiation steps will give rise to growth plates on both proximal and distal ends of the anlage. Subsequently, blood vessels and bone-building cells, termed osteoblasts, invade the cartilage anlage and drive ossification of the cartilaginous template from the mid-shaft, pursuing the progression of the growth plates (Kronenberg, 2003).

Endochondral ossification has been studied extensively for more than two centuries, generating a vast amount of information on this process. Nevertheless, mechanisms that grant each bone with its distinctive shape are missing. A hallmark of the unique morphology of each long bone is the superstructures that protrude from its surface known as bone eminences, such as the greater and deltoid tuberosities, greater and lesser trochanters, etc., and condyles, such as the distal lateral and medial epicondyles of the humerus. One of the main functions of bone superstructures is to provide an attachment point for tendons and ligaments, which transmit force from the contracting muscles to the skeleton (Lessa et al., 2008; McHenry and Corruccini, 1975; Polly, 2007).

Interestingly, it was previously shown that the site where tendon attaches to bone is formed by a unique set of progenitors that co-express Sox9 and scleraxis (Scx) (Blitz et al., 2013; Sugimoto et al., 2013). Moreover, it was demonstrated that cells that contribute to the bony side of the attachment do not descend from the cells that create the bone shaft anlage, but are specified after them. Finally, it was demonstrated that these unique progenitors are specifically regulated by the TGFβ and BMP signaling pathways (Blitz et al., 2009, 2013). In addition to providing a mechanism for the development of tendon-to-bone attachment site, these findings also offer an alternative, modular model for long bone development. According to this model, one set of Sox9⁺ cells forms the cylindrical anlage of the future bone shaft, which will serve as the bone substructure, whereas a second pool of Sox9⁺/Scx⁺ cells will be added onto this substructure to give rise to the different superstructures.

Current literature provides ample evidence in support of the modular model of long bone development, which is based on lineage tracing, temporal initiation and cell differentiation of the Sox9⁺/Scx⁺ progenitor cells. However, this model is still missing the mechanisms that regulate the early patterning of the Sox9⁺/Scx⁺ progenitor cells, such that superstructures will form on the developing long bone at the right location, and with the correct shape and size. In this study, we examine the regulatory mechanisms underlying these aspects of modularity in skeletal development. Three-dimensional (3D) reconstruction of different long bones and their protruding superstructures allowed us to demonstrate modularity in long bone development. Comparative transcriptomic analysis and loss-of-function assays showed that GLI-Kruppel family member (Gli3) globally regulates the spatial organization of the Sox9⁺/Scx⁺ progenitors in both forelimb and hindlimb, whereas homeobox a11 (Hoxa11) and Hoxd11 regulate patterning of distal superstructures and pre-B cell leukemia homeobox 1 (Pbx1) and Pbx2 regulate proximal superstructure progenitors. Overall, we provide cellular evidence for the existence of a patterning mechanism that involves both global and regional regulation and highlight some of the genes that facilitate the patterning of superstructures along the developing long bones.

To understand better the scope of modularity in long bone morphogenesis, we performed pulse-chase cell-lineage experiments by crossing either Sox9-CreER^T2 or collagen type II alpha 1 (Col2a1)-CreER^T2 transgenic reporter mice with Rosa26-tdTomato reporter mice (Madisen et al., 2010; Nakamura et al., 2006; Soeda et al., 2010). Pregnant females were administered a single dose of 0.03 mg/g tamoxifen/body weight at either embryonic day (E) 10.5 or E11.5. Previously, we showed that this time window allows the exclusive labeling of cells of the cylindrical bone substructure, leaving the superstructure cells unlabeled (described in detail by Blitz et al., 2013; Eyal et al., 2015). To map comprehensively the unlabeled superstructure cells in different skeletal elements, we have established a 3D imaging pipeline that includes clearing of the labeled limbs optically, light-sheet fluorescence microscopy and, finally, reconstruction of obtained images into a 3D object (Treweek et al., 2015; Yang et al., 2014). To visualize the superstructures better, whole-mount limbs were immunostained for COL2A1.

As seen in Fig. 1A,B (and Movies 1-3), in E14.5 Sox9-tdTomato⁺ or E15.5 Col2a1-tdTomato⁺ limbs, all superstructures of different long bones were tdTomato negative, including greater, lesser and deltoid tuberosities, olecranon, and greater, lesser and third trochanters. Interestingly, in addition to bone eminences, we observed various condyles and sesamoid bones that were also tdTomato negative, including distal medial and lateral humeral epicondyles, medial tibial condyle, patella, and medial and lateral fabella. These results demonstrate the generality of the modular process of long bone morphogenesis in the limb.

To demonstrate this modularity further, we examined sagittal and coronal sections of E13.5 fore- and hindlimbs from ScxGFP transgenic embryos (Pryce et al., 2007) that were stained using antibodies against SOX9. As seen in Fig. 1C, we observed an abundance of SOX9⁺ and Scx⁺ double-positive cells at anatomical locations of various future superstructures along the shafts of different bones. Taken together, these results indicate the extensive contribution of Sox9⁺/Scx⁺ superstructure progenitors to the 3D morphology of cartilaginous anlagen, highlighting the modularity of this process both temporally and from lineage analyses.

An immediate implication of the modular model is the existence of a mechanism that regulates the patterning of the Sox9⁺/Scx⁺ progenitors to position the superstructure correctly along the bone substructure. To uncover this mechanism, we first searched the literature for mutations leading to observable superstructure patterning abnormalities. We found that the skeletons of Gli3-null embryos, previously known as extra toes (Xt), exhibit a variety of superstructure abnormalities (Johnson, 1967). To characterize these defects, we crossed Gli3 heterozygous mice and examined skeletal preparations from E17.5 embryos. As seen in Fig. 2A-I, both fore- and hindlimbs of Gli3^null embryos exhibited numerous dysplastic superstructures, including lesser and deltoid tuberosity, olecranon, patella, and additional ectopic sesamoids. Notably, the aberrant patterning of the deltoid tuberosity ranged in severity from mildly abnormal morphology to complete separation from the humeral shaft (compare left and right forelimbs in Fig. 2F and 2G, respectively) (Table S1). We also examined Prx1-Cre;Gli3^fl/fl embryos, in which Gli3 was conditionally knocked out (cKO) in limb mesenchyme (Blaess et al., 2008; Logan et al., 2002). Similar to Gli3^null mutants, limb-specific Gli3 cKO also resulted in diverse abnormally patterned superstructures, such as induction of supernumerary sesamoids around the knee, dysplastic medial tibial condyle and dysplastic ulnar olecranon and coronoid processes (Fig. 2J-O).

To study the possibility that Gli3 regulates the patterning of superstructure progenitors, we sought to examine the distribution of Sox9⁺/Scx⁺ cells in Gli3^null limbs. For that, we crossed heterozygous Gli3 mice with Gli3^+/−; ScxGFP transgenic reporter mice and stained sagittal sections of E13.5 fore- and hindlimbs against SOX9 to highlight chondrocytes. As seen in Fig. 3A-F, the spatial distribution of Sox9⁺/Scx⁺ progenitors in Gli3^null limbs was abnormal. Whereas in control limbs these progenitors were patterned in juxtaposition to the humeral shaft (Fig. 3A, dotted line), in the mutant these cells were laterally spread and less condensed (Fig. 3D, dotted line). To understand how the abnormal distribution of Sox9⁺/Scx⁺ progenitors translates into superstructure abnormality, we examined the spatiotemporal differentiation of deltoid tuberosity precursors in E13.5 and E16.5 sagittal sections from control and Gli3 KO forelimbs. Sections were stained using antibodies against SOX9 and COL2A1, a marker for differentiated chondrocytes (Fig. 3G-L). As expected, at E13.5 we repeatedly observed the lateral spreading of deltoid tuberosity precursors (Fig. 3H,I, dotted lines). Moreover, we noticed that the degree of cellular lateralization varied both in pattern and in distance from the humeral shaft (Fig. 3H,I, white bars). At E16.5, the abnormal deltoid tuberosity morphologies were clearly evident and ranged from dysplastic to separate deltoid tuberosity (Fig. 3K and 3L, respectively). We concluded that Gli3 is necessary for correct spatial organization of deltoid tuberosity Sox9⁺/Scx⁺ precursors and, thus, for superstructure patterning.

The finding that Gli3 regulates superstructure patterning raises the question whether Gli3 acts autonomously within the Sox9⁺/Scx⁺ progenitors. To directly examine whether Gli3 operates autonomously within the Sox9⁺/Scx⁺ progenitor population, we blocked the expression of Gli3 in Scx- or Sox9-expressing cells by crossing Scx-Cre or Sox9-Cre mice to Gli3^fl/fl mice. As seen in Fig. 4, skeletal preparation of both Scx-Cre;Gli3^fl/fl (Fig. 4A-B′) and Sox9-Cre;Gli3^fl/fl cKOs (Fig. 4C-D′) revealed hypoplasia of the deltoid tuberosity and lateral fabella. In addition, we observed dysplasia of the patella in Sox9-Cre;Gli3^fl/fl mutants (Fig. 4D′). This finding clearly indicates an autonomous role of Gli3 in superstructure patterning; yet, the milder phenotypes that we observed in the skeletons of both Scx-Cre;Gli3^fl/fl and Sox9-Cre;Gli3^fl/fl cKOs relative to the Gli3^null mice suggest a potential additional, non-autonomous effect of Gli3 in superstructure patterning.

Although Gli3 had a global effect on superstructure patterning, we suspected that regional regulatory circuitry also exists, such that patterning of specific superstructures could be regulated independently. To search for genes potentially involved in this mechanism, we developed a protocol for prospective isolation of Sox9-tdTomato⁺/ScxGFP⁺ progenitors by fluorescence-activated cell sorting (FACS). Moreover, to enable distinction between global and regional patterning regulators, we designed our system to allow separate isolation of proximal and distal cell populations. To this end, we crossed Sox9-CreER^T2;Rosa26-tdTomato reporter mice with ScxGFP transgenic reporter mice. Pregnant females were administered a single dose of 0.03 mg/g tamoxifen/body weight at E12.0 (Fig. 5A). Next, we harvested E13.5 Sox9-tdTomato⁺/ScxGFP⁺ embryos and micro-dissected their forelimbs at the shoulder level, above the presumable deltoid tuberosity, and at flanking sides of the elbow (Fig. 5Bi,Bii). Finally, the dissected segments were collected into individual tubes, homogenized and then prospectively isolated by FACS (Fig. 5Biii). For each proximal and distal segment, we collected three biological repeats from three separate litters. Each sample contained 7-9% Sox9-tdTomato⁺/ScxGFP⁺ progenitors of total cell population, of which 10,000 cells were collected for transcriptome analysis (Fig. 5C). Further information regarding controls, gate settings and instrument configuration is found in Fig. 5D and in Materials and Methods.

Comparison between the transcriptomes of proximal and distal Sox9-tdTomato⁺/ScxGFP⁺ progenitors revealed 561 differentially expressed genes with at least two-fold change in expression levels between segments (Table S2). We further analyzed 40 of these genes that were considered most likely to be involved, based on the difference in expression levels, enrichment in relevant biological function and known roles in embryonic development (Fig. 5E,F). Notably, 16 of these 40 candidate genes (listed in Table 1) have previously been implicated in superstructure development, providing a strong validation for the effectiveness of our strategy.

Our transcriptome analysis revealed that Hoxa11 and Hoxd11 were highly expressed in distal Sox9-tdTomato⁺/ScxGFP⁺ progenitors and scarcely expressed in proximal cells. Intriguingly, a previous study has demonstrated that a Hoxa11/Hoxd11 compound mutation (Hox11^aadd) resulted in diverse distal forelimb abnormalities, such as the formation of a detached olecranon resembling a sesamoid bone, whereas the proximal limb segment, including the deltoid tuberosity, was not affected (Koyama et al., 2010). Nevertheless, the mechanism underlying the abnormal olecranon development remained to be elucidated. Therefore, we proceeded to study this process in Hox11^aadd mutant embryos.

To that end, we crossed double-heterozygous Hox11^AaDd mice and harvested embryos at E13.5 and E17.5 (Wellik and Capecchi, 2003). Examination of skeletal preparations from E17.5 control and mutant embryos validated the abnormal development of the olecranon (Fig. 6A,A′). Next, we validated the differential expression of Hoxd11 by fluorescent in situ hybridization (FISH) assay on sagittal sections of E13.5 forelimbs. In agreement with our transcriptome analysis results, Hoxd11 was highly expressed distally in proximity to the elbow, but scarcely expressed in proximal parts of the limb (Fig. 6B). Furthermore, Hoxd11 expression was downregulated in Hox11^aadd mutants (Fig. 6B′). The differential expression of Hoxa11 at this stage has been previously demonstrated (Swinehart et al., 2013).

To study the regulatory roles of Hoxa11 and Hoxd11 in distal patterning of superstructure progenitors, we examined the patterning of the olecranon in control and mutant embryos. For this, we crossed double-heterozygous Hox11^AaDd mice and stained forelimb sections from E13.5 double-homozygous Hox11^aadd embryos using antibodies against SOX9 and COL2A1, such that the substructure cells were expected to be Sox9⁺/Col2a1⁺ and the undifferentiated superstructure progenitors to be Sox9⁺/Col2a1⁻. Results showed that the olecranon progenitors failed to organize in the typical pattern observed in control limbs (Fig. 6C-E′). Specifically, the cells were scattered in a proximolateral direction away from the ulnar shaft (Fig. 6C,C′; dotted lines). Sections from E17.5 Hox11^aadd mutants illustrated that the olecranon precursors differentiated at a distance from the ulnar shaft, creating a distinct skeletal element resembling a sesamoid bone (Fig. 6F,F′). These results suggest that Hoxa11 and Hoxd11 are necessary for correct localization of the olecranon progenitors and, thus, for their patterning. Moreover, the proper patterning of superstructure progenitors on the proximal side supports our hypothesis that Hoxa11 and Hoxd11 act as regional distal regulators (Fig. S1A,A′).

Another promising candidate gene from our transcriptome analysis was Pbx1. In contrast to Hoxa11 and Hoxd11, Pbx1 was highly expressed in proximal Sox9-tdTomato⁺/ScxGFP⁺ progenitors compared with distal cells. Indeed, Pbx1^null mutation has been shown to affect only proximal elements, such as the deltoid tuberosity, whereas the distal limb segment, including the olecranon, was not affected (Selleri et al., 2001). We therefore performed the same analysis on Pbx1^null;ScxGFP mutants. After validation of the deltoid tuberosity phenotype in skeletal preparations of E15.0 control and mutant embryos (Fig. 6G,G′), we analyzed Pbx1 expression by immunostaining sagittal sections of E13.5 forelimbs against PBX1 and COL2A1. In agreement with the transcriptome analysis results, PBX1 was highly expressed proximally in proximity to the deltoid tuberosity and was less prominent in the distal forelimb (Fig. 6H). As expected, PBX1 expression was completely ablated in Pbx1^null mutants (Fig. 6H′).

Next, forelimb sections from Pbx1^null;ScxGFP⁺ E13.5 embryos were stained using antibodies against SOX9. As seen in Fig. 6I-K, the Sox9⁺/Scx⁺ deltoid tuberosity precursors were specified in both control and mutant limbs; yet, in the mutant they were organized abnormally and scattered laterally. Finally, examination of sections immunostained against SOX9 and COL2A1 at E15.5 showed that in Pbx1^null mutants, the deltoid tuberosity precursors differentiated into a distinct sesamoid-like cartilaginous element parallel to the humeral shaft (Fig. 6L,L′). These results indicate that Pbx1 is necessary for correct localization of the Sox9⁺/Scx⁺ deltoid tuberosity precursors and, thus, for their patterning. Importantly, the patterning of superstructure progenitors at the distal region was independent of Pbx1 regulation (Fig. S1B,B′).

Together, these results indicate that superstructure patterning is regulated regionally through early modulation of the spatial organization of Sox9⁺/Scx⁺ progenitors by Hoxa11, Hoxd11 and Pbx1.

Pbx2 is a paralog of Pbx1 that is expressed throughout the limb mesenchyme. Previously, it was shown that Pbx1 and Pbx2 act in a dosage-dependent manner during proximal limb development (Capellini, 2006). This led us to examine whether Pbx2 co-regulates proximal superstructure patterning with Pbx1. To avoid early embryonic lethality, we ablated one allele of Pbx2 on the background of a limb-specific knockout of Pbx1 (Prx1-Cre;Pbx1^fl/fl;Pbx2^+/−) and compared the phenotype with both wild-type (Prx1-Cre⁻;Pbx1^+/+;Pbx2^+/+) and Pbx1 cKO (Prx1-Cre;Pbx1^fl/fl) embryos (Ficara et al., 2008; Selleri et al., 2004). To this end, sections from E13.5, E15.5 and E16.5 control and mutant forelimbs were immunostained using antibodies against SOX9 and COL2A1 (Fig. 7A-C″). At E13.5, we observed abnormal patterning of deltoid tuberosity precursors in both Prx1-Cre;Pbx1^fl/fl and Prx1-Cre;Pbx^fl/fl;Pbx2^+/− mutant embryos. The precursors were laterally scattered (Fig. 7A′,A″), comparable to the deltoid tuberosity precursors of Pbx1^null limbs (Fig. 6I′). By E15.5, the deltoid tuberosity precursors of Prx1-Cre;Pbx1^fl/fl embryos had fully differentiated, whereas in Prx1-Cre;Pbx1^fl/fl;Pbx2^+/− we observed delayed differentiation (Fig. 7B-B″). Notably, in Prx1-Cre;Pbx1^fl/fl embryos the humerus-deltoid tuberosity boundary was populated by SOX9⁺ cells that had not been detected in Prx1-Cre;Pbx1^fl/fl;Pbx2^+/− embryos. Moreover, the gap between the humeral shaft and deltoid tuberosity precursors was roughly twice as large in Prx1-Cre;Pbx1^fl/fl;Pbx2^+/− limbs than in Prx1-Cre;Pbx1^fl/fl limbs. Finally, by E16.5, the deltoid tuberosity of Prx1-Cre;Pbx1^fl/fl embryos has attached to the humeral shaft, whereas the deltoid tuberosity of Prx1-Cre;Pbx1^fl/fl;Pbx2^+/− mutants remained separated (Fig. 7C-C″), as in Pbx1^null mutants (Fig. 6L′). These diverging morphologies were further validated by skeletal preparations made of limbs taken from each mutants (Fig. 7D-D″). Taken together, these results suggest that PBX genes act together in regulating the patterning of proximal superstructures, possibly by regulating the level of lateralization in a dose-dependent manner.

The finding that superstructure patterning can be regulated both globally and regionally led us to hypothesize that these two programs interact with one another. To examine this possibility, we produced a compound mutant mouse carrying limb-specific knockout of Pbx1 and Gli3. Specifically, on the background of a limb-specific knockout of Pbx1 (Prx1-Cre;Pbx1^fl/fl) we ablated either one or two alleles of Gli3 (Prx1-Cre;Pbx1^fl/fl;Gli3^fl/+ or Prx1-Cre;Pbx1^fl/fl;Gli3^fl/fl, respectively). Then, sections from E13.5 and E16.5 control and mutant forelimbs were immunostained using antibodies against SOX9 and COL2A1 (Fig. 8A-B″). At E13.5, deltoid tuberosity precursors in the Prx1-Cre;Pbx1^fl/fl;Gli3^fl/fl embryos were laterally scattered and noticeably less condensed than in Prx1-Cre;Pbx1^fl/fl;Gli3^fl/+ embryos (Fig. 8A-A″). The phenotypic difference between genotypes became more pronounced by E16.5. Whereas ablation of a single Gli3 allele resulted in a dysplastic deltoid tuberosity that was attached to the humeral shaft, ablation of both alleles resulted in a detached deltoid tuberosity (Fig. 8B-B″). These diverging morphologies were further validated by skeletal preparations of limbs from each mutant group (Fig. 8C-C″). Importantly, whereas individual limb-specific knockout of either Prx1-Cre;Gli3^fl/fl (data not shown) or Prx1-Cre;Pbx1^fl/fl (Fig. 7D′) resulted in dysplastic deltoid tuberosity, the compound Prx1-Cre;Pbx1^fl/fl;Gli3^fl/fl double knockout mutant displayed a detached deltoid tuberosity, as in both Gli3^null (Fig. 2G) and Pbx1^null (Fig. 6G′) embryos. These results suggest that, similar to the compound Pbx1/Pbx2 mutation, genetic interactions also occurred between Pbx1 and Gli3 in a dose-dependent manner, highlighting the possibility of interactions between the global and regional regulatory programs.

The abnormal superstructure development we observe in the different mutant mouse models could be the result of several possible mechanisms that affect the Sox9⁺/Scx⁺ progenitors, including patterning, proliferation, differentiation and cell death. Having already demonstrated abnormal patterning of the Sox9⁺/Scx⁺ progenitors in the different mutants, we next sought to study whether proliferation rate contributed to the observed phenotype. To this end, we analyzed incorporation of 5-bromo-2-deoxyuridine (BrdU) into Sox9⁺/Scx⁺ cells, focusing on the hypoplastic deltoid tuberosity. Gli3 and Prx1-Cre;Pbx1;Pbx2 pregnant mice were injected subcutaneously with BrdU. Incorporation of BrdU was allowed for a period of 2 h, after which the pregnant mice were sacrificed. Incorporated BrdU was detected by immunostaining using antibodies against BrdU and SOX9 on sections taken from E13.5 Gli3 and Prx1-Cre;Pbx1;Pbx2 control and mutant forelimbs. Next, the total SOX9⁺ and SOX9⁺/BrdU⁺ cells were counted in the region of the greater and deltoid tuberosities and the fraction of proliferating cells was calculated for both control and mutant embryos (Fig. S2). Our results showed that in Gli3 mutants, the center of proliferation appeared to shift inwards from the lateral side of the superstructures to the main shaft (Fig. 9A′,B′; yellow arrows demarcate center of proliferation). However, quantification showed no significant difference in the fraction of proliferating cells. In contrast, we observed a mild but significant reduction of proliferating cells in both Prx1-Cre;Pbx1^floxed and Prx1-Cre;Pbx1^floxed;Pbx2^het mutant embryos (Fig. 9D′,E′) compared with control embryos (Fig. 9C′) (Fig. S2).

The detached deltoid tuberosity could be the consequence of apoptosis. To examine the possible involvement of cell death in this phenotype, we performed a terminal deoxynucleotidyl transferase dUTP nick end labeling (TUNEL) assay on sections taken from E13.5 Gli3 and Prx1-Cre;Pbx1;Pbx2 control and mutant forelimbs. As a control, we examined the interdigital space, which undergoes massive cell death at E13.5 (Fig. 9F, inset). Results showed no cell death within the superstructure precursors in either control or mutant embryos (Fig. 9F-J; superstructures are demarcated by dotted lines).

Previously, we showed that BMP4 is necessary for the differentiation of superstructure progenitors (Blitz et al., 2009). To address the possibility that BMP signaling is reduced in our mutants, we examined the activation of SMAD1/5/8 (Retting et al., 2009) in the Sox9⁺/Scx⁺ progenitors. Sections taken from E13.5 Gli3 and Prx1-Cre;Pbx1;Pbx2 control and mutant forelimbs that were immunostained using antibodies against pSMAD1/5/8 and COL2A1 revealed different patterns of pSmad1/5/8 activity in control and mutant embryos (Fig. 9K-O), yet the abnormal patterns correlated with the mispatterning of the Sox9⁺/Scx⁺ progenitors in these mutants. This finding suggests that the aberrant Smad activation is secondary to the progenitor mispatterning. Overall, we concluded that although cell death did not contribute to the observed phenotypes, reduction in cell proliferation might have had some effect.

The musculoskeletal system acts as a system of levers and pulleys to create locomotion. Thus, one way to achieve variation in locomotor strategies is by changing the location of the connection between lever and pulley. Shifting the position of a bone superstructure, to which a muscle is connected by a tendon, along the bone shaft modifies the pulling force vector of that given muscle and, thereby, facilitates different types of locomotive capabilities. Variations in the positioning of superstructures and their effect on locomotion have been well documented (Archer et al., 2011; Milne and O'Higgins, 2012; Polly, 2007; Salton and Sargis, 2008, 2009). However, the mechanism that produces these variations was unknown.

An early step in the development of endochondral bones is mesenchymal condensation and the formation of the cartilaginous template that prefigures the future bone. Here, we identify components of the genetic program that regulate the subsequent step of this process, namely the patterning and condensation of superstructure precursors. Moreover, we demonstrate that the mechanism of superstructure patterning involves both global and regional modules, highlighting the modularity in long bone development. Specifically, by combining comparative transcriptomic analysis with cross-reference to existing literature, we identified a list of candidate genes that might be involved in this genetic program. Further analysis of several of these candidate genes, namely, Gli3, Pbx1, Pbx2, Hoxa11 and Hoxd11, established their involvement in regulating superstructure patterning. Importantly, Gli3 perturbation affected superstructures throughout different parts of the limb skeleton, whereas mutations in Pbx1 and Pbx2 or in Hoxa11 and Hoxd11 affected either proximal or distal superstructures, respectively. Interestingly, we found that the global and regional components of the genetic program are interconnected, as compound mutations of the proximal regulator Pbx1 and the global regulator Gli3 led to increased phenotype severity, compared with mutations in either gene alone.

In addition to these genes, our transcriptomic analysis identified 12 more candidate genes that were reportedly involved in long bone development. Interestingly, many of these genes regulate either proximal or distal limb segments, or alternatively, establish the proximodistal axis itself. For example, we identified additional HOX genes, such as Hoxa9, Hoxd9, Hoxa10 and Hoxd10, which were shown to regulate the patterning of proximal superstructures, such as the deltoid tuberosity (Wellik and Capecchi, 2003; Xu and Wellik, 2011). Other candidate genes, such as Shox2 or Meis1, operate locally either as downstream effectors of HOX genes (Neufeld et al., 2014) or as co-factors of proximal regulators, such as Pbx1 (Mercader et al., 2009), respectively. We also identified genes that control the formation of the proximodistal axis, such as Cyp26B1, which modulates the levels of retinoic acid (Yashiro et al., 2004) or Irx3, Irx5 and Hand2, which modulate the activity of the Shh/Gli3 signaling pathway (Galli et al., 2010; Li et al., 2014). Together, these reports further support the modular model for superstructure patterning and suggest the involvement of at least two hierarchical programs in this process. Higher is the program that governs the establishment of the proximodistal axis and, further downstream, local programs adjust the development of each individual segment. Finally, whereas our transcriptomic analysis was designed to highlight differentially expressed genes along the proximodistal axis, it lacked information regarding globally expressed genes, which could potentially regulate superstructure patterning, such as Gli3, which was expressed ubiquitously in both proximal and distal Sox9⁺/Scx⁺ progenitors. Additional evidence for the existence of other global regulatory genes is given by a recent report demonstrating that Tbx3 globally regulates patterning of forelimb superstructures, namely the greater and deltoid tuberosity and olecranon (Colasanto et al., 2016).

Although our work sheds light on the genetic mechanism that regulates superstructure patterning, the cellular mechanism underlying the patterning process is still missing. A clue for its nature is given by the observations that in Gli3, Pbx and Hox11 mutant embryos, superstructure precursors spread abnormally away from the developing bone substructure and display delayed differentiation. This may imply the existence of a mechanism controlling the migration of Sox9⁺/Scx⁺ progenitors to a specific condensation site. Once this program is perturbed, cells migrate to the wrong site, where they form an abnormal bony element. Alternatively, it is possible that Sox9⁺/Scx⁺ chondroprogenitors are specified from selected cells already present at the designated superstructure locations. In that case, the mechanism would involve activation of a chondrogenic program in specific cell subpopulations by extrinsic or intrinsic signals at specific spatiotemporal positions. In either scenario, it appears that the mechanism that decides where superstructure development should take place involves the Gli3, Pbx and Hox genes. Our finding that ablation of Gli3 specifically in Sox9⁺ and/or Scx⁺ cell lineage resulted in hypoplasia of selected superstructures clearly suggests that Gli3 regulates the Sox9⁺/Scx⁺ progenitors autonomously. However, the differences we observed in the severity of the phenotype between the Gli3-null and the different conditional alleles may suggest a non-autonomous contribution as well. An alternative explanation for the differences in severity is differences in efficiency between the various Cre lines in targeting the Gli3 locus.

A previous work in which the expression of both Shh and Gli3 was blocked proposed that the function of Shh and Gli3 in limb skeletal patterning is limited to refining autopodial morphology (Litingtung et al., 2002). Our work clearly expands the involvement of Gli3 in skeletal patterning by demonstrating its role as a global regulator of superstructure patterning. Interestingly, both Shh loss of function and Shh/Gli3 double knockout resulted in apparently normal deltoid tuberosity, suggesting that the effect of Gli3, at least in the proximal domain of the limb, is independent of Shh.

The patterning of superstructures involves control over their size. Our finding of reduced progenitor cell proliferation in Prx1-Cre;Pbx1^fl/fl and Prx1-Cre;Pbx1^fl/fl;Pbx2^+/− mutant embryos might suggest that the genetic program regulates superstructure morphology by coordinating patterning and cell proliferation. In contrast, we present evidence that size control does not involve apoptosis of Sox9⁺/Scx⁺ progenitors.

Another level of coordination that is required is between the skeleton and the musculature. Interestingly, several reports have demonstrated that mispatterning of specific muscles is coupled with aberrant superstructures. For example, ablation of Tbx3 results in abnormal patterning of stylopod musculature, which attaches to dysplastic olecranon and deltoid tuberosity (Colasanto et al., 2016), whereas in Hox11^Aadd or Hox11^aaDd compound mutants, abnormal patterning of zeugopod musculature is accompanied by abnormal olecranon patterning (Swinehart et al., 2013). Taken together, these results strongly suggest that muscle and superstructure patterning may be regulated in a coordinated manner.

Previously, we reported that sesamoid bones, which are small auxiliary bones that are highly variant in anatomical location, size and number, also originate from Sox9- and Scx-positive precursors (Eyal et al., 2015; Eyal et al., 2019). Interestingly, we demonstrated that the Sox9- and Scx-positive precursors of some sesamoid bones, such as the lateral fabella, were patterned such that they developed away from the long bone substructure (Eyal et al., 2019). Those results are in line with our current evidence showing the formation of sesamoid-like bones following perturbation in either Gli3, PBX or Hox11 genes. Moreover, they further highlight the degree of plasticity stemming from the modular model, which offers an evolutionary mechanism to induce not only local variations in bone morphology by patterning the superstructures, but also formation of new auxiliary bones without having to rewrite the entire skeletogenic program. We therefore propose that bone superstructures and sesamoid bones be placed under one category of cartilage elements that share a common origin and developmental program.

To conclude, our results provide insight into the formation, patterning and development of long bone superstructures and their impact on bone morphology. Furthermore, they provide strong evidence in support of the modular model of skeletogenesis and demonstrate the level of modularity during long bone morphogenesis in terms of cellular origins, genetic regulation, morphology and skeletal assembly of superstructures, namely their ability to attach or detach from the substructure. Specifically, we show that whereas Sox9⁺ progenitors establish the bone substructure (Fig. 10A), superstructure anlagen are formed and assembled modularly onto the substructure by Sox9⁺/Scx⁺ progenitors (Fig. 10B). We show that these unique Sox9⁺/Scx⁺ progenitors contribute not only to bone eminences but also to various condyles and sesamoid bones and that their patterning involves both global and regional regulatory modules that include Gli3, Pbx and Hox genes (Fig. 10B). Following their specification and patterning, these superstructures differentiate and become integral to the developing long bones (Fig. 10C,D). Importantly, by performing a comparative transcriptomic analysis, we were able to highlight additional candidate genes that may be implicated in superstructure patterning. Further studies of these molecular players will provide a better understanding of the regulatory network at play in superstructure formation and long bone morphogenesis.

All experiments involving mice were approved by the Institutional Animal Care and Use Committee (IACUC) of the Weizmann Institute. For all timed pregnancies, plug date was defined as E0.5. For harvesting of embryos, timed-pregnant females were euthanized by cervical dislocation. The gravid uterus was dissected out and suspended in a bath of cold PBS and the embryos were harvested after removal of the placenta. Tail genomic DNA was used for genotyping by PCR. All experiments were performed in at least three biological repeats, i.e. on three embryos from separate litters.

Sox9-CreER^T2 mice (Soeda et al., 2010) were received from the laboratory of Haruhiko Akiyama, Kyoto University, Kyoto, Japan. ScxGFP,Scx-Cre and Sox9-Cre transgenic mice were obtained from Ronen Schweitzer, Shriners Hospital for Children Research Division, Portland, OR, USA.

The generation of Col2A1-CreER^T2 (Nakamura et al., 2006), Sox9-CreER^T2 (Soeda et al., 2010), ScxGFP (Pryce et al., 2007), Gli3^null (Johnson, 1967), Prx1-Cre (Logan et al., 2002), floxed-Gli3 (Blaess et al., 2008), Pbx1^null (Selleri et al., 2001), floxed-Pbx1 (Ficara et al., 2008), Pbx2^null (Selleri et al., 2004), Hox-a11^null (Hsieh-Li et al., 1995), Hox-d11^null (Davis and Capecchi, 1994) and Rosa26-tdTomato (Madisen et al., 2010) mice has been described previously.

To create Gli3, Hox11^aadd and Pbx1 mutant mice, animals heterozygous for the mutations were crossed; heterozygous embryos were used as a control. To create Prx1-Gli3 or Prx1-Pbx1 mutant mice, floxed-Gli3 or floxed-Pbx1 mice were mated with Prx1-Cre-Gli3 or Prx1-Cre-Pbx1, respectively. To create Sox9-Gli3 or Scx-Gli3 mutant mice, floxed-Gli3 mice were mated with Sox9-Cre-Gli3 or Scx-Cre-Gli3, respectively. To create Prx1-Pbx1^flox-Pbx2^+/− mutant mice, floxed-Pbx1-Pbx2^+/− mice were mated with Prx1-Cre-Pbx1^flox mutant mice. To create Prx1-Pbx1-Gli3 mutant mice, floxed-Pbx1-Gli3 mice were mated with Prx1-Cre-Pbx1-Gli3 mutant mice. As a control, Prx1-Cre-negative embryos were used.

For genetic lineage analysis and FACS experiments, either Col2A1-CreER^T2 or Sox9-CreER^T2 mice were crossed with Rosa26-tdTomato reporter mice. Induction of Cre recombinase was performed at various pregnancy stages by administration of 0.03 mg/g tamoxifen/body weight in corn oil by oral gavage (stock concentration was 5 mg/ml).

Cartilage and bone in whole mouse embryos were visualized, after skinning and disemboweling, by staining with Alcian Blue (Millipore Sigma, A5268) and Alizarin Red S (Millipore Sigma, A5533) and clarification of soft tissue with 3% potassium hydroxide and 100% glycerol (McLeod, 1980).

For preparation of paraffin sections, embryos were harvested at various ages, dissected and fixed in 4% paraformaldehyde (PFA)/PBS at 4°C overnight. After fixation, tissues were dehydrated in an ethanol series to 100% and embedded in paraffin. The embedded tissues were sectioned at a thickness of 7 µm using a microtome (LM2235, Leica) and mounted onto slides.

For preparation of O.C.T-embedded sections, embryos were harvested at various ages, dissected and fixed in 1% PFA/PBS at 4°C overnight. Fixed embryos were then dehydrated in 30% sucrose overnight at 4°C. Next, samples were dissected, soaked in O.C.T (Tissue-Tek) for 30-60 min and then frozen in O.C.T. Frozen samples were sectioned at a thickness of 10 µm using a cryostat (CM3050S, Leica) and mounted onto slides.

FISH on paraffin sections were performed using Hoxd11 digoxigenin (DIG)-labeled RNA probe (Shwartz and Zelzer, 2014). Probe was detected using anti-DIG-POD (1:300, 11207733910, Roche), followed by Cy3-tyramide labeled fluorescent dyes, according to the instructions of the TSA Plus Fluorescent Systems Kit (PerkinElmer). Finally, slides were counter-stained using DAPI (1:1000, D9542, Millipore Sigma).

For immunofluorescence staining for SOX9 and COL2A1, 7 μm-thick paraffin sections of embryo limbs were deparaffinized and rehydrated in water. Antigen was then retrieved in 10 mM citrate buffer (pH 6.0), boiled and cooked for 10 min in a microwave oven. In order to block non-specific binding of immunoglobulin, sections were incubated with 7% goat serum, 1% bovine serum albumin (BSA) dissolved in PBT (PBS+0.1% Tween 20) for 60 min at room temperature. Following blockage, sections were incubated overnight at 4°C with primary anti-SOX9 antibody (1:150, AB5535, Millipore Sigma). Then, sections were washed in PBT and incubated with Cy3-conjugated secondary fluorescent antibodies (1:100, 711-165-215, Jackson ImmunoResearch) for 60 min at room temperature. After staining for SOX9, slides were washed in PBT and fixed in 4% PFA at room temperature for 10 min. Then, slides were incubated with proteinase K (Millipore Sigma, P9290), washed and post-fixed again in 4% PFA. Next, sections were washed and incubated overnight at 4°C with primary anti-COL2A1 antibody (1:50, II-II6B3, the Developmental Studies Hybridoma Bank). The next day, sections were washed in PBT and incubated with Cy2-conjugated secondary fluorescent antibodies (1:100; 715-225-150, Jackson ImmunoResearch) for 60 min at room temperature. Slides were mounted with Immu-mount aqueous-based mounting medium (Thermo Fisher Scientific).

For immunofluorescence staining for SOX9 and ScxGFP, 10 μm-thick cryostat sections of embryo limbs endogenously labeled for ScxGFP were used. SOX9 immunofluorescence staining was performed as described above, but with omission of the antigen retrieval step, using primary SOX9 antibody and secondary Cy3-fluorescent antibodies.

For immunofluorescence staining for PBX1 and COL2A1, paraffin sections were used and antigen retrieval and blockage of non-specific binding were performed as described above. Endogenous peroxidase was quenched by incubation in 2% H₂O₂/PBS for 30 min at room temperature. Staining for PBX1 was performed using primary anti-PBX1 antibodies (1:400, C-4342, Cell Signaling Technology) followed by HRP-conjugated secondary antibodies (1:1000, 711-035-152, Jackson ImmunoResearch) and Cy3-tyramide labeled fluorescent dyes, according to the instructions of the TSA Plus Fluorescent Systems Kit (PerkinElmer). Subsequently, staining against COL2A1 was performed as described above.

For immunofluorescence staining for pSMAD1/5 and COL2A1, 7 μm-thick paraffin sections of embryo limbs were deparaffinized and rehydrated in water. For antigen retrieval, slides were incubated with proteinase K (Millipore Sigma, P9290) for 20 min at room temperature, washed and post-fixed again in 4% PFA. Endogenous peroxidase was quenched by incubation in 2% H₂O₂/PBS for 30 min at room temperature. Next, staining for pSMAD1/5 was performed using primary anti-pSMAD1/5/9 antibodies (1:200, CST-13820, Cell Signaling Technology) followed by biotin-conjugated secondary antibodies (1:100, 711-065-152, Jackson ImmunoResearch) and HRP-conjugated streptavidin (1:200, NEL750001EA, PerkinElmer). Finally, detection was performed using Cy3-tyramide labeled fluorescent dyes, according to the instructions of the TSA Plus Fluorescent Systems Kit (PerkinElmer). Subsequently, staining against COL2A1 was performed as described above.

To examine cell proliferation, pregnant females were injected subcutaneously with 100 µg/g BrdU/body weight dissolved in PBS (stock concentration was 10 mg/ml) (B5002, Millipore Sigma). After 2 h of exposure, pregnant mice were scarified and the embryos were harvested, genotyped, fixated and embedded in paraffin. To detect BrdU-incorporated cells, immunostaining using primary anti-BrdU antibodies (1:200, MCA2060, Bio-Rad) and secondary Cy3-conjugated fluorescent antibodies (1:100; 712-165-153, Jackson ImmunoResearch) was performed on 7 μm-thick paraffin sections. Antigen retrieval using 10 mM citrate buffer (pH 6.0) and counterstaining for SOX9 were performed as described above. Quantification of BrdU-stained cells was performed using sections collected from three embryos from three different litters. From each embryo, at least two sections were used for statistical analysis (n≥6). In each section, the number of Sox9⁺ and Sox9⁺/BrdU⁺ stained cells were counted within the area of the greater and deltoid tuberosities. The fraction of proliferative cells was calculated for each group and the difference between control and mutant embryos was tested using Student's t-test. Statistical significance was determined as P≤0.05.

To examine cell death, we performed a TUNEL assay using ApopTag plus peroxidase in situ apoptosis kit (S7101, Millipore Sigma) according to the user guide supplied by the manufacturer, except for the use of DAB as the main peroxidase substrate. As a substitute, we used Cy3-tyramide labeled fluorescent dyes, according to the instructions of the TSA Plus Fluorescent Systems Kit (PerkinElmer). Finally, counterstaining was performed using DAPI (1:1000, D9542, Millipore Sigma).

For whole-mount imaging, samples were first cleared and immunostained using the PACT-decal technique (Treweek et al., 2015; Yang et al., 2014). Briefly, either Col2A1-CreER^T2- or Sox9-CreER^T2 mice were crossed with Rosa26-tdTomato reporter mice. Following tamoxifen administration at E11.5 (Col2A1-CreER^T2) or E10.5 (Sox9-CreER^T2), as described above, embryos were harvested at E15.5 or E14.5, respectively. Embryos were dissected and fixed in 1% PFA/PBS at 4°C overnight. Next, samples were washed in PBTx (PBS+0.1% Triton X-100+0.01% sodium azide) at room temperature, then embedded into a hydrogel of 4% (wt/vol) acrylamide in 1× PBS with 0.25% thermal initiator 2,2′-azobis[2-(2-imidazolin-2-yl)propane] dihydrochloride (VA-044, FUJIFILM-Wako). The hydrogel was allowed to polymerize at 37°C for 5 h. The samples were removed from the hydrogel, washed in PBTx, and moved to 10% sodium dodecyl sulfate (SDS) with 0.01% sodium azide, shaking at 37°C for 3 days, changing the SDS solution each day. Samples were washed four times with 1× PBTx at room temperature over the course of 24 h and fixed in 4% PFA at room temperature for 10 min. Then, samples were incubated with proteinase K (Millipore Sigma, P9290) for 30 min with shaking at 37°C, washed and post-fixed again in 4% PFA. To detect COL2A1, samples were first incubated with 5% goat serum dissolved in PBTx at 37°C overnight in order to block non-specific binding of immunoglobulin. Next, samples were incubated with primary COL2A1 antibodies (1:150, II-II6B3, the Developmental Studies Hybridoma Bank) in 2.5% goat serum/ PBTx shaking at 37°C for 3 days, changing the antibody solution each day. Samples were washed four times with 1× PBTx at room temperature over the course of 24 h. Next, samples were incubated with secondary Cy3 or Cy2 antibodies (1:150, 715-165-150 or 715-225-150, respectively, Jackson ImmunoResearch) in 2.5% goat serum/PBTx with shaking at 37°C overnight. Samples were washed again with four changes of 1× PBTx, and the refractive index (RI) of the sample was brought to 1.45 by submersion in a refractive index matching solution (RIMS) prepared by dissolving 35 g of Histodenz (Millipore Sigma, D2158) in 0.02 M phosphate buffer (74% wt/vol), and shaking gently at room temperature overnight. Finally, samples were embedded in 1% low gelling Agarose (Millipore Sigma, A9414) in PBS in a glass capillary, submerged in RIMS and stored in the dark at room temperature until imaging.

The cleared samples were imaged with a Zeiss Lightsheet Z.1 microscope. For each limb, a low-resolution image of the entire limb was taken with the 20× Clarity lens at a zoom of 0.36. Light-sheet fusion of images was done if necessary (Zen software, Zeiss). Tile stitching and 3D-image reconstruction was performed using Arivis Vision4D (Arivis) and Imaris (Bitplane) Software.

Flow cytometry analysis and sorting were performed on a BD FACS AriaIII instrument (BD Immunocytometry Systems) equipped with a 488, 407, 561 and 633 nm lasers, using a 70 µm nozzle, controlled by BD FACS Diva software v8.0.1 (BD Biosciences), at the Weizmann Institute of Science Flow Cytometry Core Facility. Further analysis was performed using FlowJo software v10.2 (Tree Star).

For collection of cells, Sox9-CreER^T2-tdTomato;ScxGFP mice were crossed with Rosa26-tdTomato;ScxGFP reporter mice. Embryos were harvested at E13.5 following tamoxifen administration at E12.0, as described above. Forelimbs were dissected and suspended in cold PBS using 15 ml tubes.

To extract cells from tissues, PBS was replaced with 1 ml pre-heated 0.05% trypsin (0.25% Trypsin EDTA solution A, Biological Industries) diluted in DMEM media (Thermo Fisher Scientific) and incubated for 15 min at 37°C, with gentle agitation every 5 min. Tissues were then dissociated by vigorous pipetting using 1 ml tips. Next, 4 ml of DMEM supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin solution (Biological Industries) was added and cell suspensions were filtered into a 15 ml tube using a syringe and a 40 µm filter net. Finally, tubes were centrifuged at 1000 rpm (90 g) for 7 min, supernatant was removed and cells were resuspended in 1 ml of cold PBS and used immediately for FACS.

Single-stained GFP and tdTomato control cells were used for configuration and determining gate boundaries. Live cells were gated by size and granularity using FSC-A versus SSC-A and according to staining with propidium iodide (PI, 1 µg/ml) and DAPI (1 µg/ml). FSC-W versus FSC-A was used to further distinguish single cells. In addition, unstained, GFP-stained only and tdTomato-stained only cells were mixed in various combinations to verify that the analysis excluded false-positive doublets.

Purified RNA from FACS-isolated Sox9⁺/Scx⁺ cell samples was used for library preparation according to a standard MARS-Seq protocol (Jaitin et al., 2014). Library quality was analyzed by a 2200 TapeStation instrument (Agilent Technologies, data not shown). The experiment produced libraries of high quality and sufficient quantity. Libraries were subsequently sequenced by Illumina NextSeq 500.

The analysis was performed as described by Köster and Rahmann (2012). Briefly, a unique molecular identifier (UMI) sequence present in the R2 read was inserted in the read name of R1 sequence file using a python script. Cutadapt (Martin, 2011) was used to trim low quality and poly A/T and adapter sequences (-a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -a “A{100}” -a “T{100}” --times 2 -q 20 -m 30). Sequences were mapped using TopHat2 to mouse genome build mm10 (Kim et al., 2013). UMI information was integrated into the BAM files as tags, using a python script. The BAM file was converted to SAM format using Samtools (Li et al., 2009). Duplicate reads were marked and filtered based on having the same UMI and mapping to the same gene, using a python script. UMI counts per gene were calculated using modified HTSeq-count script (Anders et al., 2014) and a RefSeq gtf file (downloaded from igenomes UCSC), which was modified to contain a window surrounding the 3′ UTR.

DESeq2 (Love et al., 2014) was used for normalization and to detect differentially expressed genes based on negative binomial distribution and a generalized linear regression model. The model used contained two factors: olecranon or deltoid tuberosity Sox9⁺/Scx⁺ precursors and a batch factor. Genes were considered differentially expressed if the difference in expression between at least two sample types was statistically significant after adjustment with the fdrtool package (adjusted P-value≤0.05) (Strimmer, 2008). After filtering based on expression levels (sum normalized count in all samples greater than 10 and maximum expression level higher than 5), clustering of the standardized normalized counts was carried out using click algorithm (Expander package, Ulitsky et al., 2010). Further analysis was performed using GSEA (Broad institute) and Ingenuity (IPA).

We thank N. Konstantin for expert editorial assistance, and members of the Zelzer laboratory for their advice and suggestions; R. Schweitzer for providing the ScxGFP, Scx-Cre and Sox9-Cre mice, L. Selleri for providing the Pbx1^null and Pbx2^null mice, and H. Akiyama for providing the Sox9-CreER^T2 mice; and O. Golani from the Weizmann Institute Department of Life Sciences Core Facilities for her guidance and assistance in analyzing and designing the experiments. Light sheet imaging was made possible thanks to the De Picciotto-Lesser Cell Observatory in memory of Wolf and Ruth Lesser of the MICC. We thank C. Vega from the Weizmann Institute Department of Design, Photography and Printing for designing the graphic model.