A unique stylopod patterning mechanism by Shox2-controlled osteogenesis

The body segments and appendages of vertebrates are patterned by clusters of Hox and paraHox genes that are expressed along the major body axes (Pearson et al., 2005). Mutations or even subtle changes in the expression of Hox genes cause homeotic transformation of the axial skeleton (Kondrashov et al., 2011; Pearson et al., 2005), and reduction or loss of skeletal elements in vertebrate limbs (Kondrashov et al., 2011; Zakany and Duboule, 2007). In contrast to the specific in vivo requirement of Hox gene expression for correct patterning, the in vitro binding specificity of Hox factors to DNA motifs remains controversial, raising the question whether additional machineries such as transcription co-factors exist to ensure the regional specific function of Hox genes (Mann et al., 2009). Indeed, a family of three amino acids loop extension (TALE) homeodomain proteins including Meis and Pbx subclasses has been extensively characterized as DNA binding co-factors for Hox proteins to achieve the DNA binding specificity and form a highly conserved Hox-TALE patterning system with its origin being traced back to ancestral species such as starlet sea anemone (Hudry et al., 2014; Mann et al., 2009; Parker et al., 2011; Slattery et al., 2011). However, it is still under debate whether TALE factors could fully satisfy the in vivo binding specificity of Hox proteins. The recently proposed low-affinity Hox-TALE binding motif clusters on Hox-TALE bound enhancers (Crocker et al., 2015) implies the existence of additional factors to confer sufficient in vivo binding specificity. Alternatively, instead of being the primary binding factor, Hox proteins are known to play an accessory role for the interaction of Meis transcription factors with specific enhancers. Moreover, Meis factors can even function without Hox on a large proportion of these enhancers in branchial arch (BA) patterning (Amin et al., 2015), suggesting that additional tissue-specific transcriptional mechanisms contribute to the in vivo binding specificity of enhancers with Hox and TALE factors. However, whether other co-factors exist for Hox-TALE system so far remains unknown.

In the developing vertebrate limb, bone elements form via endochondral ossification, whereas osteogenesis is preceded by the formation of cartilaginous template with Runx2⁺/Osx⁺ osteogenic precursors initially residing in the perichondium and later on invading into cartilaginous template (Long and Ornitz, 2013). Limb skeletal patterning by Hox-TALE transcriptional complexes was proposed to be essential for endochondral skeletogenesis by directly regulating cartilage template formation (Zakany and Duboule, 2007). Along the proximal to distal (PD) axis, the skeletal elements of tetrapod limbs are patterned by Hox9 to Hox13 located within the HoxA/D gene clusters. Additionally, the expression of TALE factors is also regulated by signaling pathways along the PD axis, in which context the proximal retinoic acid (RA) signaling and the distal FGF signaling antagonistically determine the proximal expression of Meis genes that marks the stylopodial segment and facilitates the nuclear localization of Pbx in the proximal limb (Cunningham and Duester, 2015; Mercader et al., 2000). Together with HoxA/D9 and HoxA/D10, Meis and Pbx provides patterning code for the stylopodial skeleton (Capellini et al., 2011; Cunningham and Duester, 2015; Penkov et al., 2013). Intriguingly, compound deletion of HoxA/D gene clusters produces considerably milder defects in the stylopodial skeleton than that in the distal zeugopodial and autopodial skeletons that are patterned by Hox10-Hox13 (Kmita et al., 2005; Raines et al., 2015), suggesting that the stylopod adopts a unique mechanism for patterning that is less dependent on HoxA/D factors.

We have shown previously that inactivation of Shox2, encoding a paired-like homeodomain transcription factor, causes developmental defects of multiple organs including the heart, palate and limb (Bobick and Cobb, 2012; Cobb et al., 2006; Espinoza-Lewis et al., 2009; Ye et al., 2015a; Yu et al., 2005,, 2007). Strikingly, Shox2 mutation causes loss of the stylopod in both forelimbs and hindlimbs, which was thought to be attributed to the direct function of Shox2 in chondrogenesis (Bobick and Cobb, 2012; Yu et al., 2007). However, an epistatic additive interaction between HoxA/D genes and Shox2 was seen in limb development (Neufeld et al., 2014), suggesting an involvement of Shox2 in the Hox-TALE patterning system.

