Cranial neural crest cell (NCC)-derived chondrocyte precursors undergo a dynamic differentiation and maturation process to establish a scaffold for subsequent bone formation, alterations in which contribute to congenital birth defects. Here, we demonstrate that transcription factor and histone methyltransferase proteins Prdm3 and Prdm16 control the differentiation switch of cranial NCCs to craniofacial cartilage. Loss of either paralog results in hypoplastic and disorganized chondrocytes due to impaired cellular orientation and polarity. We show that these proteins regulate cartilage differentiation by controlling the timing of Wnt/β-catenin activity in strikingly different ways: Prdm3 represses whereas Prdm16 activates global gene expression, although both act by regulating Wnt enhanceosome activity and chromatin accessibility. Finally, we show that manipulating Wnt/β-catenin signaling pharmacologically or generating prdm3−/−;prdm16−/− double mutants rescues craniofacial cartilage defects. Our findings reveal upstream regulatory roles for Prdm3 and Prdm16 in cranial NCCs to control Wnt/β-catenin transcriptional activity during chondrocyte differentiation to ensure proper development of the craniofacial skeleton.
Cranial neural crest cells (NCCs) give rise to chondrocytes that will form cartilaginous structures to provide the foundation for the craniofacial skeleton. The growth of these cartilage scaffolds depends on chondrocyte organization and positional orientation, which facilitates proper transitions from proliferative states towards hypertrophy and the eventual recruitment of osteoblasts to deposit bony matrices (Kimmel et al., 1998, 2010; Keller et al., 2000; Hall, 1978a; Le Pabic et al., 2014). These processes require extensive temporal and spatial regulation as alterations to the gene regulatory networks (GRNs) and signaling modules that control craniofacial chondrocyte maturation can impact the development of these structures and contribute to craniofacial defects (Biosse Duplan et al., 2016; Chai et al., 2000; Keith, 1910; Lei et al., 2016; Manocha et al., 2019; Ramaesh and Bard, 2003; Ricks et al., 2002; Svandova et al., 2020; Wilson and Tucker, 2004; Hall, 1978b).
Canonical Wnt/β-catenin signaling plays a highly dynamic role during chondrocyte differentiation. Although required early in prechondrogenic cells, β-catenin activity is downregulated as chondrocytes begin differentiation, becoming active again at late stages of chondrocyte hypertrophy to promote osteoblast differentiation and mineralization (Ben-Ze'ev and Geiger, 1998; Day et al., 2005; Hartmann and Tabin, 2000; Hill et al., 2005; Hwang et al., 2005; Ryu et al., 2002; Willert and Nusse, 1998). Abnormally high or low levels of β-catenin at inappropriate stages of chondrogenesis are inhibitory to differentiation and can cause abnormal hypertrophy, chondrocyte disorganization, stimulation of osteoblast differentiation and premature mineralization (Day et al., 2005; Hartmann and Tabin, 2000; Hill et al., 2005; Hwang et al., 2005; Ryu et al., 2002; Willert and Nusse, 1998). Although the processes governing cranial NCC-derived chondrocyte differentiation during craniofacial development have been well characterized, the upstream mechanisms driving spatial and temporal activation and repression of GRNs and signaling pathways, in particular canonical Wnt/β-catenin, remain largely unknown.
Several chromatin modifiers have important regulatory roles in cranial NCC and craniofacial development (Liu and Xiao, 2011). Among these are the PRDM (positive regulatory domain) family of lysine methyltransferases, which control gene expression by epigenetic modulation of chromatin accessibility, directly binding DNA via zinc-finger domains, or interacting with other protein complexes (Di Zazzo et al., 2013; Fog et al., 2012; Hohenauer and Moore, 2012). Human genome-wide association studies have associated single nucleotide polymorphisms in the genes encoding two of these PRDMs, PRDM3 (EVI1/MECOM) and PRDM16, with craniofacial abnormalities, including cleft lip/palate and variation in facial morphology (Jugessur et al., 2010; Shaffer et al., 2016; Li et al., 2019; Liu et al., 2012; White et al., 2021).
PRDM3 and PRDM16 are important in a variety of developmental processes in mouse and zebrafish, acting as both direct and indirect repressors or activators of specific GRNs depending on the cellular context (Baizabal et al., 2018; Bjork et al., 2010; Ding et al., 2013; Harms et al., 2015; Kajimura et al., 2009, 2008; Seale et al., 2008, 2007; Warner et al., 2007; Arai et al., 2011; Goyama and Kurokawa, 2010; Goyama et al., 2008; Sato et al., 2008). Although some genetic compensation exists between these two paralogs with their high amino acid sequence homology and similar developmental expression patterns, there is also evidence suggesting that they have independent functions (Shull et al., 2020). Both Prdm3 and Prdm16 are necessary for craniofacial skeletal development (Ding et al., 2013; Shull et al., 2020), yet the exact mechanism(s) involved in mediating cranial NCC-derived cartilage differentiation and craniofacial development are unknown.
In this study, we utilized zebrafish and mouse genetic models to dissect the molecular functions of Prdm3 and Prdm16 in chondrogenesis during craniofacial development. Based on our data, we hypothesize that Prdm3 acts as a transcriptional repressor whereas Prdm16 serves as an activator of gene expression to balance canonical Wnt/β-catenin transcriptional activity during chondrocyte differentiation. In doing so, Prdm3 and Prdm16 control proper spatial and temporal development of the vertebrate craniofacial skeleton.
prdm3 and prdm16 are required for craniofacial chondrocyte stacking and polarity
Previous studies have shown that loss of either prdm3 or prdm16 causes moderate overall craniofacial phenotypes, namely mild hypoplasia, in the developing zebrafish pharyngeal skeleton (Shull et al., 2020). We confirmed these observations by assessing the cartilage and bone phenotypes of prdm3 and prdm16 mutant 6-day-old zebrafish larvae (Fig. 1A-C′). However, high magnification of the hypoplastic cartilage elements within the craniofacial skeleton revealed significant abnormal chondrocyte organization (Fig. 1D). Unlike the stacked chondrocytes of wild-type larvae, chondrocytes in prdm3−/− and prdm16−/− mutants were highly disorganized. This was quantified by measuring the angle between adjacent chondrocytes in the direction of growth of that cartilage element (Fig. 1E,F). The more organized the cells, the closer the angle between adjacent cells to 180°. Loss of prdm3 and prdm16 caused a significant (33.4% and 34.6%, respectively) reduction of this angle (Fig. 1F). The number of cells per unit of area within these cartilage structures was increased in both prdm3−/− (∼32.3%) and prdm16−/− (∼29.6%) mutants (Fig. 1G). Previous work has demonstrated that outgrowth of pre-cartilaginous condensations requires cell-cell intercalations and changes in cell shape (hypertrophy) rather than extensive proliferation (Kimmel et al., 2010, 1998; Le Pabic et al., 2014). Loss of prdm3 and prdm16 causes no change in proliferation during early development of cartilage structures (Shull et al., 2020). As such, the increased number of cells per tissue area observed in prdm3 and prdm16 mutants likely results from the failure of these chondrocytes to intercalate properly and their inability to expand their cell shape for growth. To evaluate changes in cell intercalation and extension, live imaging from 56 hours post-fertilization (hpf) to 72 hpf was performed on wild-type, prdm3−/− and prdm16−/− zebrafish larvae crossed into the Tg(−4.9sox10:EGFP) (Carney et al., 2006) transgenic background. Live imaging revealed failure of the chondrocytes within cartilage structures to extend during development in both prdm3 and prdm16 mutants (Movies 1-3, Fig. S1A-C) and suggested that chondrocyte cell polarity was also abrogated in prdm3 and prdm16 mutants.
To assess changes in cell polarity, wild-type, prdm3−/− or prdm16−/− larvae were stained with acetylated α-tubulin to label microtubule-organizing centers (MTOCs) at 75 hpf, a time point at which chondrocyte cell-cell rearrangements and polarity have stabilized, and cells are oriented in a specific direction for growth. MTOCs in pre-chondrocytes are oriented towards the center of each condensation of the cartilage element. Over time, these cells re-orient themselves in a uniform manner along dorsal-ventral axes (Le Pabic et al., 2014). At 75 hpf, the MTOCs in the chondrocytes of the palatoquadrate in wild-type larvae were stable and localized ventrally toward the Meckel's cartilage and the jaw joint (Fig. 1H,H′, Fig. S2A-A‴). However, the MTOCs in stage-matched prdm3−/− (Fig. 1I,I′, Fig. S2B-B‴) and prdm16−/− (Fig. 1J,J′, Fig. S2C-C‴) failed to rearrange and orient themselves uniformly along this ventral axis and instead were positioned dorsally, or towards the center of the original condensation (Fig. 1I-J′, Fig. S2B-C‴). Quantification of MTOC orientation in wild-type palatoquadrates showed that 65% of chondrocytes at this stage were ventrally polarized, as denoted by MTOC puncta positioned directionally at the 180° quadrant of the cell (Fig. 1K,L,O). Conversely, loss of prdm3 and prdm16 caused a significant reduction (∼25%) in the percentage of uniformly ventrally oriented chondrocytes (180°), which corresponded to an increase in the number of chondrocytes that were instead positioned dorsally (0°) or towards the center of the cartilage element (90° or 270°) (Fig. 1M-O). These results demonstrate Prdm3 and Prdm16 are important for facilitating proper chondrocyte differentiation, including orientation, intercalation and organization in the cartilage elements that form during craniofacial cartilage morphogenesis.
Neural crest-specific function of Prdm3 and Prdm16 is required for chondrocyte organization in Meckel's cartilage of the developing murine mandible
Prdm3 and Prdm16 are expressed in the facial prominences during murine development (Fig. S3A-B′) and we have recently shown that Prdm3 and Prdm16 share some conserved functions across vertebrate species (Shull et al., 2020). To determine whether changes in craniofacial chondrocyte development were conserved in mammals with loss of Prdm3 and Prdm16, we conditionally ablated Prdm3 and Prdm16 in the murine neural crest lineage using the Wnt1-Cre driver (Chai et al., 2000; Danielian et al., 1998). Both homozygous Prdm3fl/fl;Wnt1-Cre+/Tg and Prdm16fl/fl;Wnt1-Cre+/Tg mutant animals were born at Mendelian ratios and survived to embryonic day (E) 18.5. Viability of mutant embryos past E18.5 was not assessed, but we predict that neither mutant would survive postnatally owing to either defects in heart development (Prdm3) or failure to thrive due to cleft palate defects (Prdm16) (Goyama et al., 2008; Hoyt et al., 1997; Bjork et al., 2010; Shull et al., 2020; Warner et al., 2013). Homozygous Prdm3 mutant animals developed a subtle craniofacial phenotype, namely mild anterior mandibular hypoplasia and slight defects in snout extension (Fig. 2B-B″,E). Homozygous Prdm16fl/fl;Wnt1-Cre+/Tg mutants presented with a variety of craniofacial defects, similar to those observed with loss of Prdm16 in the early epiblast, suggesting a cell-autonomous function in the neural crest (Bjork et al., 2010; Shull et al., 2020). Abnormalities included snout extension defects, anterior mandibular hypoplasia, secondary cleft palate and middle ear defects with severe hypoplasia of the tympanic rings, incus, malleus and retroarticular process of the squamosal bone (Fig. 2C-C″,E, Fig. S3C-H″).
The shared phenotype in both mutants was mandibular hypoplasia (Fig. 2A-E). Alcian Blue- and Alizarin Red-stained skeletal preparations confirmed anterior mandibular hypoplasia (Fig. 2A′-C″), which was quantified by measuring the ratio of the anterior portion of the mandible relative to the posterior length (Fig. 2D,E). This phenotype parallels some of the hypoplasia seen in prdm3 and prdm16 zebrafish mutants and suggests changes to the cellular development and maturation of the chondrocytes in cartilaginous structures (Meckel's cartilage) that support the formation of the mandible. To assess chondrocytes within the Meckel's cartilage histologically, sagittal sections were collected from control, Prdm3fl/fl;Wnt1-Cre+/Tg and Prdm16fl/fl;Wnt1-Cre+/Tg at E14.5, the stage of cartilage development at which the chondrocytes have undergone rapid growth and are starting to undergo pre-hypertrophy (Fig. 2F-H″). Safranin O and Fast Green staining of sections revealed changes to the chondrocyte organization and morphology in mutant animals compared with controls. In controls, chondrocytes were starting to swell and become pre-hypertrophic as evidenced by expansion of the intracellular space (Fig. 2F′,F″). These cells were also organized into stacks extending from the top of the cartilage structure to the bottom of the Meckel's cartilage (Fig. 2F′,F″). By contrast, Prdm3 mutant chondrocytes were smaller and tightly packed, suggesting a stall in their maturation process (Fig. 2G′,G″). Prdm16 mutant chondrocytes were unsynchronized; some appeared compacted like those in the Prdm3 mutants whereas others seemed to be undergoing accelerated pre-hypertrophy (Fig. 2H′,H″). In both mutants, the chondrocytes were disorganized (Fig. 2G″,H″).
Quantification of chondrocyte cell area and the number of cells per tissue area supported these observations of cellular changes. Loss of Prdm3 led to a significant ∼22% decrease in chondrocyte cell area (Fig. 2I) and a corresponding ∼16% increase of total number of cells per tissue area (Fig. 2J). Decreased cell area was near significance (P=0.0658), but there were no significant changes in the total number of cells with loss of Prdm16, likely because of the heterogeneric differentiation state of chondrocytes observed in these animals (Fig. 2I,J). The cellular changes observed in the developing Meckel's cartilage with loss of Prdm3 and Prdm16 suggest altered chondrocyte differentiation, which likely contributes to the anterior mandibular hypoplasia observed in these mutants. Importantly, chondrocytes within both zebrafish and mouse craniofacial structures were abnormal with loss of Prdm3 and Prdm16, suggesting a conserved function of these factors in regulation of chondrocyte differentiation and maturation in vertebrate craniofacial cartilage development.
prdm3 and prdm16 regulate global gene expression in cranial neural crest cells
To dissect the molecular mechanisms of prdm3 and prdm16 in controlling cranial NCC cartilage derivative differentiation, we analyzed the transcriptome of zebrafish cranial NCCs. Each mutant line was crossed into the Tg(−4.9sox10:EGFP) background (Carney et al., 2006). sox10:EGFP-positive cranial NCCs were isolated from wild-type, prdm3−/− and prdm16−/− embryos at 48 hpf by fluorescence-activated cell sorting (FACS) and subjected to RNA sequencing (RNA-seq). Transcriptomic analysis revealed striking differences in overall gene expression with loss of prdm3 or prdm16 (Fig. 3). In prdm3−/− mutants, an overwhelming 2688 genes were significantly upregulated, but only 189 genes were significantly downregulated (Fig. 3A,B). Conversely, in prdm16 mutants the majority of differentially expressed genes were significantly downregulated (1370), and only 279 were significantly upregulated (Fig. 3A,D). Gene ontology (GO) pathway analysis was performed separately on the genes that were upregulated in prdm3 mutants and those downregulated in prdm16 mutants. Canonical Wnt/β-catenin was identified as the top signaling pathway enriched in both the upregulated genes in prdm3 mutants and the downregulated genes in prdm16 mutants (Fig. 3B-E, Fig. S4). The opposing differentially expressed genes belonging to the canonical Wnt/β-catenin signaling axis included ctnnb1 and other signal transducers (Fzd gene family, Dvl gene family, apc) as well as downstream Wnt/β-catenin transcriptional targets (fosab, jun, tcf7, lef1) (Fig. 3B,D, Fig. S4C,D). Many of these genes encoded factors involved in the assembly of the Wnt/β-catenin enhanceosome transcriptional complex [bcl9/bcl9l, pygo1/2, smarca4a (also known as brg1), ep300a/b, arid1aa/ab, tcf7, lef1] (Fig. 3B,D, Fig. S4C,D). These opposing transcriptomic profiles in prdm3−/− and prdm16−/− mutants were validated by real-time quantitative polymerase chain reaction (RT-qPCR) on RNA from whole heads of wild-type, prdm3−/− or prdm16−/− embryos at 48 hpf (Fig. 3F-H). The trends of gene expression of canonical Wnt pathway genes (ctnnb1, fzd3b, dvl3, apc) (Fig. 3F), components of the Wnt-enhanceosome transcriptional complex (bcl9, pygo1/2, smarca4a) (Fig. 3G) and downstream Wnt/β-catenin target genes (tcf7, lef1, jun, fosab) (Fig. 3H), followed the same pattern: elevated expression with loss of prdm3 and, conversely, decreased expression with loss of prdm16.
To determine whether changes to the expression of Wnt/β-catenin signaling components identified by RNA-seq in zebrafish were also differentially expressed during formation of the mouse Meckel's cartilage, RNA was extracted from the mesenchyme of the mandibular processes of pharyngeal arch 1 at E11.5 from control, Prdm3fl/fl;Wnt1-Cre+/Tg and Prdm16fl/fl;Wnt1-Cre+/Tg embryos and RT-qPCR was performed for Wnt/β-catenin signaling components (Fig. 3I-N). Consistent with the zebrafish differential expression, loss of Prdm3 in developing murine mandibular facial tissue led to an increase in Ctnnb1 (Fig. 3I) expression, as well as members of the Wnt enhanceosome (Bcl9, Pygo2, Brg1) (Fig. 3J) and subsequent downstream Wnt/β-catenin target genes (Tcf7, Jun, Fos) (Fig. 3K), whereas loss of Prdm16 caused decreased expression of these factors (Fig. 3L-N). Although there were some differences in expression between species (i.e. decreased Apc and Lef1 expression in Prdm3 mutant mouse tissues, in contrast to elevated expression in zebrafish), the overall trends were similar and suggest that Prdm3 and Prdm16 exert opposing effects to balance Wnt/β-catenin signaling component gene expression and that this process is conserved across tetrapods and teleosts.
prdm3 and prdm16 control β-catenin stabilization and localization in craniofacial chondrocytes
Given that canonical Wnt/β-catenin signaling was the most differentially altered pathway in both prdm3 and prdm16 mutant zebrafish, and many of the differentially expressed genes associated with Wnt/β-catenin signaling were those involved in the formation of the Wnt/β-catenin enhanceosome transcriptional complex and retention of nuclear β-catenin, we next assessed changes in active β-catenin localization. Whole-mount immunofluorescence for nuclear β-catenin (phosphorylated tyrosine residue 489) was performed on wild-type, prdm3−/− or prdm16−/− zebrafish larvae at 75 hpf. Unlike wild-type chondrocytes, which have low abundance of nuclear β-catenin at this stage (Fig. 4A,D, Fig. S5A-A‴), prdm3−/− had a significant increase in the presence of nuclear β-catenin (Fig. 4B,D, Fig. S5B-B‴). Conversely, prdm16−/− had a dramatic reduction in the accumulation of nuclear β-catenin to levels significantly below that of wild-type chondrocytes (Fig. 4C,D, Fig. S5C-C‴).
To understand further canonical Wnt/β-catenin signaling in chondrocyte development in the craniofacial skeleton with loss of prdm3 and prdm16, prdm3−/− and prdm16−/− zebrafish mutants were crossed into the Wnt reporter line Tg(7xTCF-Xla.Sia:NLS-mCherry)ia5Tg (Moro et al., 2012) in combination with the Tg(fli1:EGFP) line (Lawson and Weinstein, 2002) to visualize Wnt-responsive cells in the pharyngeal arches and developing craniofacial cartilage elements (Fig. 4E-G). At 72 hpf, prdm3 mutants had increased Wnt-responsive cells in the pharyngeal arch/craniofacial elements and other structures, including the heart, otic vesicle and perichondrium (Fig. 4F, Fig. S5E-E″). Conversely, few or no Wnt-responsive cells were present in prdm16−/− larvae (Fig. 4G, Fig. S5F-F″). Together, these results validate the transcriptional changes observed in both prdm3 and prdm16 mutants and correlate with differences in expression of factors (i.e. bcl9) involved in retaining β-catenin in the nucleus, assembling the Wnt enhanceosome and driving transcription of downstream target genes. Abnormal accumulation of nuclear β-catenin in prdm3 mutant chondrocytes at later stages of chondrocyte development may explain the chondrocyte polarity, orientation, intercalation and differentiation defects, as high sustained β-catenin signaling can be inhibitory to cartilage differentiation (Hill et al., 2005; Hwang et al., 2005; Ryu et al., 2002; Sun et al., 2020). However, prdm16 mutants have the same chondrocyte differentiation defects, but instead with reduced canonical Wnt/β-catenin signaling, suggesting that low maintenance levels of nuclear β-catenin activity may be required to help facilitate chondrocyte differentiation and maturation.
prdm3 and prdm16 alter chromatin accessibility at cis-regulatory regions of Wnt/β-catenin signaling components in cranial NCCs
Previous work has implicated both Prdm3 and Prdm16 in the regulation of gene expression through methylation of lysine residues on histone H3: H3K9, a mark of repression, and H3K4, a mark associated with gene activation. Both H3K9me3 and H3K4me3 are significantly reduced in both prdm3 and prdm16 mutant zebrafish (Shull et al., 2020), suggesting global changes to the chromatin landscape. To understand how Prdm3 and Prdm16 are epigenetically regulating gene expression, specifically of canonical Wnt/β-catenin pathway components, assays for transposase-accessible chromatin paired with sequencing (ATAC-seq) were performed on FACS-sorted sox10:EGFP cranial NCCs isolated from pooled prdm3 wild-type/heterozygotes (hets), prdm3−/−, prdm16 wild-type/hets, and prdm16−/− embryos at 48 hpf (Fig. 5).
prdm3 mutant embryos exhibited a noticeable increase in open chromatin compared with matched wild-type/het samples (Fig. 5A,B, Fig. S6A,B). To quantify differences in open chromatin state, we used DiffBind to measure read depth on a union peakset generated from peaks called by Genrich in the prdm3 mutant and wild-type/het samples (Fig. 5E-G). This analysis identified 7870 peaks in the prdm3−/− embryos compared with 146 differential peaks identified from wild-type/het control siblings (Fig. 5E); 4917 regions were identified in both the prdm3−/− and wild-type/het samples (Fig. 5E). Significantly different differential peaks exhibited greater read depth in the prdm3 mutant samples compared with their wild-type controls (Fig. 5F). This result is illustrated by the greater number of statistically significant differential ATAC peaks in prdm3−/− mutants (Fig. 5G; note that volcano plot is shifted with more peaks having greater depth in the prdm3 mutant). We interpret these data to reflect a greater amount of open chromatin in the prdm3 mutants at a genome-wide scale.
In contrast, loss of prdm16 resulted in a slight decrease in open chromatin compared with matched wild-type/het samples (Fig. 5C,D, Fig. S6C,D). DiffBind analyses for prdm16 mutants identified only 702 differential peaks in prdm16−/− embryos compared with 972 in the wild-type/het controls with an overlap of 4079 peaks across both genotypes (Fig. 5H). Although few peaks reached statistical significance in the Prdm16 comparison, we note that all differential peaks do exhibit a slight decrease in the average read depth in the prdm16 mutants (Fig. 5I,J). Thus, we found an overall trend whereby read depth at ATAC peaks was generally lower in the prdm16 mutants compared with wild-type/het control samples, even if most peaks did not reach statistical significance. Importantly, this pattern is opposite to that observed for prdm3 mutants, for which differential ATAC peaks often exhibited higher read depth in prdm3 mutants (compare Venn diagrams and volcano plots in Fig. 5E-J).
ATAC peak regions with greater depth in the prdm3−/− mutants were largely associated with promoters and exons (Fig. 5K), in contrast to matched prdm3 wild-type/het controls for which peak regions were instead largely located within intronic and intergenic regions (Fig. 5K). The relatively smaller number of ATAC peaks with greater depth in the prdm16−/− mutants were mostly associated with intergenic regions (Fig. 5L). Across both mutant comparisons, 380 shared peaks were identified with greater read depth across the two mutants compared with their respective controls. Thus, approximately half of the peaks with greater chromatin accessibility in the prdm16 mutants (702 total) also had greater chromatin accessibility in the prdm3 mutants.
We next used the footprinting pipeline TOBIAS (transcription factor occupancy prediction by investigation of ATAC-seq signal) (Bentsen et al., 2020) to predict transcription factor occupancy in regions exhibiting differential chromatin states in the prdm3 and prdm16 zebrafish mutants. For these analyses, we limited our peak sets to the set of unique peaks identified by Genrich for each genotype regardless of whether DiffBind identified read depth at a peak as statistically significant or not. Thus, for these analyses we used TOBIAS to search for enriched transcription factors predicted to be uniquely bound in either the prdm3 or prdm16 mutants relative to their wild-type/het controls.
For regions unique to prdm3−/− mutants, the top 10 transcription factor families with predicted increased occupancy based on motif enrichment were associated with genes necessary for cell growth, differentiation, and developmental patterning processes (Fig. 5M, Tables S1 and S2). Several of these transcription factors are activated downstream of Wnt signaling [Ovol1, Hes7 (Her1), Tead2 (Tead family), Gbx2] whereas others (Etv5, Arnt2) recruit protein complexes to modulate chromatin architecture or facilitate gene expression through enhancers (Fig. 5M). One of the top transcription factor families predicted to be differentially bound in prdm3−/− mutant cranial NCCs was Fos/Jun. Fos/Jun is a known transcriptional target of canonical Wnt/β-catenin signaling, and our RNA-seq identified Fos/Jun as differentially expressed at the transcriptional level in both prdm3 and prdm16 mutants (Fig. 3).
The top 10 transcription factors enriched in prdm16−/− cranial NCCs included several factors that are important in neural crest cell development (Tfap2b) and formation of complexes important for facilitating chromatin remodeling [Plagl1 (Plag family)]. In addition, one of the top factors included Pparg. Prdm16 is already known to interact with Pparg to regulate brown fat adipogenesis (Kajimura et al., 2009, 2008) (Fig. 5N, Table S3). Pparg and Wnt/β-catenin can antagonistically influence expression of each other: activated Pparg decreases Wnt/β-catenin and elevated Wnt/β-catenin inhibits Pparg (Davis and Zur Nieden, 2008; Moldes et al., 2003; Takada et al., 2010). This pattern of Wnt/β-catenin expression and occupancy of Pparg correlates with the decreased Wnt/β-catenin signature observed in prdm16 mutants. Because most genes were downregulated and associated with less-accessible chromatin in prdm16−/−, the transcription factors that had a predicted decrease in binding affinity in the mutants were also assessed (Fig. 5O, Table S4). One of the top transcription factors with significantly decreased occupancy was Barx2. Barx2 has important functions in regulating cell-cell adhesions and aggregations and plays a crucial role in craniofacial development and chondrogenesis (Jones et al., 1997; Meech et al., 2005).
Because Wnt/β-catenin signaling was the most differentially regulated pathway from our transcriptomic data, we next assessed changes in chromatin accessibility at differentially expressed Wnt/β-catenin signaling components and target genes (Fig. 5P-T, Fig. S6E-I). Consistent with our transcriptomic data, loss of prdm3 led to increased open chromatin at canonical Wnt/β-catenin components, including ctnnb1 (Fig. 5P), bcl9 (Fig. 5Q) and pygo2 (Fig. 5R), as well as the downstream Wnt/β-catenin targets jun (Fig. 5S) and fosab (Fig. 5T). Conversely, chromatin accessibility of these Wnt/β-catenin targets was dramatically decreased in prdm16−/− (Fig. 5P-T). Areas of open chromatin aligned with areas of active gene expression identified by a previously published H3K27Ac ChIP-seq study of whole zebrafish embryos at 48 hpf (Bogdanovic et al., 2012). To determine whether prdm3 and prdm16 can directly bind to promoter regions of these canonical Wnt/β-catenin signaling component genes, cleavage under targets and release using nuclease paired with quantitative real-time PCR (CUT&RUN-RT-qPCR) was performed on whole zebrafish larvae at 48 hpf. Significant enrichment of Prdm3 and Prdm16 abundance was identified at putative binding sites for Prdm3 and Prdm16 near promoter regions of ctnnb1 (Fig. 5U) and the downstream Wnt/β-catenin target gene jun (Fig. 5V). Only Prdm3 had significant enrichment at fosab (Fig. 5W), suggesting differences in target genes between Prdm3 and Prdm16. There was no enrichment of Prdm3 or Prdm16 at promoter regions of bcl9 or pygo2 (Fig. S7C,D). Although this result suggests that these factors may not be direct transcriptional targets of Prdm3 or Prdm16, this does not exclude a possible protein-protein interaction influencing Wnt/β-catenin enhanceosome activity indirectly. Taken together, these results suggest that Prdm3 and Prdm16 may not only influence transcription of Wnt/β-catenin signaling components by changing chromatin accessibility near promoter regions of these genes but may also affect expression of these Wnt/β-catenin genes by binding those targets directly. These results further emphasize opposing roles of Prdm3 and Prdm16 in facilitating a balance of canonical Wnt/β-catenin signaling during cranial NCC-derived chondrocyte differentiation in the craniofacial skeleton.
Pharmacological manipulation of canonical Wnt/β-catenin signaling in prdm3−/− and prdm16−/− zebrafish mutants partially restores chondrocyte stacking defects
Our genomic data suggest that Prdm3 and Prdm16 regulate Wnt activity in opposing ways to balance Wnt/β-catenin transcriptional activity during chondrocyte differentiation. We next sought to rescue cartilage phenotypes by pharmacological manipulation of Wnt signaling in both zebrafish mutants. To rescue the effects of elevated Wnt/β-catenin signaling in prdm3 mutants, wild-type and prdm3−/− embryos were treated with either DMSO or 0.75 μM of the Wnt antagonist IWR-1 from 24 hpf to 48 hpf. IWR-1 blocks Wnt-induced β-catenin accumulation by stabilizing the Axin2 destruction complex (Chen et al., 2009). At 6 days post-fertilization (dpf), control and treated embryos were collected and stained with Alcian Blue and Alizarin Red to assess cartilage and bone (Fig. 6A-A‴). The chondrocytes in vehicle-treated prdm3 mutants were highly disorganized (Fig. 6A′). IWR-1 treatment significantly restored chondrocyte organization within mutant cartilage structures (Fig. 6A‴), as quantified by measuring the angle between adjacent chondrocytes (Fig. 6B).
To mitigate the effects of decreased Wnt/β-catenin signaling in prdm16−/− larvae, wild-type and prdm16−/− embryos were treated with either DMSO or 0.05 µM of a Wnt/β-catenin agonist, GSK-3 inhibitor XV, from 24 hpf to 48 hpf. GSK Inhibitor XV blocks Gsk3-dependent phosphorylation of β-catenin, allowing for β-catenin stabilization. Whereas chondrocytes of prdm16−/− vehicle-treated larvae exhibited stacking defects (Fig. 6C,C′), Wnt/β-catenin activation completely rescued chondrocyte defects (Fig. 6C‴,D). GSK inhibition (canonical Wnt pathway activation) also completely rescued overall cartilage phenotypes in prdm16−/− mutants, including restoration of posterior ceratobranchial cartilages (Fig. S8A). Together, these results show that chondrocyte defects in prdm3 and prdm16 mutants can be restored by rebalancing Wnt/β-catenin signaling pharmacologically in each mutant.
Combined genetic loss of both prdm3 and prdm16 rescues chondrocyte stacking defects
Although our genomic data suggest that Prdm3 and Prdm16 regulate Wnt/β-catenin activity in opposing directions, we found that single loss of prdm3 or prdm16 causes similar cartilage phenotypes: defects in chondrocyte orientation, polarity, intercalation and growth. To determine whether genetic loss of both prdm3 and prdm16 would rescue or exacerbate the craniofacial cartilage phenotypes observed in the single mutants, double-homozygous prdm3−/−;prdm16−/− mutants (and all other allelic combinations) were generated from prdm3+/−;prdm16+/− heterozygous intercrosses (Fig. 6E-H). Surprisingly, loss of both prdm3−/− and prdm16−/− rescued the craniofacial phenotypes observed in single-mutant embryos with normal posterior arch cartilage structures, and partially restored the hypoplasia observed in other cartilage elements (Meckel's cartilages, trabeculae, ethmoid plate) (Fig. 6H,H′). prdm3−/−;prdm16−/− double mutants also had a near-complete rescue of chondrocyte stacking and intercalation defects (Fig. 6H″,I). Combinatorial mutants in the allelic series did not rescue overall cartilage phenotypes or cellular stacking defects, suggesting that the rescue is dependent on complete loss of both paralogs (Fig. 6I, Fig. S8B,C). Although the near-complete rescue of cartilage defects was surprising in the double mutants, we did not fully assess the larval viability of the prdm3−/−;prdm16−/− double mutants throughout later stages of larval development, juvenile growth and adulthood. Because Prdm3 and Prdm16 also function in other tissues, this rebalancing of gene expression would be required to also occur in these tissues in order for viability to be restored. Based on the latest stage of larval growth we did assess for cartilage phenotypes (6-8 dpf), prdm3−/−;prdm16−/− double mutants phenotypically look healthy and indistinguishable from wild-type larvae, suggesting a possible rescue of viability. Future experiments will test this idea further. Together, these results emphasize the independent functional roles of Prdm3 and Prdm16 to balance GRNs and signaling modules, in particular canonical Wnt/β-catenin signaling, during craniofacial cartilage development.
In this study, we have identified functional roles for Prdm3 and Prdm16 in controlling proper chondrocyte differentiation by balancing canonical Wnt/β-catenin activity, both transcriptionally and epigenetically (Fig. 7). These two proteins have long been thought to act redundantly of each other, given their high amino acid sequence homology and similar developmental expression patterns. Contrary to this paradigm, our data highlight that each paralog has its own divergent and independent role in facilitating proper chondrocyte differentiation. Interestingly, Prdm3 and Prdm16 result from a vertebrate-specific duplication that occurred in the Gnathostomata ancestor (Vervoort et al., 2016). As such, independent functions of Prdm3 and Prdm16, at least in the case of chondrocyte maturation and bone differentiation in craniofacial development, may have coincided with their duplication event, which in turn correlates with the evolution of jawed vertebrates. In support of this idea, we hypothesize that, under normal developmental conditions during cranial NCC chondrocyte differentiation, Prdm3 functions as a traditional transcriptional repressor whereas Prdm16 acts as a non-traditional activator of gene expression to balance Wnt/β-catenin activity. Both mutants, despite their opposing transcriptomic profiles and chromatin landscape, have impaired chondrocyte orientation, polarity, intercalation and growth. Intriguingly, combined loss of both prdm3 and prdm16 dramatically rescues, rather than exacerbates, chondrocyte phenotypes. Furthermore, these results support a Prdm3- and Prdm16-dependent Wnt/β-catenin ‘Goldilocks effect’ during craniofacial chondrogenesis: there must be a precise level of Wnt/β-catenin activity during chondrocyte differentiation as too much or too little activity can have consequential effects on the proper maturation and differentiation of these cells.
Despite sharing similar craniofacial cartilage phenotypes, loss of prdm3 or prdm16 surprisingly causes drastically different changes to the transcriptome and chromatin accessibility. It will be interesting to investigate further whether these differences are due to an evolutionarily derived sub-functionalization or specialization of these proteins. It is also possible that Prdm factors have developed a sequence-level difference in specificity for their intrinsic methyltransferase activity. In addition to their methyltransferase abilities, both Prdm3 and Prdm16 commonly associate with other protein complexes to facilitate transcriptional regulation. Interestingly, there is very little overlap between known binding partners across both proteins, except for a few (Smad3, CtBP) (Palmer et al., 2001; Kajimura et al., 2009, 2008; Warner et al., 2007; Kurokawa et al., 1998). As such, the context of the protein complexes with which each Prdm is associated likely dictates how each Prdm influences gene expression and/or chromatin remodeling. Given the transcriptomic and chromatin accessibility changes we observe in each mutant, we predict that, perhaps at this developmental time point, Prdm3 is interacting in repressive complexes whereas Prdm16 is associating more strongly with co-activating complexes. Teasing apart these diverging functional characteristics and the potential for transient changes to these roles over developmental time will facilitate a better understanding of the specific molecular mechanisms of these Prdm factors in craniofacial cartilage development.
Mechanisms of transcriptional adaptation have been used to explain genetic compensation between genetic alleles derived from CRISPR gene editing. We have observed some level of genetic compensation in both prdm3 and prdm16 mutants, whereby prdm16 is modestly elevated in prdm3 mutants (Shull et al., 2020). However, prdm3 is decreased in prdm16 mutants. Interestingly, expression of a third Prdm family member, prdm1a, is increased in both prdm3 and prdm16 single mutants and combined loss of all three alleles causes drastically more severe craniofacial phenotypes (Shull et al., 2020). As such, it is possible that genetic interactions or feedback regulatory loops between these alleles could also influence the phenotypes we observe in these single mutants.
We have identified a GRN centered on balancing temporal and spatial canonical Wnt/β-catenin signaling during cartilage development that is facilitated by Prdm3 and Prdm16. We hypothesize that there are several mechanisms, both direct and indirect, that may be driving this temporal regulation. Among the differentially expressed Wnt/β-catenin signaling components, the members of the Wnt enhanceosome transcriptional complex, in particular the intermediary protein that facilitates β-catenin nuclear localization and assembly of the Wnt enhanceosome complex, bcl9/bcl9l, and its binding partner, pygo1/2, were significantly altered in both prdm3 and prdm16 zebrafish mutants. Previous work has shown that bcl9 zebrafish mutants have cranial NCC defects, including craniofacial cartilage phenotypes (Cantu et al., 2018). It would be interesting to determine whether Prdm3 and Prdm16 participate in the Wnt enhanceosome transcriptional complex through protein interactions with members of the complex, including bcl9/bcl9l, pygo1/2, or others to influence Wnt/β-catenin transcriptional activity. It could be possible that Prdm3 and Prdm16 are downstream Wnt/β-catenin transcriptional targets that would form a regulatory circuit, which could facilitate the timely gradient of Wnt/β-catenin necessary during chondrocyte differentiation and maturation during craniofacial development.
Although we have focused only on canonical Wnt/β-catenin signaling, it is possible that non-canonical Wnt/planar cell polarity (PCP) signaling is also active. Recent studies have suggested that crosstalk between the two pathways coordinates developmental processes (Navajas Acedo et al., 2019). In alignment with these studies, we did observe changes in several Wnt/PCP genes (vangl2, scrib, fat) as well as the transcription factors Jun/Fos. Although jun and fos are direct canonical Wnt/β-catenin target genes, they are also activated in response to non-canonical Wnt/PCP through JNK signaling. Because both prdm3 and prdm16 mutants share cartilage phenotypes similar to those observed in PCP zebrafish mutants, we speculate that Prdm3 and Prdm16 could also control crosstalk between both canonical and non-Wnt/β-catenin signaling activity during chondrocyte differentiation.
Finally, we define conserved functions of Prdm3 and Prdm16 across vertebrates. We show that loss of Prdm3 and Prdm16 in the murine neural crest lineage disrupts Meckel's cartilage chondrocyte organization and maturation owing to dysregulated Wnt/β-catenin signaling during early chondrogenesis in the mandibular facial processes. We hypothesize that this abnormal chondrocyte maturation could impact the development of the mandible, which could explain the anterior mandibular hypoplasia seen in these animals. It will be interesting to understand further how Prdm3 and Prdm16 regulate this potential coupling mechanism between chondrocyte maturation and subsequent bone formation in the developing mandible.
In summary, we have defined the roles for the chromatin modifiers Prdm3 and Prdm16 in cranial neural crest development and formation of the craniofacial skeleton. We show that these two seemingly functional redundant paralogs can act antagonistically independent of each other upstream of canonical Wnt/β-catenin signaling during chondrocyte differentiation to ensure proper spatial and temporal development of the vertebrate craniofacial skeleton.
MATERIALS AND METHODS
Zebrafish were maintained as previously described (Westerfield, 2007) in accordance with standard zebrafish husbandry conditions. Embryos were raised in defined Embryo Medium at 28.5°C and staged developmentally following published standards as described (Kimmel et al., 1995). The wild-type strain used was the AB line (ZIRC). The transgenic lines used were Tg(−4.9sox10:EGFP) (Dutton et al., 2008), Tg(fli1:EGFP) (Lawson and Weinstein, 2002) and Tg(7xTCF-Xla.Sia:NLS-mCherry)ia5Tg (Moro et al., 2012). These transgenic lines were crossed to the various mutant backgrounds. Zebrafish mutant lines for prdm3 and prdm16 were generated by CRISPR-based mutagenesis as previously described (Shull et al., 2020). The prdm3 and prdm16 mutant alleles used in this study (prdm3CO1005 and prdm16CO1006) are predicted frameshift mutations that interrupt the coding sequence upstream of the PR/SET domain responsible for functional histone methyltransferase activity (Shull et al., 2020). For double mutants, prdm3+/− fish were bred to prdm16+/− fish to generate prdm3+/−;prdm16+/− double-heterozygous animals, which were intercrossed to generate prdm3−/−;prdm16−/− double mutants as well as all the other resulting various allelic combinatorial animals. All experiments were completed on zebrafish embryos or larvae before sex was determined in these animals (20 dpf). The Institutional Animal Care and Use Committee of the University of Colorado Anschutz Medical Campus approved all animal experiments performed in this study and the procedures conform to NIH regulatory standards of care.
Mecomtm1mik (referred to as Prdm3fl/fl) (Goyama et al., 2008), B6(SJL)-Prdm16tm1.1Snok/J (referred to as Prdm16fl/fl) (The Jackson Laboratory) and H2afvTg(Wnt-Cre)11Rth (referred to as Wnt1-Cre+/Tg) (Danielian et al., 1998) were all maintained on the C57/Bl6 background and housed at a sub-thermoneutral temperature (21-23°C) under a 12 h light/dark cycle with water and food (PicoLab Rodent Diet 20) provided ad libitum. For timed matings, Prdm3fl/fl or Prdm16fl/fl females were bred to Prdm3fl/+;Wnt1-Cre+/Tg or Prdm16fl/+;Wnt1-Cre+/Tg males, respectively. The morning a vaginal plug was detected was considered E0.5. Embryos of matching somite numbers were used for experiments. Mice were euthanized by carbon dioxide inhalation followed by cervical dislocation as a secondary method of euthanasia. Male and female embryos were analyzed in this study and there were no sex-dependent differences in phenotypes between the two groups. Developmental stages of embryos used are indicated in the results. All embryos were stage-matched by somite counting. Mice were bred and maintained in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The protocol was approved by the University of Colorado Anschutz Medical Campus's Institutional Animal Care and Use Committee.
Fin clips, single whole embryos or single embryo larva tails were lysed in lysis buffer [10 mM Tris-HCl (pH 8.0), 50 mM KCl, 0.3% Tween-20, 0.3% NP-40, 1 mM EDTA] for 10 min at 95°C, incubated with 50 ug of Proteinase K at 55°C for 2 h, followed by inactivation of Proteinase K at 95°C for 10 min. Genotyping for prdm3 and prdm16 mutant alleles was performed as previously described (Shull et al., 2020). See Table S6 for genotyping primers.
For mice, tail clips from weanlings and tail clips or yolk sacs from embryos were lysed in DNA lysis buffer [10 mM Tris-HCl (pH 8.0), 100 mM NaCl, 10 mM EDTA (pH 8.0), 0.5% SDS] and 100 µg of Proteinase K overnight at 55°C. Genomic DNA was isolated following phenol/chloroform extraction. DNA pellets were air-dried and re-suspended in nuclease-free water. Genotyping for Prdm3, Prdm16 and Wnt1-Cre alleles was performed as previously described (Shull et al., 2020). See Table S6 for genotyping primers.
For inhibitor treatment experiments, embryos were dechorionated at 24 hpf. The clutches were divided evenly into two groups: vehicle control or inhibitor treatment groups. Prdm3 heterozygous intercrossed embryos were treated with either the Wnt inhibitor IWR-1 (Sigma-Aldrich), at a final concentration of 0.75 µM, or DMSO, in E2 embryo water. Prdm16 heterozygous intercrossed embryos were treated with either the Wnt activator Gsk inhibitor XV (CalbioChem), at a final concentration of 0.05 µM, or DMSO for vehicle control, in E2 embryo water. Inhibitor- or vehicle-containing E2 embryo water was removed after a 24-h window (at 48 hpf developmental time) and replaced with fresh E2 embryo water. At 6 dpf, larvae were collected and fixed in 2% paraformaldehyde (PFA) and subjected to Alcian Blue/Alizarin Red staining for cartilage and bone assessment.
For zebrafish, Alcian Blue (cartilage) and Alizarin Red (bone) staining was performed at room temperature as previously described (Walker and Kimmel, 2007). Briefly, 6 dpf larvae were collected and fixed for 1 h in 2% PFA. Following a 10-min rinse in 100 mM Tris (pH 7.5) and 10 mM MgCl2, larvae were incubated in Alcian Blue solution [0.04% Alcian Blue, 80% ethanol, 100 mM Tris (pH 7.5), 10 mM MgCl2] overnight at room temperature. Larvae were then rehydrated and de-stained through a gradient of ethanol solutions [80%, 50%, 25% ethanol containing 100 mM Tris (pH 7.5) and 10 mM MgCl2], then bleached for 10 min in 3% H2O2 with 0.55% KOH at room temperature, washed twice in 25% glycerol with 0.1% KOH, then stained in Alizarin Red [0.01% Alizarin Red dissolved in 25% glycerol and 100 mM Tris (pH 7.5)] for 30-45 min at room temperature. Samples were de-stained in 50% glycerol with 0.1% KOH. Whole-mount and dissected and flat-mounted specimens were mounted in 50% glycerol and imaged with LAS v4.4 software on a Leica M165 FC stereomicroscope. High-magnification images of chondrocytes were imaged with LAS v4.4 software on an Olympus BX51 WI compound microscope.
For mice, Alcian Blue and Alizarin Red staining was performed as previously described for E18.5 embryos (Wallin et al., 1994; Shull et al., 2020). Briefly, mouse embryos were harvested at E18.5 in 1×PBS. Skin and internal organs were removed, and specimens were fixed overnight in 95% ethanol at room temperature followed by an incubation in 100% acetone for 2 days at room temperature. Embryos were then stained in Alcian Blue/Alizarin Red staining solution (0.015% Alcian Blue, 0.05% Alizarin Red, 5% glacial acetic acid and 70% ethanol) for 3 days at 37°C. Stained embryos were rinsed in water before undergoing an initial clearing in 1% KOH overnight at room temperature, followed by a gradient series of decreasing KOH concentrations and increasing glycerol concentrations. Skeletal preparations were stored and imaged in 80% glycerol on a Leica M165 FC stereomicroscope with LASX v4.4 software.
Mouse embryos were collected at E14.5 in 1×PBS. The mandible was dissected and removed from the heads of the animals and fixed in 4% PFA, cryoprotected in 30% sucrose solution, embedded in OCT embedding medium, sectioned to a thickness of 8 µm and mounted onto glass slides with a Leica CM1520 cryostat. For staining, sections were brought to room temperature and rehydrated in 1×PBS before staining with Weigert's iron Hematoxylin, 0.05% Fast Green and 0.1% Safranin O and mounted with Permount (Electron Microscopy Sciences). Stained sections were imaged on an Olympus BX51 WI compound microscope with LASX v4.4 software.
FACS of neural crest cells
At 48 hpf, sox10:EGFP-positive embryos were stage-matched and selected under a fluorescence dissecting microscope. prdm3−/−;Tg(sox10:EGFP) and prdm16−/−;Tg(sox10:EGFP) single-mutant embryos were identified based on their phenotype (small domed heads and reduced eye pigmentation) at 48 hpf and extensively validated and confirmed by genotyping so that mutant embryos could be identified for FACS and subsequent sequencing experiments (Fig. S9). To prepare single-cell suspensions for FACS, 30-40 embryos of each genotype were dechorionated and rinsed in 1× DPBS (Ca2+ and Mg2+ free). Embryos were then dissociated in Accumax (Innovative Cell Technologies, AM-105) containing DNaseI (Roche Diagnostics). The samples were incubated at 31°C and agitated by pipetting every 10-15 min for 1 h to promote cell dissociation. Following the digest, a wash solution (1× D-PBS and DNaseI) was added to stop the reaction. Cells were filtered through a 70 µm nylon mesh strainer, centrifuged at 2000 rpm (376 g) for 5 min at 4°C and resuspended in FACS basic sorting buffer [1 mM EDTA, 25 mM HEPES (pH 7.0), 1% fetal bovine serum in 1× D-PBS]. Cell suspensions were stained with DAPI (1:1000) and kept on ice. GFP-positive cells were sorted on a MoFlo XDP100 cell sorter (Beckman-Coulter) and collected in 1× DPBS. Following FACS, GFP-positive cells were then processed for RNA-seq or ATAC-seq.
Following FACS, GFP-positive NCCs were centrifuged briefly at 2000 rpm (376 g) for 5 min and resuspended in TRIzol LS lysis reagent (Invitrogen/Life Technologies). Total RNA was extracted from sorted cells using chloroform extraction and the Direct-zol RNA miniprep kit (Zymo Research) according to the manufacturer's instructions. Purified RNA quality and quantity was assessed on a High Sensitivity RNA Screen Tape (Agilent Technologies) and Infinite M200pro plate reader (Tecan). cDNA libraries were generated using the Clontech Pico Library Prep Kit. Following library generation, sequencing was performed on an Illumina NovaSEQ 6000 system to a depth of ∼50 million reads. Library construction and sequencing was performed at the University of Colorado Anschutz Medical Campus Genomics and Microarray Core Facility.
RNA isolation and RT-qPCR
For zebrafish, whole heads were dissected and removed from wild-type, prdm3−/− and prdm16−/− embryos at 48 hpf. Five to seven embryos of the same genotype were pooled and lysed in TRIzol LS lysis reagent (Invitrogen/Life Technologies). Total RNA was isolated following a chloroform extraction. For mice, mandibular processes (MdPs) were dissected on ice from three independent replicates of E11.5 Prdm3fl/fl;Wnt1-Cre+/Tg, Prdm16fl/fl;Wnt1-Cre+/Tg and control (Prdm3fl/fl or Prdm16fl/fl) embryos. The overlying ectoderm was removed by digestion in 0.25% trypsin for 10-15 min on ice. MdPs were rinsed in 10% FBS for 1 min before a quick rinse in 1×PBS. MdPs were lysed in TRIzol LS (Invitrogen/Life Technologies) and total RNA was isolated from these samples using the Direct-zol RNA miniprep kit (Zymo Research) according to the manufacturer's instructions.
For RT-qPCR, total RNA was isolated as described above for zebrafish and mouse MdPs and (0.5-1.0 µg) was reverse transcribed to cDNA with SuperScript III First-Strand Synthesis cDNA kit (Invitrogen/Life Technologies) for real-time semiquantitative PCR (RT-qPCR) with primers (Table S5) and SYBR Green Master Mix (Bio-Rad). Transcript levels were normalized to the reference gene, gapdh (zebrafish) or Actb (mouse). Transcript abundance and relative gene expression were quantified using the 2−ΔΔCt method relative to control.
ATAC-seq was performed as previously described (Buenrostro et al., 2015; Liu et al., 2020). Briefly, following tissue dissociation FACS (as described above), 25,000 GFP-positive cells per genotype were pelleted at 500 g for 5 min at 4°C. The cell pellets were gently resuspended in cold lysis buffer (10 mM Tris-HCl, pH 7.5, 10 mM NaCl, 3 mM MgCl2, 0.1% v/v NP-40). Lysed cells were immediately centrifuged at 500 g for 20 min at 4°C. Cell nuclei pellets were resuspended in tagmentation mix with Nextera Tagmentation Buffer (Illumina) and Nextera Tagmentation DNA Enzyme (Illumina). Tagmentation reactions were adjusted to accommodate 25,000 cells and incubated at 37°C for 30 min, mixing at 10 min intervals. Tagmented DNA was purified using the Zymogen DNA Clean and Concentrator kit (Zymo Research). Libraries were amplified and indexed using NEBNext High Fidelity 2× PCR master mix (New England Biolabs). Following 11 cycles of amplification, libraries were purified with AmpureXP beads (Beckman Coulter), quantified with Qubit and subjected to sequencing on the Illumina NovaSEQ 6000 system at a depth of ∼50 million reads at the University of Colorado Anschutz Medical Campus Genomics and Microarray Core Facility. ATAC-seq experiments were performed in duplicate for two biological replicates per genotype.
CUT&RUN paired with RT-qPCR
CUT&RUN was performed on wild-type whole embryos as described (Skene and Henikoff, 2017). Briefly, ∼200 48 hpf wild-type zebrafish embryos were pooled and dissociated to a single-cell suspension using Accumax (Innovative Cell Technologies) and DNaseI (Roche Diagnostics) for 1 h with gentle pipetting to agitate the tissue every 10 min. Following dissociation, a wash solution containing 1×PBS and DNaseI was added to stop the reaction. Cells were passed through a 40 µm filter and counted then 500,000 cells were incubated on activated Concanavalin A conjugated paramagnetic beads (Epicypher) for 10 min at room temperature. Cell-bound beads were resuspended in antibody buffer [20 mM HEPES, pH 7.5; 150 mM NaCl; 0.5 mM Spermidine (Invitrogen); 1× Complete-Mini Protease Inhibitor tablet (Roche Diagnostics); 0.01% Digitonin (Sigma-Aldrich); 2 mM EDTA] and incubated in the corresponding antibodies (validated to work in zebrafish; see Fig. S7A,B). Although validated, we cannot rule out the possibility that these antibodies may to a lesser extent recognize other proteins of similar size to the wild-type band. Antibodies [IgG (Jackson ImmunoResearch, 111-005-003, RRID: AB_2337913), Prdm3/Evi1 (Abcam, ab28457, RRID: AB_732271), Prdm16 (antibody gifted from Patrick Seale, University of Pennsylvania, USA; R&D Systems, AF6295, RRID: AB_10717965) and H3K27ac (Cell Signaling Technology, 4353S, RRID: AB10545273)] were added to samples and incubated overnight with rotation at 4°C. The following day, cells were washed in Digitonin Buffer wash solution [20 mM HEPES, pH 7.5; 150 mM NaCl; 0.5 mM Spermidine (Sigma-Aldrich); 1× Complete-Mini Protease Inhibitor tablet (Roche Diagnostics); 0.01% Digitonin (Sigma-Aldrich)] twice and then incubated with pAG-MNase (EpiCypher) for 10 min at room temperature. Following another two washes in Digitonin Buffer, 1 µl of 100 mM CaCl2 was added to each sample and incubated at 4°C for 2 h. This digestion reaction was stopped with the addition of Stop Buffer (340 mM NaCl, 20 mM EDTA, 4 mM EGTA, 50 µg/ml RNaseA and 50 µg/ml glycogen) for 10 min at 37°C. DNA fragments were purified using a DNA Clean and Concentrator Kit (Zymo Research). For RT-qPCR, eluted CUT&RUN fragmented DNA was amplified using the NEBNext Ultra II DNA Library Prep Kit (New England Biolabs) as per the manufacturer's instructions. RT-qPCR was performed with primers designed to flank putative Prdm3 and Prdm16 binding motifs at promoter regions of Wnt/β-catenin targets genes or negative control genes with no binding sites (gapdh) and SYBR Green Master Mix (Bio-Rad). The fold enrichment for the abundance of Prdm3, Prdm16 or H3K27ac of those amplified regions was calculated and normalized relative to the IgG control and averaged across three different experiments.
Zebrafish larvae were collected at the indicated time points and fixed in 4% PFA overnight at 4°C. Following fixation, embryos were washed in 1× PBS (pH 7.3) with 1% Triton X-100 three times for 10 min at room temperature. For antigen retrieval, embryos were incubated in 1 µg/ml Proteinase K diluted in 1× PBS with 1% Triton X-100 for 20 min at room temperature. Following proteinase K treatment, embryos were incubated in 4% PFA for 15 min at room temperature, then washed three more times in 1× PBS with 1% Triton X-100 for 10 min at room temperature. Embryos were blocked at room temperature for 1 h in blocking solution containing 10% normal goat serum and 1% bovine serum albumin in 1× PBS. Samples were incubated in primary antibodies diluted (1:100) in blocking solution [anti-acetylated α-tubulin (Sigma-Aldrich, T6793, RRID: AB_477585), anti-phosphorylated Y489 β-catenin (Developmental Studies Hybridoma Bank, PY489-B-catenin, RRID: AB_732271), Rhodamine phalloidin (Thermo Fisher, R415, RRID: AB_2572408), Wheat Germ Agglutinin (Life Technologies/Invitrogen, W32466)], overnight at 4°C. Following primary antibody incubation, samples were thoroughly washed in PBS with 1% Triton X-100 before incubating with corresponding secondary antibodies overnight at 4°C, then were washed again thoroughly in 1× PBS with 1% Triton X-100 before incubation with DAPI diluted in 1× PBS for 1 h at room temperature. Samples were quickly washed in PBS with 1% Triton X-100 before being mounted in Vectashield mounting media (Invitrogen) on glass slides. Embryos were imaged on a Leica TCS SP8 confocal microscope. Images were processed using LASX software and ImageJ.
Whole-mount immunohistochemistry on mouse embryos was performed as previously described (Joyner and Wall, 2008). Briefly, E10.5 or E11.5 wild-type embryos were dissected in 1× PBS and fixed in methanol/DMSO (4:1) at 4°C overnight. Embryos were then transferred to methanol/DMSO/H2O2 (4:1:1) and incubated at room temperature for 8 h. Embryos were rehydrated at room temperature through a series of washes: 50% methanol for 30 min, 1×PBS for 30 min and PBSMT (2% nonfat instant skim milk, 0.5% Triton X-100 in 1×PBS) for 1 h. Specimens were then incubated in primary antibodies diluted 1:200 in PBSMT [Prdm3 (Abcam, ab28457, RRID: AB_732271), Prdm16 (Patrick Seale, University of Pennsylvania; R&D Systems, AF6295, RRID: AB_10717965)] overnight at 4°C with rocking. Embryos were then washed. Samples were incubated in secondary antibodies diluted 1:500 in PBSMT [HRP-conjugated goat anti-rabbit IgG (Santa Cruz Biotechnology, sc-2004, RRID: AB_631746)] overnight at 4°C with rocking. Embryos were washed again in PBSMT twice, 1 h each, at 4°C followed by three times, 1 h each, at room temperature and then rinsed and washed in PBTX (0.2% w/v bovine serum albumin, 0.5% v/v Triton X-100, 1×PBS) for 20 min. Samples were then incubated in DAB-NiCl2 (3,3′-diaminobenzidine tetrahydrochloride, NiCl2 in 1×PBS) for 30 min at room temperature. Hydrogen peroxide was added to a final concentration of 0.03% and the samples were rocked until the color developed (5 min), then post-fixed in 4% PFA. Samples were then rinsed through a methanol gradient and then incubated in BABB (benzyl benzoate/benzyl alcohol 2:1) for 10 min. Stained embryos were imaged on a Leica M165 FC stereomicroscope with LASX v4.4 software.
Zebrafish embryos (30-40) of each genotype were pooled at 48 hpf and incubated on ice for 5 min. Calcium-free Ginzberg Fish Ringer's solution (55 mM NaCl, 1.8 mM KCl, 1.25 mM NaHCO3) was added to the embryos for deyolking. After pelleting samples for 1.5 min at 5000 rpm (2348 g), samples were washed in deyolking wash buffer [110 mM NaCl, 3.5 mM KCl, 2.7 mM CaCl2, 10 mM Tris (pH 8.5)]. Samples were pelleted again for 1.5 min at 5000 rpm (2348 g). All liquid was removed and the pellet was resuspended in SDS lysis buffer [0.1% glycerol, 0.01% SDS, 0.1 M Tris (pH 6.8)] for 10 min on ice. Embryos were homogenized in lysis buffer and again briefly centrifuged at 2348 g for 5 min. Total protein concentrations were determined using the Bio-Rad Dc Protein Assay (Bio-Rad). Proteins (20 µg) were separated by SDS-PAGE (12%) and transferred to polyvinylidene difluoride membranes. Membranes were blotted using antibodies for Prdm3 (Abcam, ab28457, RRID: 732271), Prdm16 (antibody gifted from Patrick Seale Lab, University of Pennsylvania, USA; R&D Systems, AF6295, RRID: AB_10717965) and total H3 (Cell Signaling Technology, 9715S, RRID: AB_331563) diluted at 1:1000 and corresponding secondary antibodies. Chemiluminescence detection was performed with Immobilon Forte Western HRP Substrate (Millipore) on a Bio-Rad Chemidoc multiplex imager.
Zebrafish embryos were imaged on a Leica DMi8 microscope equipped with an Andor Dragonfly 301 spinning disk confocal system. Approximately 100 µm z-stacks were captured at 0.35 µm intervals every 30 min for approximately 16 h. ImageJ was used for image processing.
Quantification and statistical analysis
RNA-seq bioinformatics analysis
Following trimming and read alignment, paired-end reads were mapped to the zebrafish genome (danRer11) assembly using TopHat (Trapnell et al., 2009, 2012). Differential expression between mutant and wild type was calculated using Cufflinks (Trapnell et al., 2012). Gene expression was expressed in fragments per kilobase of transcript per million mapped reads (FPKM). Complete RNA-seq datasets are available at the Gene Expression Omnibus repository (GSE175767). Normalized counts were converted to z-scores for plotting heatmaps using the pheatmaps R package (https://cran.r-project.org/web/packages/pheatmap/index.html). Gene lists were analyzed for functional annotation using GO enrichment analysis based on the PANTHER Classification System (Mi et al., 2019a,b; Ashburner et al., 2000; Harris et al., 2004).
ATAC-seq bioinformatics analysis
Adapters and barcodes were removed from paired end reads using Cutadapt (Martin, 2011), and trimmed paired-end sequencing reads were aligned to the zebrafish genome (danRer11) using Bowtie2 (v2.4.4) with default parameters (--end-to-end, --sensitive) (Langmead and Salzberg, 2012; Langmead et al., 2019). Peaks were called for each sample using the Genrich peak calling program (https://github.com/jsh58/Genrich), with the following parameters: -r, -p0.01, -a200, -j, -e MT. Intersection peaksets were generated for each genotype (e.g. Prdm3 mutant, Prdm3 wild type, Prdm16 mutant, Prdm16 wild type) by using bedtools (v.2.30.0) or bedops (v.2.4.37) to identify the sets of peaks observed in both replicates of each genotype. DiffBind (v.3.4.0) was used to assess differentially accessible regions between genotypes and their respective sibling controls using the default parameters [false discovery rate (FDR)≤0.05] (Ross-Innes et al., 2012). For volcano plots, −log(pvalue) was used for visualization of the data. Annotation of peaks corresponding to transcription start site (TSS)/promoter, intergenic, intronic and transcription end site (TES) locations was carried out using HOMER (v4.11) annotatePeaks.pl script (Heinz et al., 2010). DeepTools (v2.0) was used to generate bigWig coverage files for visualization (normalization=counts per million) (Ramirez et al., 2016). TOBIAS was performed following the standard pipeline and workflow (Bentsen et al., 2020). Complete ATAC-seq datasets are available in the GEO repository (GSE175767).
To quantify chondrocyte organization from Alcian-stained and dissected flat mounts, high-magnification images of cartilage elements were imported into ImageJ. An angle was drawn from the center of three adjacent chondrocytes moving along the cartilage element in the direction of growth of that structure. Angles were measured and averaged across at least five cartilage elements (namely the posterior ceratobranchial cartilages) per individual. All possible combinations were measured within the field of view to eliminate any possible cell-selection bias.
To quantify chondrocyte cell orientation, individual chondrocytes in the developing palatoquadrate were divided into quadrants in ImageJ. The positioning of acetylated α-tubulin puncta was indicated as either 0°, 90°, 180° or 270° in the direction of growth of the palatoquadrate, i.e. anteriorly toward the jaw joint junction with the Meckel's cartilage. Positioning was tracked through z-stack images and the number of acetylated α-tubulin puncta in each quadrant was tabulated for each individual and normalized to the total number of cells analyzed for that individual. Total counts were collected for at least five individuals per genotype then averaged across genotypes and plotted as circular graphs. Nuclear β-catenin puncta (phosphorylated Y489 β-catenin) were quantified and tracked across ten chondrocytes in the palatoquadrate through z-stack images for one individual. The number of puncta was averaged for each individual and at least five individuals were analyzed per genotype.
For quantifying chondrocyte cell area in mouse Meckel's cartilage, the area of 200 cells across four or five sections anterior to posterior through the tissue were measured and averaged across three individuals per genotype. For cell numbers, cells were counted in a designated region of tissue area across four or five sections anteriorly to posteriorly throughout the tissue and averaged across three individuals per genotype.
Data shown are mean±s.d. from the number samples or experiments indicated in the figure legends. All assays were repeated at least three times with independent samples. P-values were determined with unpaired, two-tailed Student's t-tests.
We thank members of the Artinger laboratory for their thoughtful discussions on this project and Dr Richard Dorsky for the Tg(7xTCF-Xla.Siam:NLS-mCherry) line. We would also like to acknowledge undergraduates Oscar Yip and Joey Gerlach for their contributions to this work. We thank Kristen D'Elia and Drs Jeremy Dasen and David Schoppik at New York University for their feedback and insightful comments on this project. Special thanks to Dr David Clouthier for mouse craniofacial anatomy consultation. We thank the Denver zebrafish community, and the zebrafish and mouse facility staff for excellent animal care. The Functional Genomics and Biostatistics and Bioinformatics were funded by National Cancer Institute shared resource (P30CA046934 to the University of Colorado Cancer Center).
Conceptualization: L.C.S., K.B.A.; Methodology: L.C.S., K.B.A.; Formal analysis: L.C.S., E.S.L., H.M.K., J.C.C., K.J.; Investigation: L.C.S., E.S.L.; Resources: S.G., M.K.; Writing - original draft: L.C.S., K.B.A.; Writing - review & editing: L.C.S., E.S.L., J.C.C., K.B.A.; Supervision: K.B.A.; Funding acquisition: L.C.S., K.B.A.
Funding was provided by the National Institute of Dental and Craniofacial Research (R01 DE024034 to K.B.A. and F32 DE029099 to L.C.S.) and the National Institute of Neurological Disorders and Stroke (P30 NS048154 to the UC Anschutz Medical Campus zebrafish core facility). Deposited in PMC for release after 12 months.
The RNAseq and ATACseq datasets generated in this paper are available at Gene Expression Omnibus under accession number GSE175767.
Peer review history
The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.200082.
The authors declare no competing or financial interests.