Forelimbs (FLs) and hindlimbs (HLs) develop complex musculoskeletal structures that rely on the deployment of a conserved developmental program. Pitx1, a transcription factor gene with expression restricted to HL and absent from FL, plays an important role in generating HL features. The genomic mechanisms by which Pitx1 effects HL identity remain poorly understood. Here, we use expression profiling and analysis of direct Pitx1 targets to characterize the HL- and FL-restricted genetic programs in mouse and situate the Pitx1-dependent gene network within the context of limb-specific gene regulation. We show that Pitx1 is a crucial component of a narrow network of HL-restricted regulators, acting on a developmental program that is shared between FL and HL. Pitx1 targets sites that are in a similar chromatin state in FL and HL and controls expression of patterning genes as well as the chondrogenic program, consistent with impaired chondrogenesis in Pitx1−/− HL. These findings support a model in which multifactorial actions of a limited number of HL regulators redirect the generic limb development program in order to generate the unique structural features of the limb.
INTRODUCTION
In tetrapods, forelimbs (FLs) and hindlimbs (HLs) develop complex musculoskeletal structures within the framework of a common three-segment organization. Despite differences in structure and utility, the FLs and HLs of mice and humans both contain a proximal one-bone stylopod, an intermediate zeugopod segment with two bones, and a distal autopod with five digits. The arrangement of cartilage, bone, muscle and tendon within this framework, however, depends on the reproducible implementation of an underlying developmental program, and the differences in pattern between limbs arise from limb-specific modifications of this program, itself a derivation of the fin program of the lobe-finned fish ancestral to all tetrapods (Petit et al., 2017; Pieretti et al., 2015; Shubin et al., 1997).
Limb-specific transcription factors (TFs) are important elements of the limb program: paired-like homeodomain 1 (Pitx1) is an HL-restricted TF gene that is expressed throughout the posterior mesoderm and consequently in the early HL bud but not in FL (Lanctôt et al., 1997; Lamonerie et al., 1996). Genetics experiments in mice suggest that Pitx1 is an important upstream regulator of HL patterning: Pitx1−/− mice develop HLs that lack several key HL characteristics (Lanctôt et al., 1999; Szeto et al., 1999). This role is conserved in evolution (Chan et al., 2010). Strikingly, Pitx1−/− HLs fail to develop the load-bearing architecture of the knee: they lack a patella and have elbow-like bone contacts between the two bones of the zeugopod and the femur, which is twisted and shortened. At the molecular level, Pitx1−/− HLs are deficient in the expression of another limb type-restricted TF gene: Tbx4 (Lanctôt et al., 1999). Rescuing Tbx4 expression via a Prx1-driven Tbx4 transgene rescues some features lost in the Pitx1−/− HL, including femur length and certain qualitative muscle patterning characteristics (Ouimette et al., 2010).
In chick, leg and wing patterning has been correlated by various means with the expression of Tbx4 and the closely related FL-restricted Tbx5 (Gibson-Brown et al., 1998; Ohuchi et al., 1998; Rodriguez-Esteban et al., 1999; Takeuchi et al., 1999), but the capacity of these proteins to sufficiently determine HL versus FL morphology of the skeleton in the mouse is, however, less definitive (Ouimette et al., 2010). Pitx1 plays a prominent role in patterning the skeleton and has the capacity to generate HL-like features in FL: ectopic expression of Pitx1 in the FL of mice leads to structural changes in the FL skeleton, as well as ectopic expression of the HL-restricted TF Hoxc10 (DeLaurier et al., 2006; Minguillon et al., 2005). In humans, mutations leading to ectopic PITX1 expression in the FL cause Liebenberg syndrome, a disease in which the FLs of patients display HL-like patterning features (Spielmann et al., 2012).
The breadth and characteristics of the Pitx1-directed HL program are not well understood: what is the extent of Pitx1-dependent gene regulation and to what extent is this Pitx1-dependent gene regulatory network unique to the HL? Here, we use expression profiling of morphologically stage-matched FL and HL, combined with profiles of chromatin marks in these limbs, to illustrate the gene regulatory networks that drive FL and HL development. We also isolate the Pitx1-directed elements of the HL program through expression profiling of Pitx1−/− mice and ChIPseq of Pitx1 in E11.5 HL. We show that the programs that drive FL and HL development are very similar, marked by the expression of relatively few limb type-restricted genes. We also show that Pitx1 largely acts upon a chromatin landscape and genetic program that are common between FL and HL. Although the chromatin landscape of select Pitx1-dependent genes varies in accordance with their enriched expression in HL versus FL, the Pitx1−/− phenotype is principally defined by an expression loss of genes common to both limbs. In conjunction with the loss of anterior skeletal features in the Pitx1−/− HL, we observe direct Pitx1 targeting of chondrogenic genes in the HL, suggesting that Pitx1 generates HL-specific features by targeting common limb elements involved in chondrogenic expansion during early HL development.
RESULTS
The HL and FL programs are defined by a limited number of limb type-restricted genes
In order to examine the extent of the limb-specific developmental programs, we compared the transcriptomes of morphologically stage-matched FL and HL by RNAseq analysis. Targeting an early, comparable stage of development between FL and HL, we chose to evaluate E10.5 FL and E11.0 HL, stages at which the limb bud has completely emerged from the flank but before the bud acquires any hint of a paddle-like morphology. At this stage, the anterior and posterior margins of the bud are parallel to each other and perpendicular to the flank, which corresponds to a somite count of 33-36 somites in FL and 38-42 somites in HL, stages that will henceforth be referred to as E10.5 and E11.0, respectively (Fig. 1A). Using a false discovery rate (FDR) cutoff of 0.01, we observe 986 genes that are differentially expressed between FL and HL, although the vast majority of these differences are low-magnitude expression changes, i.e. less than 2-fold changes (Fig. 1B, Fig. S1). These low-magnitude changes, however, are not biologically spurious: for example, HLs show a 27% reduction in Gli3 expression and a corresponding 66% increase in Shh expression, and subtle changes in gene dosage of these anterior-posterior patterning genes have been previously associated with limb-specific morphological changes (Li et al., 2014).
We observe a relatively narrow network of genes that are expressed in a limb type-restricted manner, defined here as a log2 fold change greater than 3 (Fig. 1B,C). The extent of the FL-restricted gene network is effectively limited to Tbx5 and the predicted gene Gm43050, an antisense transcript present at the Tbx5 locus (Fig. 1C). In HL, the network of restricted genes is slightly more extensive, consisting of Pitx1, Tbx4, Isl1, and several of the 5′ HoxC genes. Much is known about these prominent genes of the limited HL-specific network: Pitx1 and Isl1 are genetically complementary, as both are upstream of Tbx4, although Pitx1 contributes to the development of anterior HL structures while Isl1−/− mice fail to develop posterior HL elements, such as the ischium and zeugopod (Itou et al., 2012). Hoxc10 and Hoxc11, in conjunction with their paralogs in the HoxA and HoxD clusters, are necessary in HL for the development of the stylopod and zeugopod, respectively (Wellik and Capecchi, 2003). Gm53, a long non-coding RNA located 5′ of Hoxb9, shows HL-restricted expression, although Hoxb9 itself, which was previously shown to be expressed preferentially in the leg of chick (Nelson et al., 1996), plays a role with the other Hox9 paralogs in the establishment of the posterior limb field in FL exclusively (Xu and Wellik, 2011). Rxfp1, a relaxin family peptide receptor gene that shows HL-restricted expression, is not known to be associated with limb patterning defects, although its ligand, relaxin, is associated with negative regulation of collagen turnover, TGFβ signaling, and fibrosis (Samuel et al., 2005). The presence of these known limb regulators in the differentially expressed dataset, although accompanied by a few novel limb type-restricted genes, validates the results of our transcriptomic approach.
The Pitx1 gene regulatory network partially overlaps with the network of limb-enriched gene expression
With a thorough characterization of FL versus HL gene regulatory networks in hand, revealing extensive low-magnitude expression changes and very few limb type-restricted genes, we sought to characterize the Pitx1 gene regulatory network at E11.0 to better understand its role in the HL program. We performed expression profiling of E11.0 Pitx1−/− versus wild-type (wt) HL and cross-compared this with our stage-matched FL versus HL expression profiling analysis (Fig. 2A). Using a conservative threshold (FDR<0.001) we identified 689 genes that are Pitx1 dependent. Of these, most are not expressed in a limb-specific manner: 538 genes show no preference between FL and HL at FDR<0.01 (Fig. 2B). The genes that do show limb-enriched expression can be separated into four categories: Pitx1-dependent genes with HL-enriched expression; Pitx1-dependent genes with FL-enriched expression; HL-enriched genes that show increased expression in the Pitx1−/− HL; and FL-enriched genes that show increased expression in the Pitx1−/− HL. Of these groups, the set of Pitx1-dependent HL-enriched genes and the set of FL-enriched genes that show an increase in expression in the Pitx1−/− HL define the limb-specific components of the Pitx1 gene regulatory network.
The set of Pitx1-dependent HL-enriched genes is notably marked by positive controls Pitx1 and Tbx4, but in addition to these expected genes we find Hoxa13 and several other 5′ genes of the HoxA cluster, including Evx1, Evx1os and Hottip (Fig. 2B, Table S1). Also prominently featured in this group is E130114P18Rik, which encodes a long non-coding RNA of unknown function; this gene shows a 4-fold preference in expression in HL, and a 4-fold decrease in HL expression in the absence of Pitx1. The Pitx1-dependent HL-enriched genes constitute the largest category with a limb-enriched expression profile, which also includes Tgfb3, Gtf2ird1 and, interestingly, Hand2, which is important for the establishment of limb bud polarity and the onset of Shh expression (Galli et al., 2010).
FL-enriched genes that show an increase in expression in Pitx1−/− HL represent potential components of the limb program that are normally downregulated in the endogenous HL developmental program, where Pitx1 is normally expressed. Pitx2 is notably present on this list, as are several genes that are typically expressed in anterior compartments of the wt FL, including Epha3 and Gria2 (Fig. 2B, Table S2A). The set of genes that show HL-enriched expression as well as increased expression in the absence of Pitx1 constitutes a small list that most prominently includes Isl1, a known genetic parallel of Pitx1, as well as the HL-specific Rxfp1 (Table S2B). It is likely that the genes in this list are components of the Isl1-dependent posterior HL developmental program that persists unimpeded by the loss of Pitx1.
The last set of genes, showing Pitx1-dependent and FL-enriched expression, are a peculiar group as they represent components that are more prominently featured in the FL program that are nonetheless significantly decreased in the Pitx1−/− HL (Table S3). This group contains genes known to be involved in diverse processes, from proximal limb development and chondrogenesis to, interestingly, many genes associated with the establishment of anterior domains in anterior-posterior limb patterning. Shox2, an important gene for proximal development of both HLs and FLs, is an upstream regulator of Runx2 (Cobb et al., 2006); both of these genes are FL enriched and Pitx1 dependent. Gli3 and Cdon, which have an antagonistic genetic relationship with Shh, are also present in this set, indicating Pitx1-dependent expression despite the relative suppression of these genes in wt HL versus FL (Mo et al., 1997; Probst et al., 2011; Cardozo et al., 2014). Taken together, these data highlight the complex nature of the Pitx1-dependent gene network, as many Pitx1-dependent genes in these four categories have different or even opposing functions.
Pitx1 targets a limb program that is already primed for use in FL
In order to investigate the direct targets of Pitx1, we performed ChIPseq of Pitx1 in E11.5 HL. We also performed ChIPseq of the FL-restricted TF Tbx5, with the aim of using the targets of these limb-specific TFs as seeds for the investigation of limb-specific enhancer elements. In total, we isolated 10,273 peaks in our Tbx5 dataset and over 50,000 peaks in our Pitx1 dataset. For Pitx1, however, we isolated the peaks that match a corresponding ChIPseq replicate, performed using a different antibody in the same tissue (Infante et al., 2013). This corresponding dataset shows fewer total peaks, but of the 25,025 present 71% are covered by our data. This final list of 17,783 peaks thus represents a highly reliable list of Pitx1 target sites (Fig. 3A). We find that Pitx1 and Tbx5 predominantly bind gene-distant, non-promoter regions (Fig. 3B). Interestingly, the top-scoring underlying motif for Tbx5 is not a T-box half-site (TCACACCT), but a composite T-box–Hox site that contains five bases of the T-box half-site directly abutting the six core bases of a Hox binding site (Fig. 3C). Pitx1 sites are principally enriched for the canonical TAATCC Pitx1 binding site, as well as Hox and bHLH sites; unlike Tbx5, however, these Hox sites are randomly dispersed around the principal Pitx1 motif and not in a sequence-specific orientation (data not shown).
In addition to the expression profiles of stage-matched FL and HL, we also generated a library of chromatin mark data in these limbs. Using the targets of Pitx1 and Tbx5 as seeds for a screen of possible limb-specific enhancer regions, we evaluated the profiles of chromatin modifications associated with active and repressed chromatin (Fig. 3D-G), namely acetylation of lysine 27 of histone H3 (H3K27ac) and trimethylation of lysine 27 of histone H3 (H3K27me3), respectively, in stage-matched FL and HL (Cao et al., 2002; Shlyueva et al., 2014; Cotney et al., 2012). For these analyses, only Pitx1 and Tbx5 peaks at least 2.5 kb from the nearest transcription start site (TSS) were used so as to avoid contaminating our analysis with promoter regions, which have a different chromatin organization compared with enhancers. These TSS-proximal peaks represent 16% and 13% of the total, respectively. We found that the chromatin surrounding direct Pitx1 and Tbx5 targets is in largely the same state in both FL and HL (Fig. 3D-G), comprising a bimodal H3K27ac peak and a depletion of H3K27me3, an indication of enhancer chromatin in an active state. The ratio of H3K27ac signal between FL and HL is similar at sites of Pitx1 or Tbx5 binding: in fact, between limbs and comparing 100 bp bins across the entire genome, all chromatin marks that we evaluated correlate with each other very strongly, with an average Pearson's correlation coefficient of 0.94 (Fig. S1B). In light of this similarity, it is clear that Pitx1 does not globally change the epigenetic state of the limb development program, but rather acts upon genomic targets of an epigenetic state shared between FL and HL.
A small set of Pitx1 targets and Hox loci show limb-enriched chromatin profiles
Although the overall pattern of Pitx1 targets suggests a common state between limbs, there are 129 non-promoter Pitx1 targets that show a 2-fold enrichment in H3K27ac signal in the 2 kb window spanning the Pitx1 peak summit. We isolated these sites, assigned them to their putative target genes using the GREAT tool (McLean et al., 2010), and then cross-referenced this list of putative targets with our set of HL-enriched genes at FDR<0.01 (Fig. 4A). The GREAT tool associates ChIPseq peaks with nearby genes within a 5 kb to 1 Mb dynamic genomic window. The sites most enriched for H3K27ac in HL are at the truly HL-restricted genes, namely Pitx1, Tbx4, Isl1 and the HoxC cluster, although several genes of the HoxA cluster also show HL-enriched H3K27ac at Pitx1 targets in conjunction with HL-enriched expression (Fig. 4B). This pattern corresponds with an enrichment of H3K4me3 over the 5′ side of the HoxA locus in HL, as well as a depletion of H3K27me3 over the same region (Fig. S2). This pattern of chromatin marks and HL-enriched gene expression is also apparent at the HoxD locus (Fig. S3). To varying extents, all Hox loci show a preference for HL expression in the 5′ side of the cluster and a corresponding preference for FL expression in the 3′ side of the cluster, although the boundaries at which expression becomes FL-enriched or HL-enriched vary by cluster (Figs S2-S5).
Of these HL-enriched genes near Pitx1 sites that show HL-enriched H3K27ac signal, only half are Pitx1 dependent. Pitx1 sites near Col5a1, Hoxc10 and Hand2 are shown as representative loci (Fig. 4C). Interestingly, in light of the fact that Pitx1 has been shown to be present at and drive Hoxc10 and Hoxc11 expression (Infante et al., 2013; Park et al., 2014; DeLaurier et al., 2006), Hoxc10 shows only a limited degree of Pitx1 dependence, with a 25% reduction in expression at E11.0 in Pitx1−/− HL (FDR=0.006), while Hoxc11 is not Pitx1 dependent at E11.0 and the expression of neither gene is Pitx1 dependent at E11.5 (Fig. 4B, Fig. S5A). These results reveal putative enhancers where Pitx1 action may increase active chromatin marks, but also highlight the complementary, interrelated nature of the HL-restricted gene network.
Pitx1 directly targets Sox9 and regulates the chondrogenic program in HL
In an effort to isolate the most certain and direct targets of Pitx1, we sought to find genes with Pitx1-dependent expression in HL at both E11.0 and E11.5 at a highly significant threshold, and then to focus on these differentially expressed genes that show evidence of Pitx1 binding within the associated locus. We used two strategies to correlate Pitx1 binding sites with target genes: published 4C contacts data, at loci for which it is available (Andrey et al., 2017), as well as the aforementioned GREAT tool. Using FDR<0.00001 in both the E11.0 and E11.5 Pitx1−/− expression profiles reveals 67 genes, the overwhelming majority of which show decreased expression in Pitx1−/−: 54 genes are downregulated, whereas only 13 genes show increased expression in the Pitx1−/− HL (Fig. 5). The top four Pitx1-dependent genes most frequently targeted by Pitx1 at regions of established enhancer-promoter contact are Sox9, Tcf7l2, Tbx15 and Tbx18 (Fig. S1C): Pitx1 is found at a total of 41 regions that contact these genes (Fig. 5).
To assess the dynamic gene expression changes over the course of development, we also isolated the set of genes that are differentially expressed in the Pitx1−/− HL at E11.5 that show no differential expression at E11.0 (Fig. 6A). We found that twice as many genes show decreased versus increased expression in the Pitx1−/− HL: 220 versus 102, respectively (Fig. 6A). Gene Ontology analysis shows that this larger set of Pitx1-dependent genes is enriched for genes associated with proteinaceous extracellular matrix, as well as collagens; in other words, genes associated with chondrogenesis and other tissues such as tendons, ligaments and muscle connective tissue (Fig. 6B).
We assayed the expression pattern of Sox9 in E11.5 wt and Pitx1−/− HL by in situ hybridization in order to assess the spatial changes in expression that correspond to the quantitative changes we observe in our E11.0 and E11.5 expression profiles. We see a prominent patch of Sox9 expression in the anterior proximal wt HL that is absent in the Pitx1−/− HL (Fig. 7A). We also tracked the skeletal phenotype of wt and Pitx1−/− HL, as well as wt FL, over the course of their development from E12.5 to E14.5, staining for cartilage with Alcian Blue 8GX (Fig. 7B). We see that this anterior patch of Sox9 expression corresponds to regions of the HL that develop into the robust elements of the distal femur and proximal tibia, corresponding to the regions most affected in the Pitx1−/− HL. This large anterior condensation never forms in the Pitx1−/−, and a weak pattern of condensation patterned in the basic Y-shape of the FL still lacks intensity relative to the FL skeletal phenotype (Fig. 7B). These condensations are thin and undeveloped at E13.5, and by E14.5 there is a visible kink in the developing femur of the Pitx1−/− HL, a contrast to the patch of cells that migrate down from the shoulder to form the deltoid tuberosity in the humerus (Blitz et al., 2013), and an indication that the early chondrogenic defects observed in the Pitx1−/− result ultimately in failed long-bone formation in the stylopod. Pitx1 direct targeting of Sox9 (Fig. 7C) and related chondrogenic genes at earlier stages of limb development might thus be responsible for the many skeletal defects observed in the Pitx1−/− HL at later stages.
DISCUSSION
The FL and HL of tetrapods evolved from the pectoral and pelvic fins, respectively, during the fin-to-limb transition 360 million years ago (reviewed by Clack, 2009). The role of Pitx1 in generating differences between posterior and anterior appendages has far-reaching evolutionary implications and is not unique to tetrapod limbs: the expression of Pitx1 has been associated with the presence or absence of pelvic fins in species of stickleback, which are ray-finned fish with distinct ancestry relative to the lobe-finned predecessors of modern tetrapods (Chan et al., 2010). Recent studies suggest that, between species, the early-stage limb development program is robustly conserved between divergent mammalian species (Sears et al., 2015). The classical model follows that FLs and HLs are serial homologs, although recently an alternative hypothesis has been put forth that the extensive similarities between FLs and HLs are the result of convergent evolution during the fin-to-limb transition. Regardless of the origin of the similarity between the developmental programs of FLs and HLs, however, it is clear that the responsibility of generating limb-specific features of the FL and HL programs rests with a narrow network of limb-specific elements of this program. Our results match these evolutionary constraints, with few limb type-restricted genes and a broadly conserved developmental program.
HL patterning is largely implemented by transcriptional, not epigenetic, actions
Pitx1, a prominent HL-restricted gene, can manifest its role in HL patterning by tinkering with the established and common limb program or by directly effecting skeletal pattern by regulation of the amount, timing and placement of chondrogenic condensation during limb development. It is clear, in light of our results and the abundance of data across species, that there is not an extensive, limb type-restricted network downstream of Pitx1. Our results indicate that Pitx1 binds to developmental enhancers that are ready to use in FL, which is to say that Pitx1 does not appear to engage in extensive remodeling of the chromatin state of its targets. This common state between FL and HL suggests that Pitx1 acts by direct transcriptional regulation rather than chromatin remodeling. The fact that the preponderance of genes misregulated in the Pitx1−/− HL are downregulated, as opposed to derepressed in its absence, combined with the active chromatin profile around Pitx1 target sites, suggests that Pitx1 predominantly functions to activate its transcriptional targets. Finally, the fact that Pitx1 targets a region with an HL-enriched chromatin profile near the Pitx1 gene itself might suggest autoregulation and, indeed, positive autoregulation of Pitx1 expression has been shown in the pituitary (Goodyer et al., 2003).
Transcriptional activation is usually accompanied by enhancement of active enhancer chromatin marks, but this is different from an action such as that of pioneer factors which instill chromatin accessibility where there was none. Our data do not support a significant pioneer role for Pitx1. Rather, the data indicate that Pitx1 acts as a classical TF at enhancers in an active chromatin state (Fig. 7D).
Pitx1 modifies the common limb program as a member of a narrow HL-specific network
Previous analyses suggest a limited overlap between the HL-specific components of the HL program: Isl1 partially drives Tbx4 expression early in limb development, but is not necessary for Tbx4 expression upon HL bud outgrowth, nor is it required for Pitx1 or Hoxc10 expression (Itou et al., 2012; Kawakami et al., 2011). Our analysis confirms the parallel and independent nature of Pitx1 and Isl1, as Isl1 expression is even increased in the absence of Pitx1. The prominent target of Isl1, however, which is known to be Hand2 (Itou et al., 2012), is shown here at E11.0 to be Pitx1-dependent. We do not see the complete abrogation of Shh expression in the Pitx1−/− HL, however, as is the case in Hand2−/− HL (Galli et al., 2010) and Isl1−/− HL (Itou et al., 2012), although there is a 40% reduction in Shh expression in the E11.0 Pitx1−/− HL (FDR=0.007). These results further complement transgenic experiments in which ectopic Pitx1 expression in the FL suppresses Shh signaling pathways and posteriorly expands the territory of Gli3 expression (DeLaurier et al., 2006).
At the same time, Pitx1 regulates components of the anterior patterning program, including Gli3, which are disrupted in the absence of Pitx1. Although much of the research on Shh and Gli3 focuses on digit patterning in the autopod, it was shown that Gli3−/− HLs display a reduction in femur size (Litingtung et al., 2002). The decrease in femur length in the Pitx1−/− might in part be due to the decreased expression of Gli3 in these HLs. In related work, Irx3 and Irx5 were shown to have HL-specific contributions to anterior limb pattern in a manner that is dependent on the interplay between Shh signaling and anterior patterning genes such as Gli3 (Li et al., 2014). Irx5 is not differentially expressed in Pitx1−/− HL at E11.0 or E11.5, however, and Irx3 shows only a slight decrease in expression in Pitx1−/− HL above the threshold of significance taken here (FDR=0.014). It is therefore likely that Pitx1 and Irx3/5 both contribute to the development of anterior structures in HL, but not interdependently.
Tbx15 and Tbx18 are revealed in this study as direct targets of Pitx1, with Pitx1-dependent expression in the HL, although both genes are also expressed in the FL. Tbx15 and Tbx18 are closely related members of the Tbx1 subfamily of T-box genes (Papaioannou, 2014). Tbx15 is strongly expressed in the core mesenchyme of the developing limb bud and Tbx15−/− mice show skeletal defects in both FL and HL, including reduced femur size (Singh et al., 2005). Tbx18 is expressed in the anterior-proximal limb mesenchyme (Kraus et al., 2001), although there are no morphological effects on the limb in Tbx18−/− mice (Bussen et al., 2004). In light of the fact that Pitx1 is necessary for the expression of Tbx4, Tbx15 and Tbx18, it is possible that there is a degree of functional redundancy between these genes, and that Pitx1 exerts an influence in the early HL mesenchyme by coordinating the expression of all three genes.
Our data also show that Pitx1 directly targets Sox9 (Fig. 7D), and there is a marked loss of Sox9 expression in the early Pitx1−/− HL in regions that correspond to the most prominent skeletal defects in the fully formed Pitx1−/− HL. Ultimately, bones do form in the Pitx1−/− HL, so the chondrogenic cascade is not completely disrupted and Pitx1 is not explicitly necessary for condensation and skeletogenesis in the HL. This does not preclude a role of Pitx1 expression in directly controlling the total number of cells that condense in the anterior HL. Pitx1 might alter the rate of proliferation of uncommitted mesenchymal progenitors, directly stimulate commitment to the Sox9-positive chondrogenic lineage, influence proliferation of Sox9-positive cells, or any combination of these possibilities. This direct manner of regulation might occur complementarily to, or in concert with, Pitx1-dependent regulation of T-box family genes and anterior-posterior patterning genes. Further, expression of these regulators in the early HL bud (Marcil et al., 2003) is consistent with the present interpretation.
How to make a different limb
In summary, the broad set of Pitx1 targets suggests a role as a major modulator of the limb development program. This action is implemented through transcriptional regulation (primarily activation) of a large set of target genes. These target genes are for the most part expressed in both FL and HL and therefore the effect of Pitx1 action (and of the downstream Tbx4, or of Tbx5 in FL) is to reorganize the gene expression patterns throughout the developing bud. It is this reorganized program – the sum of small changes that create a distinction in aggregate – that defines HL identity.
MATERIALS AND METHODS
Chromatin immunoprecipitation (ChIP)
All animal experimentation was approved by the IRCM Animal Ethics Review Board and followed Canadian guidelines. ChIP was performed using limb bud tissue collected from CD1 mice ordered from Charles River Laboratories or bred in the animal facilities of the IRCM. Embryonic day (E) 0.5 was designated to be noon of the day on which a plug was found. Pitx1 ChIPseq was performed with HL tissue from mice staged at E11.5, while Tbx5 ChIPseq was performed on the FL of mice staged at E10.5. Limb buds were collected in cold 1× PBS, cross-linked with 1% formaldehyde for 14 min, and chromatin was sheared using a manual sonication probe. Previously validated homemade rabbit antibodies against Pitx1 (Lamonerie et al., 1996; Tremblay et al., 1998; Marcil et al., 2003) and Tbx5 (Georges et al., 2008) were used, and all ChIPs were performed with equal amounts of Protein A and Protein G Dynabeads (Invitrogen). These TF ChIPseq samples were sequenced on an Illumina HiSeq 4000 platform, yielding 50 bp, paired-end reads, which were then mapped to mm10 using Bowtie2 version 2.2.6 (Langmead et al., 2009). Peaks were called using MACS2 with default parameters for paired-end reads (Zhang et al., 2008), and peaks with summits present in genome-wide RepeatMasker or an Encode-sourced comprehensive empirical blacklist of regions with artificially high signal were removed, as were summits whose loci contained anomalous input signal (RepeatMasker Open-3.0, http://www.repeatmasker.org) (ENCODE Project Consortium, 2012). Pitx1 summits from our data that appeared at any region within a peak defined in Infante et al. (2013) were considered to be present in both datasets. De novo motif searches were performed on peaklists using HOMER with default parameters (Heinz et al., 2010). Non-promoter sites are defined as summits greater than 2.5 kb from the nearest TSS present in the gene annotation from Encode version M11 (Ensembl 86). ChIPseq tracks for visualization were created in HOMER and displayed using Integrated Genome Viewer (IGV) (Thorvaldsdottir et al., 2013).
For the chromatin profiles of stage-matched FL and HL, samples were sequenced on the Illumina HiSeq 4000 platform, yielding 50 bp, single-end reads. These single-end read datasets were mapped to mm10 using Bowtie version 1.1.2 (Langmead et al., 2009). Genome-wide correlations were performed using EASeq (Lerdrup et al., 2016). Heatmaps, computed over a 5 kb window divided into 500 10-bp bins, as well as corresponding fillplots, were also created in EAseq.
Expression profiling by RNAseq
For stage-matched E10.5 FL versus E11.0 HL comparisons, limb tissue was collected from wt CD1 mice. Somites were counted for each embryo: 33-36 somite FLs were taken as E10.5 FLs and 38-42 somite HLs were taken as E11.0 HLs. Both E11.0 and E11.5 Pitx1−/− versus wt HL comparisons were performed in a Balb/c background. All animal experimentation was approved by the IRCM Animal Ethics Review Board and followed Canadian guidelines.
RNA samples were collected via silica gel spin-column purification. For E10.5 FL versus E11.0 HL, one replicate was prepared using a strand-specific ribosomal RNA-depleted RNAseq kit from Illumina (TruSeq Stranded Total RNA). One replicate of E11.0 Pitx1−/− versus wt HL was prepared using only the left HL, using a strand-specific ribosomal RNA-depleted KAPA kit (KAPA Stranded RNA-Seq Library) with Ribozero Gold (Epicentre), and one replicate of E11.5 Pitx1−/− versus wt HL was prepared using an unstranded, TruSeq mRNA enrichment library prep kit (Illumina). Two additional replicates of E10.5 FL versus E11.0 HL and E11.0 Pitx1−/− versus wt HL were prepared using the strand-specific ribosomal RNA-depleted KAPA kit with Ribozero Gold (Epicentre), and two additional replicates of E11.5 Pitx1−/− versus wt HL were prepared using an mRNA-enrichment kit from NEB [NEBnext poly(A) mRNA magnetic isolation] and the KAPA Stranded RNA-Seq Library.
For each comparison, three total replicates were fused for an experimental design with two conditions and two batches for each condition. All reads were mapped to the genome using STAR (Dobin et al., 2013). Transcript features were counted using featureCounts (Liao et al., 2014). All RNAseq experiments were analyzed using edgeR with robust dispersion estimates, likelihood ratio testing, and a generalized linear model (GLM), allowing us to accommodate batch effects between replicates (Robinson et al., 2010; McCarthy et al., 2012; Zhou et al., 2014).
In situ hybridization and cartilage stains
In situ hybridization was performed using anti-digoxigenin and probes generated from Sox9 cDNA plasmid digested with BamHI. Cartilage stains of E12.5 through E14.5 mouse embryos were performed by fixing embryos in Bouin's fluid, bleaching embryos in 70% ethanol/0.1% ammonium hydroxide, staining with 0.05% Alcian Blue 8GX in 5% acetic acid, and clearing embryos with 2:1 benzyl benzoate:benzyl alcohol solution, similar to Nagy et al. (2009).
Acknowledgements
We thank Isabelle Brisson, Dimitar Dimitrov and Sara Demontigny for their work in the IRCM animal care facility; Amandine Bemmo for helpful assistance with bioinformatics analyses; Mona Nemer, University of Ottawa, for the gift of Tbx5 antibody; and Évelyne Joyal and Tabasum Abdul-Rasul for secretarial work.
Author contributions
Conceptualization: S.N., J.D.; Methodology: S.N., J.D.; Formal analysis: S.N.; Investigation: S.N., M.L., A.H.S.; Resources: T.P., J.D.; Data curation: S.N., D.J.; Writing - original draft: S.N.; Writing - review & editing: S.N., J.D.; Visualization: S.N.; Supervision: J.D.; Project administration: J.D.; Funding acquisition: T.P., J.D.
Funding
This work was supported by grants (MOP-123213, CEEHRC EP1-120608) from the Canadian Institutes of Health Research (CIHR).
Data availability
ChIPseq and RNAseq data have been deposited at Gene Expression Omnibus under accession number GSE100734.
References
Competing interests
The authors declare no competing or financial interests.