Here, using our unique Shox2 allelic toolsets, we undertook a comprehensive analysis of Shox2 expression and the fate of Shox2⁺ cells in the developing limb. Hereby, we reveal an unexpected direct role of Shox2 in osteogenesis for stylopodial skeletal patterning. Our ChIP-Seq and RNA-Seq analyses demonstrate that Shox2 functions by directly regulating enhancers that are co-occupied by Hox-TALE factors to specify the stylopod that emerges in the juxtaposition of the trunk with strong Meis and Pbx gene expression and the proximal limb where Shox2 is highly expressed. Moreover, by retrospective and de novo characterization of Shox2-occupied enhancers, we demonstrate that co-occupancy of Shox2 and TALE factors represents a key feature of the enhancer grammar for proximal limb-specific enhancer activity. Our results indicate that Shox2 acts as a repressor on these enhancers and is required for modulation of cell fate choices in limb development.

We have shown previously that Shox2 is expressed in the developing proximal limb and Shox2 deficiency causes severely mispatterned stylopodial skeletal elements (Cobb et al., 2006; Yu et al., 2007). To unravel the functional mechanism of Shox2 in stylopod patterning, we first sought to comprehensively document and analyze Shox2 expression patterns and cell populations that are derived from Shox2⁺ cells. We took advantage of our recently generated Shox2^HA/+ (Wang et al., 2014; Ye et al., 2015b) and Shox2^Cre/+ (Sun et al., 2013) alleles for protein expression analysis and fate mapping.

The Shox2^HA/+ allele represents a knock-in of a Flag-2xHA-Shox2a-IRES-DsRed-pA coding cassette into the endogenous Shox2 locus (Wang et al., 2014; Ye et al., 2015a), which allows live imaging of Shox2 expression via DsRed and Shox2 protein localization using anti-HA antibodies. We chose to analyze Shox2 expression starting from embryonic day (E)12.5 at which stage the distinct anlagen of stylopodial and zeugopodial skeletal elements are formed. Live imaging revealed strong Shox2 expression in the proximal portion of the developing limb at E12.5 (Fig. 1A,B). An intensive Shox2 expression in the proximal limb, especially around the chondrogenic center of the stylopod, was observed by anti-HA staining at the same stage (Fig. 1C). Consistent with broad Shox2 expression at E12.5, fate mapping using Shox2^Cre/+ allele revealed contribution of Shox2⁺ cells to multiple connective tissue cell types, including chondrocytes, osteoblasts, adipocytes, and dermal fibroblasts (Fig. 1D-G). In contrast to the extensive labeling of osteogenic cells in the bone-forming area (Fig. 1E), Shox2^Cre/+ labeled only a portion of chondrocytes (Fig. 1D), consistent with relatively weak real-time Shox2 expression in the chondrogenic center (Fig. 1C).

As the expression analysis and fate-mapping results did not indicate a lineage etiology for the phenotype observed in Shox2 mutants and Shox2-expressing cells contribute extensively to the proximal connective tissue cells, we speculated that virtual loss of the stylopod in Shox2 mutants results from absence of Shox2⁺ cells. We therefore compounded the Rosa^mTmG allele to Shox2^Cre/− mice (null mutant) to examine the contribution of Shox2⁺ cells in the absence of Shox2.

Shox2^Cre/−;Rosa^mTmG mice were analyzed at E13.0, the latest stage the Shox2 deficient embryos appear normal prior to severe complications caused by cardiac defects (Espinoza-Lewis et al., 2009; Ye et al., 2015b). In these embryos, the Shox2⁺ lineage is labeled by mGFP, and the expression pattern of mGFP reporter did not change significantly compared with controls (Fig. 2), indicating that Shox2 is not required for the proliferation and migration of the majority of Shox2⁺ cells in the stylopod. Consistent with the weak and transient expression of Shox2 in the chondrocytes of controls (Fig. 2A-C), the chondrogenic center of the stylopodial skeleton formed in Shox2-deficient Shox2^Cre/−;Rosa^mTmG embryos (Fig. 2D-F). Interestingly, intensive Shox2 expression in the perichondrium of the stylopod of control embryos was observed at this stage (Fig. 2A). Immunostaining on the adjacent sections showed that these cells are also positive for Runx2, a molecular marker for osteoblastic precursors (Fig. 2C). However, in Shox2^Cre/−;Rosa^mTmG embryos, these Shox2⁺/Runx2⁺ cells lost Runx2 expression completely and expressed Sox9 aberrantly, indicating defective fate decision (Fig. 2B,C,E,F).

Given a potential role of the embryonic perichondrium as the osteogenic precursor pool (Colnot et al., 2004) and the extensive contribution of Shox2⁺ cells to osteoblasts/osteocytes (Fig. 1E), we asked whether Shox2⁺ osteogenic cells are derived from Shox2⁺/Runx2⁺ perichondrial cells and if the defective fate decision of early perichondrial Shox2⁺/Runx2⁺ cells could lead to aberrant osteogenic differentiation of Shox2⁺ cells at a later stage. To overcome the cardiac pacemaker defects that cause early embryonic lethality in Shox2 mutants, we made use of Nkx2-5^F/F allele that, in combination with the Shox2^Cre/−;Rosa^mTmG alleles, conditionally rescues the cardiac defects and allows Shox2 mutants to survive to the birth without interfering defects in other organs irrelevant to Nkx2-5 (our unpublished results).

