Developmental gene expression patterns are orchestrated by thousands of distant-acting transcriptional enhancers. However, identifying enhancers essential for the expression of their target genes has proven challenging. Maps of long-range regulatory interactions may provide the means to identify enhancers crucial for developmental gene expression. To investigate this hypothesis, we used circular chromosome conformation capture coupled with interaction maps in the mouse limb to characterize the regulatory topology of Pitx1, which is essential for hindlimb development. We identified a robust hindlimb-specific interaction between Pitx1 and a putative hindlimb-specific enhancer. To interrogate the role of this interaction in Pitx1 regulation, we used genome editing to delete this enhancer in mouse. Although deletion of the enhancer completely disrupts the interaction, Pitx1 expression in the hindlimb is only mildly affected, without any detectable compensatory interactions between the Pitx1 promoter and potentially redundant enhancers. Pitx1 enhancer null mice did not exhibit any of the characteristic morphological defects of the Pitx1−/− mutant. Our results suggest that robust, tissue-specific physical interactions at essential developmental genes have limited predictive value for identifying enhancer mutations with strong loss-of-function phenotypes.
Embryonic development depends on the spatial, temporal, and quantitative control of gene expression. The complex regulatory programs that orchestrate morphogenesis are specified by thousands of tissue-specific, distant-acting enhancers, which are brought into contact with their gene targets via long-range chromatin looping interactions (Tolhuis et al., 2002; Pennacchio et al., 2006; Visel et al., 2009; Fullwood et al., 2009; Noonan and McCallion, 2010; Montavon et al., 2011; Sanyal et al., 2012; de Laat and Duboule, 2013; DeMare et al., 2013; Shlyueva et al., 2014; Ghavi-Helm et al., 2014; Rao et al., 2014). These contacts occur within stable, megabase-scale topologically associated domains (TADs), many of which are conserved across tissues and species (Dixon et al., 2012; Nora et al., 2012; Phillips-Cremins et al., 2013; Lupiáñez et al., 2015). Despite their clear overall importance for development, quantifying the contribution of individual enhancers to the expression of their target genes has been difficult. Genetic knockout studies of single enhancers in the mouse have yielded a range of molecular and phenotypic effects. In a classic example, deletion of the ZRS, a long-range enhancer that controls expression of Shh specifically in the developing limb bud, results in loss of Shh expression and truncation of the mouse limb (Lettice et al., 2003; Sagai et al., 2005). However, knockouts of other enhancers near developmental genes show less severe phenotypes (Attanasio et al., 2013; Dickel et al., 2016, 2018) or no obvious phenotype at all (Ahituv et al., 2007).
One potential mechanism to account for this result is that many developmental genes are regulated by multiple, redundant enhancers with overlapping functions that may compensate for the loss of a single enhancer (Hong et al., 2008; Degenhardt et al., 2010; Frankel et al., 2010; Perry et al., 2010; Montavon et al., 2011; Cannavò et al., 2016; Osterwalder et al., 2018). This redundancy potentially serves to buffer the effects of genetic disruptions of individual enhancers. Clearly, however, the loss of some enhancers cannot be buffered; the ZRS appears to be the primary source of regulatory information for specifying Shh expression in the limb bud (Lettice et al., 2003; Sagai et al., 2005). Identifying such crucial enhancers that are likely to exhibit overt loss-of-function phenotypes when disrupted has been challenging. Evolutionary conservation, tissue-specific chromatin signatures, and activity in transgenic reporter assays have all been used to map developmental enhancers in the genome (Noonan and McCallion, 2010). However, these metrics have not provided the means to distinguish enhancer mutations with large effects from mutations that can be buffered by compensatory mechanisms. This presents a major barrier for efforts to understand the contribution of regulatory variation to human disease.
We hypothesized that topological interactions would be a useful filter, in that individual enhancers crucial to tissue-specific expression of pleiotropic developmental genes would show robust, tissue-specific interactions with those targets. In this way, topology may be used to identify not just an enhancer's targets, but also to predict the potential effects of an enhancer deletion. Supporting this hypothesis, it has been shown that the ZRS interacts with Shh in the limb bud (Amano et al., 2009; Williamson et al., 2016) and that HoxD gene expression in the developing limb is specified by an array of interacting enhancers (Montavon et al., 2011). To identify potentially essential developmental enhancers, we leveraged the results of a recent study that utilized chromatin interaction analysis by paired-end tag sequencing (ChIA-PET) targeting the SMC1A cohesin subunit to identify over 2500 cohesin-bound chromatin loops in E11.5 mouse limb (DeMare et al., 2013). From this study, we identified a robust, hindlimb-specific interaction between the pleiotropic developmental gene Pitx1 and a sequence located 132 kb away, which we termed the Pitx1 distal enhancer (PDE).
Pitx1 encodes a homeodomain transcription factor prominently expressed in the hindlimb and mandibular arch (Szeto et al., 1999; Lanctôt et al., 1999). Pitx1+/− mice are overtly normal, exhibiting only minor morphological defects at low levels of penetrance (Alvarado et al., 2011). Pitx1−/− mice show early postnatal lethality, as well as loss of the ilium, a reduction in length of the femur, a loss of the patella, and a reduced, malformed mandible (Szeto et al., 1999; Lanctôt et al., 1999). Pitx1 is hypothesized to establish hindlimb identity, and overexpression of Pitx1 in the forelimb produces a hindlimb-like morphology (Logan and Tabin, 1999). Prior studies have identified Pitx1 enhancer mutations that result in substantial morphological phenotypes in non-mammalian systems (Chan et al., 2010; Domyan et al., 2016). We predicted that loss of the PDE may also have substantial molecular and morphological effects. The PDE displays strong, hindlimb-specific H3K27ac marking, a signature of active enhancers; in the forelimb the sequence is marked by H3K27me3, which is associated with PRC2-mediated repression (Cotney et al., 2012). The PDE is also highly conserved across amniote species, suggesting that it serves an ancient role in Pitx1 regulation.
We used CRISPR/Cas9 to delete the PDE in the mouse, thus abolishing any regulatory input or topological insulation conferred by the interaction. We anticipated two potential outcomes from the loss of the PDE: the recapitulation of morphological defects observed in Pitx1 null mice, or the formation of aberrant and/or compensatory interactions between Pitx1 and other potential regulatory elements in the absence of the primary interaction. Instead, we found PDE loss to have only mild effects on Pitx1 expression in the hindlimb and no overt effects on hindlimb or mandibular morphology. No changes in enhancer activity or new long-range interactions were identified, at least at a level detectable by our methods, that might have compensated for the loss of the PDE.
Mapping Pitx1 regulatory interactions
We first characterized the topological landscape of the Pitx1 promoter in the mouse embryonic forelimb and hindlimb using circularized chromosome conformation capture (4C) (Simonis et al., 2006). We were able to recapitulate the previously described Pitx1-PDE interaction in the wild-type mouse E11.5 hindlimb. This interaction is the most prominent topological interaction we identified involving the Pitx1 promoter in the genome (FourCSeq replicate 1 P=0.00754, replicate 2 P=0.00374; see Materials and Methods; Fig. 1, Table S1). Pitx1 also interacts with the PDE at a low level in the E11.5 forelimb (FourCSeq replicate 1 P=0.03067, replicate 2 P=0.05542). These interaction frequencies correlate with the expression levels of Pitx1 and the epigenetic profile of the enhancer in the hindlimb and forelimb. In the hindlimb, Pitx1 is expressed and the PDE exhibits high levels of H3K27ac (Fig. 1). Conversely, Pitx1 is not expressed in the forelimb, and both the Pitx1 promoter and the PDE sequence are marked by H3K27me3, indicating they are repressed in this tissue (Cotney et al., 2012).
To further characterize the regulatory landscape of the Pitx1 locus, we turned to a recent study that incorporated multiple functional and comparative genomic datasets to predict enhancers likely to be active in the limb (Monti et al. (2017). This study identified 31 predicted enhancers within the TAD that encompasses Pitx1, including the PDE (see Fig. 4 and Fig. S1). To determine whether these additional enhancers interact with the Pitx1 promoter, we analyzed Capture-C interaction maps obtained from mouse embryonic limb (Andrey et al., 2017). This study identified interactions between the Pitx1 promoter and the PDE in the hindlimb at multiple developmental time points (Fig. 1). This study also identified Pitx1 promoter interactions with four additional loci (Fig. S1), which we did not detect using 4C. One of these interactions involved the promoter of a nearby gene, H2afy. Three others involved putative limb enhancers, all of which showed H3K27ac marking in E11.5 forelimb and hindlimb, including a sequence shown to drive reporter gene expression in embryonic mouse forelimb and hindlimb at this time point (Fig. S2) (Pennacchio et al., 2006). However, the PDE shows significantly greater marking in hindlimb compared with forelimb (Fig. S1) (Cotney et al., 2012), and its strong, hindlimb-specific interaction with Pitx1 has been detected by multiple independent methods, suggesting that it might substantially contribute to the hindlimb-specific expression pattern of Pitx1.
Generation and validation of PDE−/− mice
In order to disrupt the Pitx1-PDE interaction, we used CRISPR/Cas9 genome editing to generate a Cre-loxP conditional deletion allele in C57BL6/J mice (Figs S2 and S3). The locations of the loxP insertions were designed to ensure that all local SMC1A binding sites identified would be removed (DeMare et al., 2013). We first validated the deletion of the PDE sequence via PCR and Sanger sequencing (Fig. S1). We then validated the ablation of the topological interaction using 4C (Fig. 1). The loss of the PDE sequence is also evident as a lack of sequencing reads in the deleted interval in limb H3K27ac ChIP-seq and 4C datasets (Fig. 1). We did not detect any evidence of increased interaction frequency in the genomic regions flanking the deletion, supporting that removal of the PDE completely disrupts the promoter-enhancer interaction.
Although Pitx1−/− mice exhibit neonatal lethality, we did not observe any reduction in viability in constitutive PDE−/− mice. PDE+/−×PDE+/− crosses yielded wild-type, PDE+/− and PDE−/− offspring at the expected 1:2:1 ratio (Table S2). We therefore compared wild-type and constitutive PDE homozygous knockout mice in all subsequent analyses.
Comparison of Pitx1 expression in wild-type and PDE−/− mice
We next considered the effects of the PDE deletion on Pitx1 expression. We examined two tissues that are known to both express Pitx1 during embryonic development and exhibit malformations in Pitx1−/− mice: the hindlimb and the mandibular arch. To analyze Pitx1 expression, we initially compared hindlimbs and mandibular arches from litter-matched wild-type and PDE−/− mice using RT-qPCR (Fig. 2, Fig. S4). We carried out a timecourse analysis of Pitx1 expression in the hindlimb from E10.5 through E13.5. Pitx1 expression was reduced by 13-16% in PDE−/− hindlimb at all four time points. This reduction was a significant effect of the PDE deletion, and not solely attributable to litter effects or variation across developmental time points (three-way mixed effect ANCOVA F1,20=18.624, P=3.36×10−4; Fig. 2; see Materials and Methods). We also carried out pairwise comparisons at each time point, which indicated that Pitx1 expression was significantly reduced in PDE−/− hindlimb at E12.5 and E13.5 (Mann–Whitney U-test P=0.012 and P=0.021, respectively; Fig. S4) (Mann and Whitney, 1947). Loss of Pitx1 expression was more substantial in PDE−/− E11.5 mandibular arch than in hindlimb (39% expression reduction, Mann–Whitney U-test P=8.6×10−4).
We then performed RNA-seq in wild-type and PDE−/− hindlimb and mandibular arch at E11.5 to identify potential global changes in gene expression due to reduced Pitx1 levels. In PDE−/− mice, Pitx1 showed a similar reduction in expression in both tissues, as we observed by RT-qPCR (14.5% in hindlimb, 32.14% in mandibular arch; Tables S3, S4). However, the change of expression in each tissue was not statistically significant after multiple testing correction [hindlimb DESeq Benjamini–Hochberg corrected (BH) P=0.9997, mandibular arch DESeq BH P=0.5946 (Anders and Huber, 2010)]. Expression of Tbx4, a known downstream regulatory target of Pitx1 in the hindlimb (Logan and Tabin, 1999), was reduced by 4.68% (P=0.295, BH P=0.9997). Our DESeq analysis did not identify any genes to be significantly altered in expression between wild type and PDE−/− after multiple testing correction in either tissue (at a threshold of BH P<0.05).
To gain insight into potential subtle effects of Pitx1 reduction of expression, we conducted gene ontology (GO) enrichment analyses of all genes exhibiting expression changes with a non-adjusted P-value <0.05 (271 hindlimb genes, 234 mandibular arch genes; Table S5). Genes that met this relaxed threshold in either tissue were associated with GO categories related to embryonic development and gene regulation. However, no ontologies reached a BH adjusted P-value <1.5×10−7 (Huang et al., 2008). This further suggests that the reduced dosage of Pitx1 did not cause substantial changes to regulatory networks or developmental processes in either the limb or the mandibular arch.
To determine if loss of the Pitx1-PDE interaction resulted in changes in the spatial distribution of Pitx1 expression in the hindlimb and mandibular arch, we visualized Pitx1 expression using whole-mount in situ hybridization on litter-matched wild-type and PDE−/− E11.5 embryos. In agreement with previous studies (Szeto et al., 1999; Lanctôt et al., 1999), our in situ analysis detected Pitx1 expression localized to the mandibular arch and hindlimb in wild-type embryos (Fig. 3). We observed the same pattern of expression in PDE−/− embryos, indicating that the spatial distribution of Pitx1 expression at E11.5 is not substantially altered by loss of the PDE (Fig. 3, Fig. S5).
Activity of the PDE in a transgenic enhancer assay
To evaluate the tissue-specific enhancer activity of the PDE independent of its interaction with the Pitx1 promoter or its genomic context, we used a transgenic mouse reporter assay (Pennacchio et al., 2006). We tested an 8.3 kb interval within the PDE that included the region of H3K27ac enrichment in hindlimb and all PhastCons conserved elements within the H3K27ac peak (Siepel et al., 2005). At E11.5, this sequence drove strong, reproducible lacZ reporter expression in the mandibular arch (four of seven positive embryos; Fig. S6). However, only three embryos showed reporter gene expression in the limb mesenchyme, and in inconsistent patterns, indicating that the PDE does not drive robust gene expression in the limb when removed from its endogenous location. These findings, together with the greater reduction of Pitx1 expression in mandibular arch versus hindlimb in PDE−/− embryos, support the conclusion that the PDE only contributes moderately to Pitx1 expression in the hindlimb despite its strong interaction with the Pitx1 promoter.
Comparing epigenetic signatures of regulatory element activity between wild-type and PDE−/− mice
Deletion of the PDE sequence removed the major interacting partner of the Pitx1 promoter we could detect in the hindlimb. We hypothesized that this loss would result in an unoccupied chromatin anchor point at the Pitx1 promoter, leaving it potentially accessible to other regulatory interactions that might compensate for loss of the primary enhancer. Such compensatory events might be detected as changes in the observed interaction frequency and epigenetic profile at other loci within the TAD containing Pitx1 (Dixon et al., 2012). We detected three 4C interactions within the Pitx1 TAD that reached reproducible, nominal significance in PDE−/− hindlimbs but not in wild type (Fig. 4, Fig. S1, Table S1). However, there is little evidence that these provide novel or compensatory regulatory information to Pitx1. Two of these regions did not involve potential enhancers based on H3K27ac enrichment, while one involved an interaction detected by Capture-C in wild-type hindlimb (Fig. S1) (Andrey et al., 2017). Our results therefore do not support widespread or robust changes in regulatory interactions arising due to loss of the PDE (Fig. 4). As a caveat, we note that the limited sensitivity of 4C does not allow us to rule out the formation of transient chromatin interactions or interactions restricted to a subset of cells in the limb bud.
To identify potential compensatory increases in enhancer activity (independent of topology and the limitations of 4C), we used ChIP-seq to profile H3K27ac in E11.5 hindlimbs derived from wild-type and PDE−/− mice. Compensatory enhancer activation events may appear as new H3K27ac peaks, or as increases in H3K27ac at active enhancers. We also profiled H3K27me3 in the E11.5 PDE−/− and wild-type hindlimb to identify any potential changes to repressed sequences in the PDE−/− background. H3K27ac signatures were very similar in the wild type and PDE−/−, both within the Pitx1 TAD and across the genome (Fig. 4). Genome-wide, we identified 11 H3K27ac and 33 H3K27me3 regions showing statistically significant differential marking between the wild-type and PDE−/− samples (DESeq BH P<0.01; Tables S6, S7). The Pitx1 gene body itself exhibited a significant overall loss of acetylation of ∼10% in PDE−/− mice (DESeq BH P=4.33×10−5), consistent with the loss of expression that we observed. Two other significant changes in H3K27ac signal were found on chromosome 13, although neither was located within the Pitx1 TAD. These comprised a 15% gain and a 10% loss 4.5 Mb and 25 Mb away, respectively, from Pitx1. Neither of these sites showed evidence of significant interactions with the Pitx1 promoter in our 4C data in either wild-type or PDE−/− hindlimb. The other interacting enhancers identified by Capture-C did not show significant changes in H3K27ac levels. Thus, we found no clear evidence of regulatory compensation at Pitx1 in E11.5 PDE−/− hindlimb.
Morphological analysis of PDE−/− mice
We also evaluated the impact of the Pitx1 PDE deletion on skeletal morphology. Previous studies have described the complete loss of the ilium and patella in Pitx1−/− mice (Szeto et al., 1999; Lanctôt et al., 1999). These mice also exhibit a shortened hindlimb stylopod and a malformed mandible. We examined skeletal morphology in 20 PDE−/− mice at E18.5 using Alizarin Red/Alcian Blue staining, and none showed any of the gross anatomical defects observed in Pitx1−/− mice (Fig. 5, Figs S7 and S8). Both the overall morphology and length of the hindlimb stylopod and mandible appeared unchanged in PDE−/− mice.
To detect potentially subtle reductions of the hindlimb we measured the length of the stylopod and zeugopod in litter-matched E18.5 pups. In our comparisons, we considered the ratio of the length of the stylopod and zeugopod (S/Z) to normalize for any natural variation in overall embryo size. We measured each limb segment both by the length of ossified bone alone as well as the combined length of the bone and associated cartilage. Across multiple litter comparisons, the average reduction of the S/Z ratio was under 2%. This reduction was not statistically significant using either metric (ossified bone two-way mixed effect ANOVA F1,18=2.32, P=0.145; combined bone and cartilage ANOVA F1,18=0.032, P=0.859; Table 1).
In this study, we disrupted a major topological interaction partner of the Pitx1 promoter in the mouse hindlimb. Based on the robustness and tissue specificity of the interaction, we predicted that the PDE would provide substantial quantitative, spatial or temporal regulatory input to Pitx1, and that PDE−/− mice might exhibit profound molecular and morphological phenotypes due to loss of Pitx1 expression in the hindlimb. However, our findings show that the deletion of the PDE and the resulting loss of the interaction does not result in reduced viability or recapitulation of the morphological phenotypes of Pitx1−/− mice (Szeto et al., 1999; Lanctôt et al., 1999). Although we did not observe major spatial or temporal changes in Pitx1 expression in PDE−/− mice, our findings do support that the PDE provides quantitative regulatory input to Pitx1. The loss of the PDE results in reduced Pitx1 expression in both the developing hindlimb and mandibular arch, as well as reduced H3K27ac signal at the Pitx1 promoter. These findings support that the PDE does enhance Pitx1 expression in both tissues, although the more modest reduction of Pitx1 expression in hindlimb due to PDE loss and our transgenic reporter assay results indicate that the PDE is more active in the mandibular arch than in the hindlimb. However, loss of the PDE is not sufficient to substantially alter hindlimb or jaw morphology in our model, indicating that the development of these structures is robust to moderate reduction in Pitx1 dosage. It is possible that PDE−/− mice exhibit hindlimb or jaw defects at very low penetrance, or subtle phenotypes that might be revealed by detailed morphometric analyses that are beyond the scope of our study (Attanasio et al., 2013).
Compensatory regulatory interactions that buffer Pitx1 expression against the loss of the PDE might account for the moderate molecular effects of disrupting the PDE-Pitx1 interaction. These compensatory interactions might involve recruitment of novel enhancers in the Pitx1 locus or upregulation of enhancers already in use, such as the other interacting enhancers detected by Capture-C (Andrey et al., 2017). We detected no evidence of robust novel interactions or enhancer upregulation in hindlimbs of PDE−/− mice. It is possible that transient or weak compensatory interactions might be occurring that our 4C methods are not sensitive enough to detect. However, the absence of a novel robust interaction similar to the PDE-Pitx1 interaction in PDE−/− mice, coupled with the lack of quantitative changes in H3K27ac marking in the Pitx1 TAD, suggest that Pitx1 regulatory architecture is largely unchanged by disruption of the PDE. Moreover, any compensatory mechanisms that might exist are necessarily incomplete, as deletion of the PDE has a clear negative effect on Pitx1 expression levels.
The PDE might thus provide only a moderate contribution to Pitx1 regulation, despite its prominent interaction with the Pitx1 promoter. Other regulatory elements in the Pitx1 locus, which might be partially redundant with the PDE, must determine the spatial, temporal and quantitative expression of Pitx1. Regulatory redundancy is an established feature of many developmental genes (Hong et al., 2008; Degenhardt et al., 2010; Frankel et al., 2010; Perry et al., 2010; Montavon et al., 2011; Shlyueva et al., 2014; Cannavò et al., 2016). One recent study supports that combinatorial disruption of multiple redundant enhancers is required to yield molecular or developmental phenotypes for such genes, although this study did not leverage interaction data to select candidate enhancers for deletion as we do here (Osterwalder et al., 2018). Our 4C analysis did not reveal other robust, long-range interactions that might be redundant with the Pitx1-PDE interaction in the hindlimb. However, the interacting enhancers previously detected by Capture-C might nevertheless contribute to Pitx1 regulation. These enhancers are marked by H3K27ac in both forelimb and hindlimb, but additional hindlimb-specific mechanisms might restrict their impact on Pitx1 expression to that tissue. Notably, a 44 kb deletion associated with feathered feet in pigeons and reduced Pitx1 expression in the hindlimb spans the pigeon ortholog of one of these enhancers (Fig. S1) (Domyan et al., 2016). H3K27ac profiles in the hindlimb also suggest that additional enhancers might be located near the Pitx1 promoter itself. We detected high levels of H3K27ac marking across the Pitx1 gene body, including intronic marking that might point to proximal enhancers (Fig. 1). Interactions with such proximal elements cannot be distinguished from background in 4C, given their location near the Pitx1 promoter viewpoint. Other distal enhancers within the Pitx1 TAD might also act redundantly with the PDE (Fig. 4) by interacting with the Pitx1 promoter at levels below the sensitivity of chromosome conformation capture methods. These enhancers might be active at earlier time points than those interrogated here or in other studies, and disrupting those elements might have greater effects on Pitx1 expression.
In light of our results, we speculate that the PDE may nevertheless be crucial for Pitx1 expression and normal development in ways that cannot be measured in a laboratory setting. The PDE is highly conserved, suggesting it serves an essential function. Redundant enhancers are thought to be maintained by selection in part because they confer a degree of regulatory robustness that buffers developmental gene expression against perturbation (Hong et al., 2008; Perry et al., 2010; Cannavò et al., 2016; Osterwalder et al., 2018). Embryonic development in wild populations is subject to environmental pressures and insults not present in experimental systems. The primary regulatory contribution of the PDE could be to stabilize Pitx1 expression against variation during development. In the absence of the PDE, Pitx1 expression might exhibit higher levels of sensitivity to environmental factors or genetic variation, and potentially show much greater changes than we can detect in our model. One approach to address this question in a future study would be to cross PDE null mice with Pitx1 heterozygous knockout mice, and determine if the resulting compound heterozygotes exhibit more severe phenotypes due to the greater reduction in Pitx1 dosage (Osterwalder et al., 2018).
Our results have implications for using topology maps to identify crucial regulatory interactions for developmental genes. A primary motivation for generating large-scale maps of topological interactions in the human genome is to discover noncoding variants that contribute to disease by disrupting long-range regulation and perturbing gene expression (Sanyal et al., 2012). However, as the PDE illustrates, robust, tissue-specific interactions may not serve to predict enhancer mutations that have strong effects in model systems. Although it has become commonplace to globally map putative regulatory interactions in multiple biological contexts, we still lack the means to distinguish essential versus redundant interactions a priori. This will require large-scale genetic disruption of many interactions, individually and in combination, in order to empirically measure the distribution of effects and characterize the functional diversity of regulatory interactions in the genome.
MATERIALS AND METHODS
Generating PDE conditional knockout alleles using genome editing
All animal work was performed in accordance with approved Yale IACUC protocols. Genetically modified mice were generated at the Yale Genome Editing Center (Yang et al., 2013). sgRNAs were selected for minimal off-target effects based on a CRISPR Design Tool score of >70 (http://crispr.mit.edu/) and the absence of target sites with three or fewer mismatches on the same chromosome (chr13) as the targeted sequence. C57BL/6J zygotes were injected with embryo microinjection buffer (10 mM Tris-HCl pH 7.4, 0.25 mM EDTA) containing two sgRNAs each at 50 ng/μl, Cas9 RNA at 100 ng/μl and two corresponding single-stranded DNA donor sequences both at 100 ng/μl filtered at 0.22 μm before injection (all sequences listed in Table S8). The Cas9 and sgRNA in vitro transcription templates were produced via PCR using a px330 plasmid into which the sgRNAs had been cloned as described (Cong et al., 2013) as a template, placing the sequences under control of the T7 promoter. The resulting DNA products were subjected to in vitro transcription using the MEGAshortscript T7 Kit (Thermo Fisher Scientific #AM1354) for the sgRNAs and mMESSAGE mMACHINE T7 ULTRA Transcription Kit (Thermo Fisher Scientific #AM1345) for Cas9. The transcribed RNAs were then purified with the MEGAclear Transcription Clean-Up Kit (Thermo Fisher Scientific #AM1908). Primers are described in Table S8.
Resulting F0 mice were backcrossed to wild-type C57BL/6J. To delete the PDE sequence, these mice were then crossed with an actin-Cre C57BL/6J mouse line provided by the Yale Genome Editing Center. All mice used in our analysis were from the F3 generation or later. The PDE knockout line has been deposited as cryopreserved sperm at the Mutant Mouse Resource and Research Center (MMRRC) (ID 43534).
4C analyses in hindlimbs from wild-type and PDE−/− mice
Hindlimbs and forelimbs (100 of each) were dissected from E11.5 C57BL/6J embryos. The tissue was crosslinked and processed in accordance with a protocol adapted from Naumova et al. (2012) and van de Werken et al. (2012) (associated 4C primers are listed in Table S8 and the 4C protocol is provided in the supplementary Materials and Methods). The resulting libraries were prepped and sequenced (1×75 bp) on an Illumina HiSeq 2500. For sequencing, the samples were indexed and multiplexed over two lanes. 30% PhiX control DNA (Illumina) was introduced into each sequencing run to mitigate Illumina sequencing artifacts arising due to the low complexity of PCR-amplified 4C libraries.
After sequencing, the 5′ and 3′ primer sequences were removed from the raw reads using cutadapt (Martin, 2011) (v1.4.1) (cutadapt --discard-untrimmed -g $firstprimer -n 10 -m 10 -O 10) and the trimmed reads were then aligned to the mm9 reference genome using Burrows-Wheeler Aligner (bwa) (Li and Durbin, 2009) (v0.7.10) (bwa mem -t 4). These aligned reads (ranging from 7.5-11 million per replicate) were then used to build a statistical model as described (Klein et al., 2015) modified to incorporate aligned sequence from both ends of paired-end reads. Raw reads were converted to wiggle tracks using BEDTools (Quinlan and Hall, 2010) (v2.17.0) for visualization.
PDE transgenic enhancer assay
We amplified an 8388 bp region [chr13:56056917-56065304 (mm9)] within the PDE with a nested PCR strategy using the primers shown in Table S8. The product was cloned into an Hsp68-lacZ reporter vector using Gibson Assembly (NEB #E2611). Generation of transgenic mice and embryo staining were carried out as previously described (Pennacchio et al., 2006).
RT-qPCR analysis of Pitx1 expression
Hindlimb RNA was collected from five 1:1 litter-matched pairs of mice at four developmental time points: E10.5, E11.5, E12.5 and E13.5 (five litters from each time point, one wild type and one PDE−/− from each litter, 40 total mice considered). RNA was purified from hindlimb pairs from each individual using the miRNeasy Micro Kit (Qiagen #74004) for E10.5 samples and miRNeasy Mini Kit (Qiagen #74106) for all others. RNA quality was measured using an Agilent RNA 6000 Pico chip; all samples analyzed had RIN values ≥8. We used 100 ng total RNA to prepare cDNA from each sample using the Invitrogen Superscript III Reverse Transcription Kit (#18080-051). The resulting cDNA was used in qPCR using Power SYBR Green Mastermix (Thermo Fisher Scientific #4367659) (qPCR primers listed in Table S8). All samples were run in sets of triplicates reporting the Ct values of Pitx1 and Hprt1, which was used as an internal reference. Reductions in gene expression were calculated by comparing ΔΔCt values of PDE and wild-type samples.
Global transcriptome analyses in hindlimbs and mandibular arch from wild-type and PDE−/− mice
RNA was collected from hindlimb bud pairs of four litter-matched E11.5 mice (two wild type, two PDE−/−) generated from a PDE+/−×PDE+/− cross. The mandibular arch analysis RNA was derived from three 1:1 pairs of litter-matched E11.5 mice (three litters, one wild type and one PDE−/− from each). RNA was purified using the Qiagen miRNeasy Kit (#74106). Prior to sequencing, RNA quality was analyzed on an Agilent 2100 RNA Pico Chip; all samples submitted for sequencing had RIN values ≥8. RNA was library prepped and sequenced using standard Illumina protocols (2×75 bp) on an Illumina HiSeq 2500 (mandibular arch) and Illumina HiSeq 4000 (hindlimb). Raw reads were aligned to the mm9 transcriptome (GRCm38.p4) using TopHat2 (Kim et al., 2013) (v2.0.9) and Gencode vM.7. Statistical analysis of expression changes was conducted using DESeq2 (Anders and Huber, 2010) using default settings and a false discovery rate (FDR) of 0.1.
Analysis of Pitx1 expression by in situ hybridization
Global profiling of H3K27ac and H3K27me3 in hindlimbs from wild-type and PDE−/− mice
We performed ChIP-seq for each epigenetic mark using chromatin derived from 180 wild-type and 180 PDE−/− hindlimbs (30 per replicate, three replicates per mark) from more than ten litters. Hindlimb tissue was crosslinked, pooled, and sonicated. The pools of each genotype were then split into three ChIP replicates for each epigenetic mark and processed independently for the remainder of the protocol. ChIP-seq of H3K27ac (2 μg/ChIP, Active Motif #39155) and H3K27me3 (5 μg/ChIP, Active Motif #ab4729) was performed as described (Cotney et al., 2012; Cotney and Noonan, 2015). Samples were sequenced (2×75 bp) on an Illumina HiSeq 2500. To control for batch effects, all H3K27ac samples were multiplexed and sequenced on a single lane. The same method was used for H3K27me3 samples. Raw reads were aligned to the mm9 reference genome using Bowtie 2 (v2.2.3) and peaks were called using MACS (Zhang et al., 2008) (v1.4.2).
For quantitative analysis, H3K27ac or H3K27me3 peaks that were reproducible in all three replicates for each genotype were merged. Peaks were considered reproducible if they overlapped by at least 1 bp (identified with the default settings of the intersectBed command in BEDTools). The relative enrichment of H3K27ac and H3K27me3 in these peaks was quantitatively compared using DESeq2 and an FDR of 0.1.
Analyses of skeletal morphology in E18.5 wild-type and PDE−/− mice
E18.5 skeletons were stained with Alcian Blue and Alizarin Red as described (Nagy et al., 2003). Litter-matched mice from six litters were used to obtain 30 wild-type and 24 PDE−/− independent hindlimbs. Bone and cartilage of the stylopod and zeugopod portions of the hindlimb were imaged under a stereo microscope (Leica S6 D) and the segment length measured independently and blinded to genotype by two individuals using Adobe Photoshop CC (2017.0.0). Both left and right limbs were used for these studies.
ANCOVA and ANOVA analyses
To dissect the main effects of the factors contributing to the observed differences in RT-qPCR and limb growth measurements, we analyzed these datasets in comprehensive modeling frameworks. For the RT-qPCR dataset, we employed a three-way mixed effect ANCOVA with ΔΔCT measurements as the response variable, ‘genotype’ as a two-level fixed effect factor, ‘embryonic age’ (days) as a covariate, and ‘litter’ as a random effect factor that is nested within ‘embryonic age.’ For the limb measurements, we used a two-way mixed effect ANOVA with stylopod-to-zeugopod ratio (both accounting for bone only and for bone and cartilage together) as the response variable, ‘genotype’ as a two-level fixed effect factor, and ‘litter’ as a random effect factor. In this way, we were able to isolate the effects of ‘genotype,’ our factor of interest, from the potential effects of other contributing factors.
We thank Jacques Drouin for providing the Pitx1 in situ probe. We also thank the Yale Center for Research Computing for high-performance computing support, the Yale Center for Genome Analysis for sequencing services, and the staff of the Yale Genome Editing Center for generating PDE null mice.
Conceptualization: R.S., J.P.N.; Methodology: R.S., E.V.D., S.D.W., T.N.; Validation: R.S.; Formal analysis: R.S., S.U.; Investigation: R.S., A.A.K., D.E., T.N.; Writing - original draft: R.S., J.P.N.; Writing - review & editing: A.A.K., D.E., S.U., S.D.W., J.P.N.; Visualization: A.A.K.; Supervision: J.P.N.; Project administration: J.P.N.; Funding acquisition: J.P.N.
This work was supported by a grant from the National Institute of General Medical Sciences (GM094780) to J.P.N., a Ruth L. Kirschstein National Research Service Award Individual Postdoctoral Fellowship (F32 GM106628) from the National Institutes of Health to D.E., a fellowship from the Deutsche Forschungsgemeinschaft (UE 194/1-1) to S.U., and a National Science Foundation Graduate Research Fellowship (DGE-1122492) to A.A.K. R.S. was supported by National Institutes of Health training grants T32 GM007223 and T32 HD007149. Deposited in PMC for release after 12 months.
All associated raw sequence data have been deposited at the Sequence Read Archive (BioProject ID: PRJNA388175).
The authors declare no competing or financial interests.