At E14.0, GFP⁺/Osterix⁺ cells could be clearly detected in the inner layer of perichondrium in the stylopod of Shox2^Cre/+;Nkx2-5^F/F;Rosa^mTmG embryos (Fig. 3A). At E17.5, the GFP⁺/Osterix⁺ osteoblast/osteocytes were abundantly localized in the diaphysis of the stylopodial skeleton (Fig. 3B). Together with the restricted Shox2 expression in the perichondrium at E14.5 (Fig. 3E-H) and the complete lack of Shox2 expression in the stylopodial skeleton at E17.5 (data not shown), our observations indicate a direct contribution of the Shox2⁺ perichondrial cells to the osteogenic lineage. However, in Shox2-deficient embryos (Shox2^Cre/−;Nkx2-5^F/F;Rosa^mTmG) at E14.0 and E17.5, the Osterix⁺ layer of perichondrial cells and the GFP⁺/Osterix⁺ osteogenic cells that were otherwise found abundantly in the diaphysis of control mice, were absent (Fig. 3C,D). The requirement of Shox2 for the formation of osteogenic cells appears to occur at the early fate determination phase. This is supported by the observations that at E14.5, Shox2 expression becomes predominantly restricted to the outer layer of the perichondrium that is characterized by relatively low levels of Runx2 expression (Fig. 3E-H) and was thought to have a less definitive role in osteogenesis (Swinehart et al., 2013).

It was proposed that skeletal patterning in the limb by Hox genes is largely exerted by their expression in the perichondrium (Swinehart et al., 2013; Villavicencio-Lorini et al., 2010). Given that Shox2 deficiency causes a specific loss of perichondrial osteogenic cell type signature and leads to a depletion of osteogenic cells, we hypothesized that Shox2 functions to control stylopodial skeleton patterning through maintaining the fate of perichondrial osteogenic precursors that initially provide instructive cues for chondrocyte differentiation, and later differentiate to osteoblasts for bone formation. In supporting this hypothesis, we found that Shox2⁺ cells contribute specifically to the osteogenic population of the stylopodial skeleton but not that of the zeugopodial skeletal elements (Fig. S1), correlating well with the specific loss of the stylopod in Shox2 mutants.

To determine if Shox2 expression in the osteogenic cells is essential for stylopod patterning, we utilized Col3.6-Cre (Liu et al., 2004) to inactivate Shox2 specifically in osteogenic cells (Fig. 3I). In Col3.6-Cre;Shox2^F/F mice, the stylopodial elements of both forelimb and hindlimb were severely shortened (Fig. 3J), mimicking the limb defects observed in Shox2-null mutants (Fig. 3K). Together with the fact that cartilage-specific deletion of Shox2 by Col2a1-Cre did not produce a similar limb phenotype (Fig. 3K; Fig. S1), our results demonstrate a pivotal role for Shox2 in the osteogenic lineage essential for stylopodial skeleton patterning.

To reveal the regulatory molecular framework behind this unique Shox2 function in osteogenic differentiation in the stylopod, we used the Shox2^HA allele that allows efficient immunoprecipitation of endogenous HA-tagged Shox2 (Fig. 4A). We performed anti-HA ChIP-Seq using developing limbs from Shox2^HA/+ and wild-type control embryos at E12.5 and E13.5. Resulting peaks with various levels of fold enrichment were successfully verified by ChIP-qPCR (Fig. S2).

We first assigned the Shox2 binding peaks to the nearby coding genes and assessed the gene ontology (GO). The GO terms for skeletogenesis are highly enriched for nearby genes flanking Shox2 binding peaks, indicating a pivotal role of Shox2 in osteogenesis despite its broad expression in progenitor cells of multiple connective tissue cell types (Fig. 4B). Our analysis showed that Shox2 interacts predominantly with distal regulatory elements that are highly conserved among vertebrates as shown by plotting based on PhastCons (Fig. 4C). We next plotted Shox2 binding signals with the DNaseI hypersensitivity sequencing (DHS) data generated by the Encode project (Yue et al., 2014). Shox2 binding signals and DHS signals from the developing limb (E11.5), and adult heart and lung were categorized based on relative signal intensity using k-mean clustering into four groups. As shown in Fig. 4D, the genome-wide binding intensity of Shox2 at E12.5 and E13.5 is generally similar, suggesting a persistent function of Shox2. Moreover, the majority of Shox2 binding sites colocalize with DHS sites accessible specifically in the embryonic limb, which are enriched for skeletogenic-related GO terms (Fig. S3). Therefore, our identified limb-specific Shox2 target regions likely contain the unique set of Shox2-controlled cis-regulatory modules that is required for the transcriptional control of skeletal patterning genes in the stylopod.

To further characterize if these Shox2-bound elements represent active enhancers, we plotted the published histone modifications H3K4me1, H3K27ac and H3K27me3 ChIP-Seq data obtained from embryonic limbs (Attanasio et al., 2014) versus the Shox2 binding peak summits. Indeed, the active enhancer markers H3K4me1 and H3K27ac are generally enriched around Shox2 binding sites (Fig. 4E). In addition, H3K27me3, which is frequently associated with inactive enhancers, appears to be mildly enriched around summits of Shox2 binding sites, although it is clearly absent from the center of Shox2 binding peaks (Fig. 4E), suggesting that a set of Shox2 binding enhancers is in a poised state. Given the association of Shox2 occupancy with limb-specific DHS sites and enhancer marks, we sought to determine if the occupancy of Shox2 can be used as a feature to isolate limb-specific enhancers. We exploited the vista enhancer database (Visel et al., 2007) and found that only 19% of active enhancers validated through transgenic enhancer assays contain limb-specific activity (Fig. 4F). After filtering the active enhancers using Shox2 occupancy as a criterion, the percentage of the limb-specific active enhancers increased to 50%. Consistent with Shox2 expression in the developing proximal limb, many of the enhancers occupied by Shox2 displayed proximal specific activity in the limb, as exemplified in Fig. 4G. These observations indicate that the enhancers occupied by Shox2 and located near genes involved in skeletogenesis integrate the transcriptional program specifically required for stylopod formation.

We next sought to search for transcription co-factors of Shox2 that co-operate in and underlie its unique function in stylopodial skeletogenesis. We first performed de novo motif discovery on Shox2 binding peaks in the embryonic limb. As shown in Fig. 5, in contrast to results discovered using the same methods in hindlimb ChIP-Seq top peaks of Pitx1 (Infante et al., 2013), which is also a paired-like homeodomain transcription factor, the top motifs discovered in the Shox2 binding sites are highly related to motifs for Hox-TALE factors, including those similar to the Hox9-TALE motif (Huang et al., 2012), Meis motif and Pbx-TALE motif. These analyses suggest an association of Shox2 with Hox-TALE factors.

Given that Hox9, Pbx and Meis are all transcription factors that are highly abundant in the proximal limb (Capellini et al., 2011; Zakany and Duboule, 2007) and act potentially as determinants for accessible chromatin landscape, the outcome of the motif discovery on Shox2-occupied sites may be biased by the availability of transcription factor-accessible sites that are already established by Hox9, Pbx and Meis. To determine the specific association of Shox2 with Hox-TALE factors, we further conducted Shox2 ChIP-Seq on the developing palate where Shox2 is specifically expressed in the anterior domain overlapping with the future bony hard palate (Gu et al., 2008; Yu et al., 2005) and performed motif discovery similarly. The occupancy of Shox2 in the palate shows distal enhancer binding properties and high relevance to skeletogenesis similar to that in limb (Fig. 6A,B). TALE-related motifs were also discovered in the Shox2-bound peaks in the palate (Fig. 5). However, despite the similar association of Shox2 with TALE-Hox factors in the limb and palate, the top peak signals of the limb and palate do not overlap (Fig. 6C), suggesting that Shox2, along with TALE-Hox factors, is involved in different sets of machineries for skeletogenesis in the embryonic limb and palate, respectively.

To understand the collaborative regulatory manner of Shox2 and Hox-TALE, we performed ChIP-Seq of Pbx, a key TALE family homeodomain factor, using E12.5 limb tissue. The top motifs discovered from the Pbx-occupied sites (Fig. 5) are similar to that discovered from the Shox2-occupied sites, and GO analysis indicates a similar function for Pbx in controlling cis-regulatory modules near skeletogenic genes in the developing limb (Fig. 6D,E). Furthermore, the motif for Meis was abundant in the peak center of Pbx ChIP-Seq data, consistent with the fact that Pbx and Meis often exist in the same transcription complex and share the same binding sites (Choe et al., 2009; Moens and Selleri, 2006). Interestingly, a discrepant genome-wide binding pattern between Pbx in the limb and that of Pbx and Meis in the brachial arches (BAs) (Amin et al., 2015) was also discovered (Fig. 6F), indicating that Pbx, similar to Shox2, is involved in a set of transcription programs unique for the skeletal patterning of the limb. A unique genome-wide co-occupation between Shox2 and Pbx was further revealed by heat-map plotting (Fig. 6G) with the binding signals of Hand2 (Osterwalder et al., 2014) and Runx2 (Wu et al., 2014), which are crucial for skeletal formation. Moreover, aggregate plots (Fig. 6H) indicate that the binding peaks of Shox2 colocalize with those of Pbx in the embryonic limb (this study) and Meis (Amin et al., 2015), suggesting that Shox2 and TALE factors can function in the same transcription complex.

Interaction of Shox2 and Pbx binding peaks resulted in about 9000 co-occupied elements that cluster around skeletogenic-related genes (Fig. 6I). Moreover, multiple Shox2 and Pbx co-occupied sites were found around such genes, for example, flanking Tbx18 and Runx2 (Fig. 6J), suggesting that these putative regulatory elements function additively on their putative target gene as a ‘regulatory archipelago’ – a phenomenon that was found for the regulation of Hox genes (Montavon et al., 2011). Remarkably, within the cluster of Shox2-Pbx co-occupied enhancers downstream of Runx2, an enhancer that is located within the same topologically associated domain (TAD) as Runx2 [inferred from published Hi-C data (Dixon et al., 2012)] showed proximal limb specific activity in the transgenic enhancer assay (Fig. 6K). This observation suggests that the regional control of Runx2 expression by a Shox2-TALE-Hox mechanism contributes to stylopodial skeleton patterning.

Given a strict requirement of Shox2 for Runx2 expression in the stylopod and the co-occupancy of Shox2 and Pbx on the regulatory elements of Runx2, we initially hypothesized that Shox2 and TALE-Hox factors would be collectively required for the expression of key genes for stylopod patterning. This notion seemed to be supported by the hypothesis that Meis and Pbx genes are essential for stylopodial skeleton patterning based on their strong proximal limb expression (Capellini et al., 2011; Tabin and Wolpert, 2007). To determine if an epistatic additive output of Shox2 and TALE factors contributes to stylopod patterning, we first analyzed the Shox2⁺ domain in relation to that of Meis gene expression. As shown in Fig. 7A-C, similar to previous reports (Mercader et al., 2000), Meis1/2 is expressed most strongly in the proximal limb adjacent to the trunk, and the expression intensity decreases distally. To our surprise, the expression domain of Meis1/2 appears partially overlapped with, but complementary to the Shox2⁺ domain that increases distally in the stylopod (Fig. 7A-C). A closer examination on the expression of Meis and Pbx proteins and Shox2 in the developing stylopod also identified a partial exclusive expression of both Meis and Pbx proteins from Shox2-expressing cells surrounding the cartilage primordium (Fig. 7D-K). These observations together suggest that the Shox2⁺ perichondrial mesenchymal cells acquire a distinct distal identity within the stylopod, and Shox2 may enforce this identity by repressing the expression of TALE transcription factors. Consistently, Shox2 colocalizes with both Pbx and Meis transcriptional regulators on putative cis-regulatory elements near Meis and Pbx genes (Fig. S4), suggesting that Shox2 exerts its repressive function on Meis and Pbx via a self-regulatory machinery.

To further understand the epigenomic change underlying the unique identity of Shox2⁺ perichondrial mesenchymal cells, we compared the Pbx/Meis ChIP-Seq data obtained from embryonic trunk (Penkov et al., 2013) with our Shox2/Pbx ChIP-Seq data. As summarized in the Venn diagram (Fig. 7L), the Shox2/Pbx binding sites in the limb overlap very little with the Pbx/Meis binding sites in the trunk, consistent with the notion that the Shox2⁺ perichondrial cells have a distinct chromatin landscape accessible to TALE factors.

The cis-regulatory elements for Meis genes have not been characterized in vivo. Although Meis1/2 are expressed in both the trunk and proximal limb, we reasoned that if the occupancy of Shox2/TALE transcription factors underlies the trunk to limb epigenomic change and that Shox2 represses the expression of Meis and Pbx genes directly by acting on their putative enhancers, then the strong co-occupancy of Shox2 and Pbx in the limb and the weak occupancy of Meis and Pbx in the trunk would serve as a principle to identify limb-specific enhancers for Meis. As shown in Fig. 7M, a potential enhancer for Meis1 showed prominent co-occupancy of Shox2 and Pbx in the limb, but mild occupancy by Meis and Pbx in the trunk. Transgenic enhancer assay confirmed its enhancer activity in the developing limb (Fig. 7N). Interestingly, this site does not show identifiable H3K27ac ChIP-Seq signal in the limb (Fig. S4) (Shen et al., 2012), suggesting that co-occupancy by context-specific transcription factors could serve as a more sensitive approach for identifying tissue-specific enhancers. In addition, this enhancer showed H3K27ac enrichment in differentiated osteoblast-like cell line (Fig. S4) (St John et al., 2014), suggesting that it is an osteogenic enhancer for Meis1 and that Shox2 represses this enhancer to prevent Meis1 expression in perichondrial osteogenic lineage cells. In support of this notion, in the absence of Shox2, an elevated expression of Meis was observed in the distal region of the humerus (Fig. 7O), and elevated Meis and Pbx expression was also detected in the Shox2⁺ perichondrial mesenchymal cells (Fig. S5). Shox2 appears to establish a trunk to limb transition status (Fig. 7P). It was reported previously that excessive expression of Pbx (Gordon et al., 2010) or Meis (Mercader et al., 2009) inhibited osteogenesis. Shox2 thus probably functions to antagonize the inhibitory effect of Meis and Pbx for stylopodial skeletogenesis. This hypothesis is also in line with the fact that excessive proximal RA signaling that is upstream of Meis could also cause failure of proper skeletogenesis in the developing limb, including the stylopod (Cunningham and Duester, 2015).

The likely repressive function of Shox2 on the expression of Meis and Pbx genes prompted us to further understand the general function of Shox2 in transcriptional regulation of its putative target genes. We accordingly performed RNA-Seq to profile the transcriptome of GFP⁺ cells sorted from E12.5 Shox2^Cre/−;Rosa^mTmG limb and compared results with those from Shox2^Cre/+;Rosa^mTmG mice. The list of differentially expressed genes (DEGs) was intersected to the list of genes with Shox2-bound peaks in their putative regulatory elements to identify Shox2-targeted genes. Subsequent analysis of the intersected DEGs and gene-peak association by BETA, a program designed for generalizing the function of transcription factors (Wang et al., 2013), showed that Shox2 functions primarily as a repressor (Fig. 8A). Additionally, the same analysis using Shox2/Pbx co-bound peaks as input also revealed a repressive function for Shox2 (data not shown), suggesting that TALE-Hox machinery is also involved in the repressive action of Shox2. GO analysis further demonstrated that the putative Shox2 repressed targets are more relevant to limb development compared with targets activated by Shox2 (Fig. 8B). Indeed, many putative direct repressive targets of Shox2, such as Prrx1, Osr2, Gria2, Epha7, Rsrc1 and Bmp4 (Tables S1 and S2), were found to show elevated/ectopic expression in the absence of Shox2 in our previous studies by microarray and in situ hybridization on E10.5 and E11.5 limbs (Vickerman et al., 2011; Yu et al., 2007). Similar to Runx2 and Tbx18 (Fig. 6J), Shox2-TALE co-occupied sites were also identified as clusters around their associated genes (Fig. S4), supporting the notion of a convergent function of Shox2-Hox-TALE machinery toward their target genes.

We next investigated whether Shox2 represses target enhancers directly by using a transgenic enhancer assay. First, we used an enhancer containing a Shox2-Pbx-Meis co-occupied element downstream of Bmp4, which revealed complementary activity of the transgenic enhancer compared with the expression domains of Shox2 in developing limbs (Fig. 8C), suggesting a repressive role of Shox2 on this enhancer. We next created a stable transgenic line (mm741-lacZ) from a Shox2/Pbx1/Meis co-bound enhancer in intron 3 of the Rsrc1 in the Rsrc1-Shox2 locus (Fig. 8D), given that Rsrc1 is a robust repressive target of Shox2 and that Shox2 also shows self-repressive activity (Fig. S6). Similar to its human ortholog hs741 (Fig. 5G), this enhancer displayed strong limb-specific activity (Fig. 8E). Compounding this transgenic enhancer-reporter allele to Shox2^−/− background caused significant elevation of the enhancer activity in the proximal limb (Fig. 8E; N=3). A direct repressive activity of Shox2 toward the Shox2-Rsrc1 locus is further supported by the co-occupancy of Shox2, Pbx and Meis on an enhancer of Shox2 itself (Rosin et al., 2013) and a strong repressive activity of Shox2 on this enhancer in luciferase reporter assay in vitro (Fig. 8F). Interestingly, many genes we identified as putative repressive targets of Shox2 were shown or known previously to be broadly expressed in mesenchymal progenitor cells of the early developing limb (Vickerman et al., 2011; Yu et al., 2007) and less widely expressed in the osteogenic lineage, raising the possibility that the repressive action of Shox2 helps to steer the progenitor cells away from a non-osteogenic cell fate.

In this study, we identified a Shox2⁺ cell population specifically required for the patterning of the stylopodial skeleton, as summarized in Fig. 9. This population of Shox2⁺ perichondrial mesenchymal cells appears to instruct the maturation and hypertrophic differentiation of chondrocytes in the cartilaginous mode and contributes to the osteogenic lineage at later stage, because in the Shox2^−/− limb, the stylopodial cartilage mode fails to undergo hypertrophy and lacks bone formation (Cobb et al., 2006; Yu et al., 2007), in contrast to the milder limb phenotype in mice carrying cartilage-specific inactivation of Shox2 (Bobick and Cobb, 2012; and this study). In the perichondrial mesenchymal cells, Shox2 acts to determine the fate of osteogenic lineage, as it is required for the expression of key osteogenic genes including Runx2. On the other hand, Shox2 represses a panel of genes, as demonstrated by integrated analysis of our RNA-Seq and ChIP-Seq data. Among the repressed genes are Meis and Pbx, both of which contribute to correct proximal skeleton patterning. Together with the fact that Shox2, Meis and Pbx show broad genome-wide co-occupation on osteogenic genes, the repressive activity of Shox2 toward Meis and Pbx assists to establish a Shox2 dominant form of Shox2/TALE occupancy on Shox2/TALE co-occupied regulatory elements and favors pro-osteogenic transcription output essential for the instructive role of the perichondrial cells towards cartilaginous mode patterning as well as osteogenic cell fate (Fig. 9).

Understanding the grammar for tissue-specific enhancers is essential for deciphering the subregion-specific transcription output. In addition to the limb enhancers characterized by approaches utilizing whole-limb tissue (Visel et al., 2009), the enhancers for the zone of polarizing activity (ZPA) and apical ectodermal ridge (AER) were also characterized by microdissection with H3K27ac ChIP-Seq (VanderMeer et al., 2014). However, the enhancer grammar underlying the PD patterning of the limb and whether the PD specificity can be recapitulated by enhancer elements are yet to be addressed. We took the advantage of the specific expression domain of Shox2 in the proximal limb and its essential function for proximal limb patterning to characterize putative functional enhancers with proximal limb-specific activities. By computational and comparative analyses on our Shox2 ChIP-Seq data obtained from the limb and palate as well as Pbx ChIP-Seq on the limb with existing ChIP-Seq data, we found that the regulatory elements for Shox2/TALE co-occupation probably represent the key features for the proximal limb-specific enhancers. This notion is further supported by a de novo discovered proximal limb enhancer for Meis1 (Meis1^En) based on Shox2/TALE co-occupancy.

Although Meis and Pbx are strongly expressed in the trunk region (Mercader et al., 2000; Penkov et al., 2013) and our study also revealed an expression gradient of Meis from the trunk to the distal end of the stylopod, the Meis1^En identified in this study is weakly occupied by Meis and Pbx in the trunk region (Penkov et al., 2013), suggesting that the alternative use of the limb enhancer that is repressed by Shox2 underlies the trunk to limb (proximal to distal) expression gradient of Meis. Notably, a large population of Shox2/Pbx co-occupied sites is not (or is only weakly) accessible to Pbx and Meis in the trunk region (Fig. 7L), suggesting that the Shox2/TALE co-occupied elements that are specifically accessible in the limb represent the epigenomic change required for the specification of the proximal limb mesenchyme. This highlights the importance of Shox2 in specification of the stylopod that emerges in juxtaposition of the trunk, with its strong Meis/Pbx expression and the proximal limb, where Shox2 is highly expressed. However, because of the closely located binding centers of Shox2 and Pbx genome wide (Fig. S7), and the lack of a Pbx de novo motif – a phenomenon that is typical for cooperative binding of transcription factors (Jolma et al., 2015) – Shox2-TALE co-occupancy appears to be competitive. Moreover, since Hox9/Hox10 proteins were thought to interact with Meis/Pbx to pattern the proximal limb (Capellini et al., 2011; Zakany and Duboule, 2007), the fact that Shox2 interacted directly with Hoxa9 in co-IP assay (Fig. S8) further suggests that on the Shox2/TALE co-occupied enhancers, Shox2 may also act as an alternative co-factor for Hox9/Hox10 to exert a different transcription output. Interestingly, in the palate that is derived from the Hox-free first BA (Amin et al., 2015), Shox2 displays a distinct pattern of genome-wide occupancy, implying that, in the absence of Hox factors, Shox2-TALE factors regulate a unique set of palate-specific modules.

Interestingly, a large population of Shox2/TALE co-occupied sites is also occupied by Gli3 (Vokes et al., 2008), an early anterior-posterior limb patterning factor and an effector for Shh signaling, in the developing limb bud (Fig. S9). In light of the fact that Hox and Shh show substantial epistatic cross talks and that Gli3 is essential for proximal limb skeleton development (Barna et al., 2005), the co-occupation of Gli3 and Shox2/TALE suggests that such cross talk is realized by converged transcription output on the same enhancers.

The fact that Shox2 mutation results in a specific ablation of the perichondrial Shox2⁺/Runx2⁺ osteogenic precursors and subsequently mispatterned cartilaginous mode in the stylopod highlight a crucial role for osteogenic precursors in chondrogenic differentiation. This conclusion is supported by phenotypic reminiscence in mice carrying Shox2 inactivation in the osteogenic lineage. Notably, the expression of Shox2 in many connective tissue cell types and its later expression in perichondrial outer layer cells are reminiscent of that of Hox11 (Swinehart et al., 2013), which patterns the zeugopodial skeleton elements. Together with the fact that the perichondrial expression of Hox13 was also proposed to pattern the autopodial skeleton (Villavicencio-Lorini et al., 2010), our study indicates that the patterning can be exerted by controlling the specification and differentiation of the perichondrial osteogenic precursor cells to simultaneously modulate the patterning of cartilaginous mode and bone formation. Interestingly, in the HoxA/D compound mutant, the expression of Shox2 in the embryonic perichondrium, but not chondrocytes, is severely affected (Neufeld et al., 2014), indicating that Shox2 itself is a key effector for Hox-TALE patterning system in the stylopod and further supports the idea that Hox-TALE patterns skeletal elements by commissioning perichondrial cells. Given the instructive role of perichondrial cells on chondrocyte differentiation and maturation (Ohbayashi et al., 2002), further dissection of Shox2-specific targets in the osteogenic lineage would further highlight key factors involved in such patterning function.

The Shox2^HA, Shox2^+/−, Shox2^LacZ, Shox2^Cre, Shox2^Flox, Rosa^mTmG, Nkx2-5^Flox, Col2a1-Cre and Col3.6-Cre alleles used in this study were described previously (Cobb et al., 2006; Liu et al., 2004; Ovchinnikov et al., 2000; Sun et al., 2013; Wang et al., 2014; Ye et al., 2015b; Yu et al., 2005). The animal experiments in this study were approved by The Tulane University Institutional Animal Care and Use Committee. Transgenic animal work performed at the Lawrence Berkeley National Laboratory (LBL) was reviewed and approved by the LBL Animal Welfare and Research Committee. Transgenic mouse assays were performed in Mus musculus FVB strain mice. Sample sizes were selected empirically based on our previous experience of performing transgenic mouse assays for >2000 total putative enhancers. Mouse embryos were only excluded from further analysis if they did not carry the reporter transgene or if they were not at the correct developmental stage. As all transgenic mice were treated with identical experimental conditions, and as there were no groups of animals directly compared in this section of the study, randomization and experimenter blinding were unnecessary and not performed.

Embryos from timed pregnant females were fixed in ice-cold 4% paraformaldehyde (PFA), embedded in paraffin, and sectioned at 8 µm for immunohistochemistry analyses. Antibodies used in the study are described in supplementary Materials and Methods. Alcian Blue/Alizarin Red skeletal staining was performed using established protocol, as described previously (Bobick and Cobb, 2012; Yu et al., 2007).

Immunoprecipitation (IP), co-immunoprecipitation (Co-IP), chromatin-immunoprecipitation (ChIP) and western blotting were performed as described previously (Ye et al., 2015b). ChIP-Seq and associated analysis are described in Supplemental Information. Antibodies used for ChIP and ChIP-Seq are described in supplementary Materials and Methods.

The GFP⁺ proximal embryonic limbs were dissected out from E12.5 Shox2^Cre/+;Rosa^mTmG and Shox2^Cre/−;Rosa^mTmG embryos under a fluorescence dissecting scope and subjected to digestion by a cocktail of collagenase I, II and IV, followed by a brief trypsin treatment. Suspended cells were sorted based on GFP signal and subsequently subjected to RNA extraction and RNA sequencing. The associated analysis is described in supplementary Materials and Methods.

All experiments were repeated at least three times. Quantification results are presented as mean±s.d., and statistical analysis was conducted using Student's t-test. For qPCR, reactions for each sample were also performed in triplicate. P<0.05 was considered significant.

The LHB-A-Luc plasmid used in this study was kindly provided by Jessica Rosin currently at Seattle Children's hospital.