The Gq/11 family of Gα subunits is necessary and sufficient for lower jaw development

Endothelin A receptor (Ednra) and peptide ligand Endothelin 1 (Edn1) are essential for patterning the embryonic tissue that gives rise to lower jaw structures (Clouthier et al., 1998; Yanagisawa et al., 1998b; Miller et al., 2000; Ozeki et al., 2004). During development, Ednra is expressed in migratory and post-migratory cranial neural crest cells (NCCs) (Kurihara et al., 1995; Clouthier et al., 1998; Yanagisawa et al., 1998a; Miller et al., 2000; Nair et al., 2007; Ruest and Clouthier, 2009), which are multipotent progenitors that give rise to all bone and cartilage structures of the face (Le Douarin and Kalcheim, 1999; Chai et al., 2000; Tang and Bronner, 2020). Post-migratory cranial NCCs in pharyngeal arch 1 (or the mandibular portion of arch 1 in mammals) are stimulated by Edn1 secreted by the adjacent ventral ectoderm (Clouthier et al., 1998; Yanagisawa et al., 1998a; Miller et al., 2000), inducing gene expression that establishes two distinct patterning domains for lower jaw structures: a ventral patterning domain that gives rise to the mandible and Meckel's cartilage and an intermediate patterning domain that gives rise to the jaw joint in zebrafish and middle ear structures in mammals (Clouthier et al., 2000; Miller et al., 2003; Ozeki et al., 2004; Jeong et al., 2008; Ruest and Clouthier, 2009; Talbot et al., 2010; Barske et al., 2016; Askary et al., 2017; Tavares et al., 2017). Edn1/Ednra also antagonizes the effects of Nr2f nuclear receptors and Jagged/Notch, two factors that establish upper jaw identity of cranial NCCs in the dorsal domain of the first and second arches (Zuniga et al., 2010; Barske et al., 2016, 2018; Teng et al., 2017). Thus, Edn1/Ednra plays the dual role of establishing lower jaw identity of cranial NCCs while also establishing the boundary separating the lower jaw and upper jaw patterning domains.

In mice and zebrafish, attenuation of the Edn1/Ednra signaling pathway diminishes intermediate and ventral patterning gene expression, causing expansion of dorsal patterning gene expression into the intermediate and ventral domains. This results in a homeotic transformation of lower jaw structures into upper jaw-like structures (Miller et al., 2000; Ozeki et al., 2004; Ruest and Clouthier, 2009; Pritchard et al., 2020; Kanai et al., 2022). Conversely, ectopic activation of the Edn1/Ednra signaling pathway in the dorsal patterning domain (or the maxillary portion of arch 1 in mammals) results in homeotic transformation of upper jaw structures into lower jaw-like structures (Sato et al., 2008b; Alexander et al., 2011; Zuniga et al., 2011; Tavares and Clouthier, 2015; Barske et al., 2018; Kurihara et al., 2023). In either case in which Edn1/Ednra signaling is lost or ectopically activated, misexpression of patterning genes that are positively or negatively regulated by Edn1/Ednra changes the positional identities of cranial NCCs. While the Edn1/Ednra-regulated genes have been extensively characterized, the intracellular signaling events connecting Ednra activation to gene regulation are not fully understood.

Ednra signaling is mediated by heterotrimeric G proteins, complexes composed of a Gα subunit and a Gβ/Gγ obligate heterodimer (Gβγ) (Gilman, 1987). Ligand-bound Ednra facilitates the exchange of GDP for GTP on the Gα subunit, causing Gα and Gβγ to dissociate and interact with their respective signaling effectors (Oldham and Hamm, 2008). Ednra can couple with all four Gα subunit family members – Gq/11, G12/13, Gs and Gi/o – to varying extents depending on the cell and tissue type (Aramori and Nakanishi, 1992; Khac et al., 1994; Kawanabe et al., 2002; Wettschureck and Offermanns, 2005; Inoue et al., 2019; Okashah et al., 2019; Masuho et al., 2023). In mice, Gq and G11 proteins are encoded by Gnaq and Gna11, respectively. Removal of Gnaq and Gna11 using conventional knockout alleles results in midgestational lethality prior to the development of craniofacial structures, although embryos with one wild-type allele for Gnaq (Gnaq^+/−;Gna11^−/−) or Gna11 (Gnaq^−/−;Gna11^+/−) can survive to term (Offermanns et al., 1998). Notably, Gnaq^−/−;Gna11^+/− embryos exhibit defects that appear to be confined to the jaw joint and middle ear structures, while Gnaq^+/−;Gna11^−/− embryos display no craniofacial defects (Offermanns et al., 1998). This suggests that, although both Gnaq and Gna11 are required for early development, Gnaq is specifically required for patterning the intermediate domain. Similarly, embryos with a combination of conventional knockout alleles for Gna11 and NCC-specific conditional knockout alleles for Gnaq (Gnaq^flox/flox;Gna11^−/−;P0-Cre) exhibited craniofacial defects that appeared to be confined to proximal lower jaw structures, including fusion of the jaw joint and hypoplasia of the proximal mandible, tympanic ring, malleus and incus (Dettlaff-Swiercz et al., 2005). In contrast, the entire mandible undergoes a homeotic transformation into a maxilla-like structure in mice lacking Edn1, Ednra or Ece1, a gene that encodes an essential enzyme in the Edn1 biosynthetic pathway (Yanagisawa et al., 1998b; Ozeki et al., 2004; Ruest and Clouthier, 2009). These differences suggest that there are Gq/11-dependent and -independent subdomains within the mandibular arch, with patterning of the Gq/11-independent domain mediated by a different Gα family member (Sato et al., 2008a). Arguing against this model, mice homozygous for knockout alleles of Gs, G12/13 or Gi/o do not exhibit patterning defects in distal lower jaw structures (Jiang et al., 1998; Dettlaff-Swiercz et al., 2005; Wettschureck and Offermanns, 2005; Plummer et al., 2012; Lei et al., 2016). Thus, the question remains whether patterning is mediated solely by Gq/11 or a combination of Gq/11 and other Gα family members.

This question has direct significance for human health. Perturbations to the Ednra signaling pathway are the underlying causes of congenital syndromic disorders, such as mandibulofacial dysostosis with alopecia (Gordon et al., 2015; Kurihara et al., 2023), oro-oto-cardiac syndrome (Pritchard et al., 2020) and auriculo-condylar syndrome (Rieder et al., 2012; Gordon et al., 2013; Marivin et al., 2016; Kanai et al., 2022). Thus, clarification of the G protein family or families mediating Edn1/Ednra-dependent jaw patterning could help elucidate additional signaling pathway components that are linked to human craniofacial differences.

In this study, we examined the role of Gq/11 in craniofacial patterning using multiple orthogonal approaches in zebrafish. Using a Gq/11 small molecule inhibitor YM-254890 (YM) (Takasaki et al., 2004; Nishimura et al., 2010), we demonstrate that YM-treated zebrafish embryos exhibit craniofacial phenotypes and gene expression changes indicative of Ednra signaling attenuation. In addition, we created mutant alleles for three genes that encode zebrafish Gq/11 proteins, gnaq, gna11a and gna11b, and show that triple homozygous mutant larvae exhibit craniofacial phenotypes nearly identical to those observed in edn1^−/− mutant larvae. Furthermore, the craniofacial phenotypes in edn1^−/− mutant embryos can be rescued by inducing expression of a constitutively active form of Gq during craniofacial patterning. Taken together, these results suggest that Gq/11 is necessary and sufficient for establishing intermediate and ventral patterning domains and the subsequent development of all Ednra-dependent skeletal elements of the lower jaw in zebrafish.

To begin our analysis of the G protein signaling pathway downstream of Ednra in zebrafish craniofacial development, we examined the effect of treating zebrafish embryos with YM, a naturally derived compound from Chromobacterium species that inhibits a subset of Gq/11 proteins with high selectivity and potency (Takasaki et al., 2004; Nishimura et al., 2010). YM binds to a hydrophobic cleft in mammalian Gq in a region spanning the GTPase (Ras-like) and α-helical domains, stabilizing the allosteric mechanism that promotes GDP release, the rate-limiting step for Gα activation (Nishimura et al., 2010; Flock et al., 2015). This inhibitory effect is dependent on a network of hydrophobic interactions between YM and eight key residues that are present in Gq, G11 and G14, but not in G15 or other Gα family members (Fig. 1) (Nishimura et al., 2010; Onken et al., 2018). In zebrafish, these eight residues are highly conserved in Gq, G11 and G14, weakly conserved in G15 and highly divergent in other Gα family members (Fig. 1), suggesting that YM should display highly specific inhibitory effects towards the Gq/11 signaling pathway in zebrafish, except for G15. Because the expression of G15 is restricted to sensory cells (Ohmoto et al., 2011; Oka and Korsching, 2011), it is unlikely that G15 signaling activity would interfere with our analysis of craniofacial development.

To test whether YM treatment would produce craniofacial phenotypes, wild-type embryos were incubated with YM or solvent control [dimethyl sulfoxide (DMSO)] from 16 to 36 h post-fertilization (hpf), the approximate window for Edn1/Ednra-dependent craniofacial patterning (Miller et al., 2000, 2003; Walker et al., 2006; Kimmel et al., 2007; Ruest and Clouthier, 2009; Vieux-Rochas et al., 2010; Alexander et al., 2011; Zuniga et al., 2011; Barske et al., 2016, 2018; Meinecke et al., 2018), and then grown to 6 days post-fertilization (dpf) and processed for bone and cartilage staining. We tested a wide range of YM concentrations, from 1 nM to 100 µM. To ensure that the DMSO concentration in the embryo media did not exceed 1% (the upper limit for DMSO that does not cause developmental defects in zebrafish; Hoyberghs et al., 2021), 100 μM YM was the highest concentration we could use (see also Materials and Methods). Compared to control embryos (Fig. 2A,C,E,H), only 100 μM YM produced craniofacial phenotypes, which included an open mouth and uninflated swim bladder (Fig. 2B,D), all indications of defective jaw function (Neuhauss et al., 1996). Additional Edn1/Ednra pathway-associated defects were further visible in flat-mounted preparations of dissected viscerocrania, including hypoplasia of the palatoquadrate and symplectic cartilages, fusion of the jaw and hyomandibular joints, and fusion of the opercle and branchiostegal ray dermal bones (Fig. 2F,I). Defects to at least one of these structures were found in 93% (26/28) of wild-type larvae treated with YM, with considerable variation in the laterality and number of skeletal elements affected per larvae (Fig. 2K).

To determine whether the phenotypes resulting from YM treatment represented a partial or complete loss of Edn1/Ednra signaling, we first generated a new edn1 allele (edn1^co3009) lacking the majority of the coding sequence for the Edn1 bioactive peptide (Fig. S1). The original edn1 mutant allele sucker results from an amino acid substitution (p.Asp8Val) in the bioactive portion of Edn1 (Miller et al., 2000) (Fig. S1). While amino acid substitutions at p.Asp8 (e.g. p.Asp8Asn, p.Asp8Ala) result in approximately 99% less vasoconstrictor activity than that observed for wild-type Edn1 in vitro (Huggins et al., 1993), the activity of Edn1 encoded by the sucker allele has never been formally tested. Larvae homozygous for the edn1^co3009 allele exhibited absence and severe hypoplasia of ventral structures (Fig. 2G,J) and recapitulated the sucker craniofacial phenotype (Fig. S1), indicating that the sucker and edn1^co3009 craniofacial phenotypes are caused by a complete loss of Edn1 activity. In contrast, the phenotypes of YM-treated larvae more closely resembled phenotypes caused by partial loss of Edn1/Ednra signaling (Walker et al., 2006; Kimmel et al., 2007; Miller et al., 2007; Talbot et al., 2010). These data suggest that YM-induced defects are caused by partial inhibition of the Edn1/Ednra signaling pathway.

To determine the developmental time during which Gq/11 activity is required for lower jaw development, embryos were treated with 100 µM YM across four overlapping time intervals (16-24, 16-36, 24-48 or 36-48 hpf) and then analyzed for craniofacial skeleton defects at 6 dpf. Craniofacial defects were observed only in embryos treated from 16 to 36 hpf, or between 24 and 48 hpf (Fig. 2L), with the percentage of affected larvae (Fig. 2L) and frequency and severity of affected skeletal elements (Fig. S2) being similar for both time intervals. These results suggest that Gq/11 activity is required for patterning between 24 and 36 hpf. This window of sensitivity to YM is similar to the window of sufficiency for Edn1/Ednra signaling in craniofacial patterning in zebrafish and mice (Miller et al., 2000, 2003; Walker et al., 2006; Kimmel et al., 2007; Ruest and Clouthier, 2009; Vieux-Rochas et al., 2010; Alexander et al., 2011; Zuniga et al., 2011; Barske et al., 2016, 2018; Meinecke et al., 2018).

The loss of Edn1/Ednra signaling causes characteristic changes to gene expression in NCCs of the first and second pharyngeal arches, including downregulation of lower jaw patterning genes and upregulation of upper jaw patterning genes in the intermediate and ventral domains (Clouthier et al., 1998, 2000; Thomas et al., 1998; Miller et al., 2000, 2003; Ozeki et al., 2004; Ruest et al., 2004; Walker et al., 2006, 2007; Nair et al., 2007; Sato et al., 2008a; Tavares et al., 2012, 2017; Barske et al., 2016, 2018; Askary et al., 2017). To determine whether YM caused similar changes to gene expression in cranial NCCs as does partial loss of Edn1/Ednra signaling, we performed single-cell RNA sequencing (scRNA-seq) on cranial NCCs isolated from Tg(sox10:mRFP;fli1a:EGFP) embryos treated with DMSO or 100 µM YM from 16 to 36 hpf (Fig. 3A). Fluorescence-activated cell sorting (FACS) of sox10:mRFP;fli1a:EGFP double-positive cells allows purification of cranial NCC populations in zebrafish embryos (Askary et al., 2017; Mitchell et al., 2021). Subsequent analyses with Seurat identified 12 cell clusters with unique transcriptional profiles (Fig. 3B, Fig. S3, Table S1). Using known marker genes (Miller et al., 2003; Askary et al., 2017; Barske et al., 2018; Mitchell et al., 2021; Fabian et al., 2022; Stenzel et al., 2022), NCC populations were identified from the frontonasal region (akap12b⁺, alx4a⁺), anterior pharyngeal arches 1 and 2 (dlx2a⁺, dlx5a⁺) and posterior pharyngeal arches 3-7 (prdm1a⁺) (Fig. 3B,D). The NCCs of pharyngeal arches 1 and 2 were distributed across three clusters, with one cluster representing the ventral domain cells of both arches 1 (hand2⁺) and 2 (hand2⁺, hoxb2a⁺), and two separate clusters representing dorsal domain cells of arch 1 (pou3f3b⁺, Hox−) and arch 2 (jag1b⁺, hoxb2a⁺). Although we did not find a dedicated cluster for intermediate domain cells, we used known marker genes to identify the approximate population of intermediate domain cells nested within the ventral and dorsal clusters (Fig. 3C, Fig. S4). The remaining clusters appeared to consist of NCC-derived pigment cells (cluster 11: dct⁺, tyrp1b⁺) and hematopoietic derivatives (cluster 9: hbbe3⁺, hbae3⁺; cluster 10: myb⁺, lyz⁺); these clusters were not examined further in this study.

To examine the effects of YM on gene expression, we performed differential expression analysis between equivalent clusters in YM-treated samples relative to control samples (Table S2). Differentially expressed genes with adjusted P-values less than 0.05 were found in six out of eleven clusters – frontonasal, ventral arches 1 and 2, dorsal arch 1, dorsal arch 2, posterior arches 3-7 and cluster 5 (unknown identity) (Fig. S5, Table S2). Overall, a total of 36 and 34 unique, non-overlapping genes were upregulated or downregulated, respectively, in YM-treated samples. Of the aberrantly upregulated genes, dorsal patterning genes including, hey1, her6, her9 and nr2f5, were present in the ventral and dorsal arch clusters (Fig. 3G, Fig. S5, Table S2) (Zuniga et al., 2010; Barske et al., 2016, 2018; Askary et al., 2017). The majority of aberrantly downregulated genes were in the ventral arches 1 and 2 and dorsal arch 1. In the ventral arch cluster, more than half of the reduced genes were Ednra targets. Notably, the intermediate patterning genes dlx4a and dlx4b (Fig. 3E,H-I′, Fig. S5) were downregulated to a greater extent than ventral patterning genes such as hand2 and dlx5a (Fig. 3F,J-K′, Fig. S5). These gene expression changes are similar to changes observed in animal models with partial loss of Edn1/Ednra signaling (Walker et al., 2006, 2007; Miller et al., 2007; Sato et al., 2008a; Ruest and Clouthier, 2009; Tavares et al., 2012). Furthermore, the relatively modest downregulation of ventral patterning genes compared with intermediate patterning is consistent with the unperturbed ventral domain-derived skeletal structures in YM-treated larvae (Fig. 2F,I). These results indicate that YM partially attenuates the Edn1/Ednra signaling pathway in cranial NCCs through the inhibition of Gq/11.

Because we could not use a YM concentration greater than 100 µM, we could not distinguish whether the YM-generated phenotype was due to incomplete inhibition of Gq/11 activity or whether Gq/11 activity was only required in the intermediate domain. To address this question, we tested the effects of YM on edn1^+/− zebrafish embryos, with the supposition that these embryos would have reduced expression of Ednra-dependent genes similar to Edn1^+/− mouse embryos (Vieux-Rochas et al., 2010) and thus would be sensitized to the effects of YM. Embryos generated from edn1^+/− mating pairs were treated with 100 µM YM or DMSO from 16 to 36 hpf, then processed for skeletal preparations at 6 dpf. DMSO treatment did not cause defects to any skeletal elements in either wild-type or edn1^+/− larvae (Fig. 4A,B). In YM-treated wild-type larvae, skeletal defects were confined to intermediate domain structures, with the ventral portion of the Meckel's cartilage and ceratohyal unaffected (Fig. 4D,G). However, in YM-treated edn1^+/− larvae, skeletal defects expanded to include ventral domain structures, shown by hypoplasia of the Meckel's cartilage and ceratohyal (Fig. 4E,G). Further, YM treatment resulted in more skeletal elements affected per individual in edn1^+/− larvae compared to wild-type larvae (Fig. 4H). YM treatment had no effect on the phenotype of edn1^−/− larvae (Fig. 4C,F). The increased severity and frequency of defects elicited by this pharmacogenetic interaction (Fig. 4G,H) suggests that YM can only partially inhibit Gq/11 activity in wild-type zebrafish embryos, whereas the inhibitory effect of YM on Gq/11 activity can be augmented in a edn1^+/− genetic background. These findings suggest that Gq/11 is required for patterning beyond the intermediate domain.

To interrogate more comprehensively the role of Gq/11 in facial development, we mutated genes encoding Gq and G11 proteins using CRISPR/Cas9. In zebrafish, Gq is encoded by gnaq, and G11 is encoded by two paralogs, gna11a and gna11b. Although gnaq expression in the embryonic head was undetectable by scRNA-seq (Fig. S6), quantitative PCR (qPCR; Fig. S7) or colorimetric in situ hybridization (Fig. S8), we proceeded to target the gnaq allele given the importance of Gnaq in mouse lower jaw development (Offermanns et al., 1998). Of the two G11 genes, gna11b expression was relatively higher than gna11a in scRNA-seq, qPCR and colorimetric in situ hybridization analyses (Figs S6-S8). Other Gq/11 family members, including gna14, gna14a and gna15.1-15.4, were not considered because they exhibit cell type-specific expression in sensory organs and hematopoietic lineages (Ohmoto et al., 2011; Oka and Korsching, 2011). We designed CRISPR/Cas9 guide RNAs to target exon 3 or 4 of the genes (Fig. S9A), similar to targeting strategies used for Gnaq and Gna11 knockout alleles in mice (Offermanns et al., 1998). For all three genes, we identified a frame-shift mutation that is predicted to result in a premature stop codon and non-functional protein product (Fig. S9B). Similar to wild-type larvae (Fig. 2), larvae and adults that were homozygous for a single mutant allele (gnaq^−/−, gna11a^−/− or gna11b^−/−) (Fig. S9) or triple heterozygous for all mutant alleles (gnaq^+/−;gna11a^+/−;gna11b^+/−) (Fig. 5I) did not exhibit overt craniofacial phenotypes and were viable and fertile. However, we observed a spectrum of craniofacial phenotypes in larvae generated from gnaq^+/−;gna11a^+/−;gna11b^+/− mating pairs (Fig. 5). A total of 336 larvae were processed for skeletal preparations at 6 dpf, followed by analysis of genotypes and phenotypes (Table S3). Flat-mounts of skeletal preparations for selected genotypes are shown in Fig. 5A-I. Quantification of the severity and frequency of defects for selected genotypes is shown in Fig. 5J. The most severe defects were observed in gnaq^−/−;gna11a^−/−;gna11b^−/− (Fig. 5A,J), gnaq^+/−;gna11a^−/−;gna11b^−/− (Fig. 5B,J) and gnaq^+/+;gna11a^−/−;gna11b^−/− (Fig. 5C,J) larvae, which exhibited defects in all Edn1/Ednra-dependent skeletal elements and resembled the edn1^−/− phenotype (Fig. 2G,J). This suggests that gna11a and gna11b are required for zebrafish lower jaw development, but gnaq is not. Larvae with one wild-type allele of either gna11a (gnaq^−/−;gna11a^+/−;gna11b^−/−) (Fig. 5D) or gna11b (gnaq^−/−;gna11a^−/−;gna11b^+/−) (Fig. 5E) exhibited phenotypes that were all less severe and more variable than that of gnaq^−/−;gna11a^−/−;gna11b^−/− larvae (Fig. 5A,J). Larvae with one wild-type allele of gna11a and gna11b (gnaq^−/−;gna11a^+/−;gna11b^+/−) (Fig. 5F) exhibited further reduced phenotypic severity and variation (Fig. 5J). Importantly, the addition of a wild-type gnaq allele to the above gna11a and gna11b allelic combinations – gnaq^+/−;gna11a^+/−;gna11b^−/− (Fig. 5G), gnaq^+/−;gna11a^−/−;gna11b^+/− (Fig. 5H) and gnaq^+/−;gna11a^+/−;gna11b^+/− (Fig. 5I) – had no effect on phenotypic severity and variation (Fig. 5J). Taken together, these results suggest that gna11a and gna11b are necessary for the development of Edn1/Ednra-dependent craniofacial structures in a gene-dosage dependent manner, whereas gnaq appears to be dispensable.

To determine whether Gq/11 signaling activity is sufficient for lower jaw development, we examined whether we could rescue the edn1^−/− mutant phenotype by expressing a constitutively active form of Gq, Gq-Q209L (Fig. 6) (Kalinec et al., 1992). The Q209L amino acid substitution prevents GTP hydrolysis and maintains Gq in an active signaling state independently of receptor activity (Kleuss et al., 1994). Gq and G11 proteins are functionally redundant (Malbon, 2005; Wettschureck and Offermanns, 2005), and the constitutively active mutants Gq-Q209L and G11-Q209L activate the same signaling pathways in the context of uveal melanoma (Annala et al., 2019; Lapadula et al., 2019). In order to control the timing of induction for Gq activity, we made a transgenic line that expresses Gq-Q209L under the regulation of a hsp70l heat shock-inducible promoter (Halloran et al., 2000) (hsp70l:Gq-Q209L). First, we determined whether heat shock induces Gq-Q209L expression in the craniofacial mesenchyme by examining changes to ERK1/2 (Mapk3/Mapk1) phosphorylation levels, an established readout for Gq/11 activity (Van Raamsdonk et al., 2009). Immunohistochemistry of coronal sections using an antibody for phosphorylated ERK1/2 (Randlett et al., 2015) confirmed that ERK1/2 phosphorylation levels were elevated in the craniofacial mesenchyme of hsp70l:Gq-Q209L embryos but not in sibling control embryos (Fig. S10). Next, embryos generated from edn1^+/− and edn1^+/−;hsp70l:Gq-Q209L mating pairs were heat-shocked at 16 hpf and then processed for fluorescence in situ hybridization at 28 hpf or skeletal preparations at 4 dpf. We performed skeletal preparations at 4 dpf instead of 6 dpf because a small percentage of hsp70l:Gq-Q209L-positive larvae exhibited early embryonic lethality starting at 5 dpf due to substantial cardiac edema (Fig. S11). In edn1^+/+ or edn1^−/− larvae (sibling controls), heat shock did not affect development of the craniofacial skeleton (Fig. 6A,B). In heat-shocked edn1^−/−;hsp70l:Gq-Q209L embryos, the majority of Edn1/Ednra-dependent structures were restored in 100% (9/9) of larvae (Fig. 6D,G). The palatoquadrate, symplectic, ceratohyal, Meckel's cartilage and the hyomandibular joint were restored at the highest frequency, whereas jaw joint restoration occurred at a lower frequency (Fig. 6G). We also observed homeotic transformations of upper jaw structures to lower jaw-like structures in heat-shocked edn1^+/+;hsp70l:Gq-Q209L (Fig. 6C,F,H) and edn1^−/−;hsp70l:Gq-Q209L (Fig. 6D,H) embryos. The pterygoid process of the palatoquadrate, typically a short, pointed structure (Fig. 6A,B,E), was transformed into a long, broad piece of cartilage resembling Meckel's cartilage (Fig. 6C,D,F). The dorsal portion of the hyomandibular cartilage was also hypoplastic and misshapen (Fig. 6,D,F). Similar transformations to the pterygoid process and hyomandibular cartilages have been reported in larvae with ectopic Ednra signaling activity in the dorsal domain of the first and second arches (Kimmel et al., 2007; Alexander et al., 2011; Zuniga et al., 2011; Barske et al., 2018).

These morphological changes were preceded by changes to patterning gene expression (Fig. 7). In 28 hpf edn1^−/− embryos, the expression of dlx5a, a transcription factor essential for lower jaw specification, was severely reduced in the first and second pharyngeal arch mesenchyme (Fig. 7B,B′). In edn1^−/−;hsp70l:Gq-Q209L embryos that were heat-shocked, dlx5a expression was not only restored in the ventral and intermediate domains of the first and second arches, but also expanded into the dorsal domain (Fig. 7F,F′). Similar expansion of dlx5a was observed in edn1^+/+;hsp70l:Gq-Q209L embryos (Fig. 7E,E′). Changes were also observed for nr2f5, an essential upper jaw specification factor (Barske et al., 2018) (Fig. 7C-D′,G-H′). In edn1^+/+ embryos, nr2f5 expression was confined to the dorsal and intermediate domains of the first and second arches (Fig. 7C,C′), whereas expression expanded into the ventral domains of the first and second arches in edn1^−/− embryos (Fig. 7D,D′). In edn1^+/+;hsp70l:Gq-Q209L and edn1^−/−;hsp70l:Gq-Q209L embryos, however, nr2f5 was reduced overall (Fig. 7G-H′). These data indicate that Gq is sufficient for induction of lower jaw patterning genes and subsequent development of lower jaw structures.

In this study, we addressed a long-standing question regarding the role of the Gq/11 family and potentially other Gα family members in lower jaw patterning (Sato et al., 2008a). Using the zebrafish model, we specifically interrogated the role of Gq/11 with three orthogonal approaches. We showed that inhibition of Gq/11 activity with a small molecule inhibitor or mutant alleles of gna11a and gna11b produced craniofacial phenotypes similar to those observed in edn1^−/− larvae and larvae injected with morpholinos targeting ednraa and ednrab (Miller et al., 2000; Nair et al., 2007). We also showed that transgene-mediated expression of a constitutively active Gq protein, Gq-Q209L, can rescue lower jaw structures in edn1^−/− embryos and can also transform the pterygoid process of the palatoquadrate to a Meckel's cartilage-like structure. Together, these results suggest that Gq/11signaling activity is necessary and sufficient to activate gene expression programs in cranial NCCs that establish lower jaw identity in zebrafish.

Our genetic analysis with mutant alleles has also revealed differential requirements for Gq/11 genes in zebrafish craniofacial development. Of the four members of the Gq/11 gene family – Gq, G11, G14 and G15 – we generated mutant alleles for Gq (gnaq) and G11 (gna11a and gna11b) only. Larvae that were homozygous mutants for both G11 genes (gnaq^+/+;gna11a^−/−;gna11b^−/−) recapitulated the edn1^−/− mutant phenotype (Figs 5C and 2G), suggesting that G11 genes are required for lower jaw patterning, whereas the Gq, G14 and G15 genes are not. The dispensability of G14 and G15 in lower jaw patterning is consistent with studies that have reported that expression of G14 and G15 genes is restricted to sensory organs and hematopoietic lineages (Ohmoto et al., 2011; Oka and Korsching, 2011). The dispensability of Gq was surprising, however, given its role in mouse lower jaw development (discussed further below). One of the most surprising findings from this work was the phenotypic discrepancy between zebrafish and mice lacking Gq and G11. The phenotype of gna11a^−/−;gna11b^−/− mutant zebrafish larvae are identical to edn1 mutants (Miller et al., 2000), with both exhibiting severe hypoplasia of the Meckel's cartilage, fusion of the jaw joint and loss of the symplectic and ceratohyal cartilages. In contrast, craniofacial phenotypes in mice do not overlap between Edn1/Ednra pathway mutants and Gq/11 mutants. In embryonic day (E) 18.5 Edn1^−/−, Ednra^−/− and Ece1^−/− mouse embryos, the entire mandible undergoes a homeotic transformation into a maxilla-like structure (Kurihara et al., 1995; Clouthier et al., 1998; Yanagisawa et al., 1998b). However, in E18.5 Gnaq^flox/flox;Gna11^−/−;P0-Cre mouse embryos, lower jaw defects were limited to the proximal two-thirds of the mandible (Dettlaff-Swiercz et al., 2005) and more closely resembled zebrafish embryos treated with YM (Figs 2 and 3) or reduced Gq/11 gene dosage (Fig. 5) than Edn1^−/− or Ednra^−/− mouse embryos (Clouthier et al., 1998; Ozeki et al., 2004). One explanation for the partial lower jaw defect in Gnaq^flox/flox;Gna11^−/−;P0-Cre embryos is that gene recombination in first pharyngeal arch NCCs induced by P0-Cre is less efficient compared to the NCC deletion strain Wnt1-Cre (Chen et al., 2017). Therefore, the residual levels of Gq expression and activity in cranial NCCs of Gnaq^flox/flox;Gna11^−/−;P0-Cre embryos may be sufficient to pattern the ventral domain but not the intermediate domain. This is consistent with numerous studies, including this one, that have reported that cranial NCCs in the intermediate domain are highly sensitive to perturbations that diminish Ednra signaling activity, such as treatment with YM (Figs 2 and 3) or an Ednra antagonist (Ruest and Clouthier, 2009), reduction of Gq/11 gene dosage (Fig. 5) (Offermanns et al., 1998), or reduction of Edn1 levels with an Edn1 morpholino (Miller and Kimmel, 2001) or a mutant allele of furina (Walker et al., 2006), an enzyme in the Edn1 biosynthetic pathway. This sensitivity of cranial NCCs in the intermediate domain has been attributed to a morphogen gradient model for Edn1, which posits that extracellular concentrations of Edn1 are highest near the source, the ventral ectoderm, and diminishes towards the dorsal regions of the arch (Kimmel et al., 2003). This question regarding the role of Gq/11 in mouse lower jaw development could be resolved with Wnt1-Cre, which has been shown to induce sufficient levels of gene recombination in cranial NCCs in Ednra^flox/flox embryos (Ednra^flox/flox;Wnt1-Cre) to recapitulate the Ednra^−/− phenotype (Ruest and Clouthier, 2009).

The zebrafish mutant alleles for Gq/11 genes revealed two additional differences between zebrafish and mice. First, mice and zebrafish embryos have different genetic requirements for Gq and G11 genes. In mice with conventional Gq/11 knockout alleles, embryos with one wild-type allele for Gna11 (Gnaq^−/−;Gna11^+/−) exhibited lower jaw defects, whereas embryos with one wild-type allele for Gnaq (Gnaq^+/−;Gna11^−/−) did not, indicating that Gnaq is primarily responsible for driving lower jaw patterning (Offermanns et al., 1998). In zebrafish, gna11a and gna11b were required for lower jaw patterning, but gnaq was dispensable (Fig. 5). One possible explanation for this discrepancy is that gnaq does in fact contribute to lower jaw patterning in zebrafish, albeit to a lesser extent, with the loss of gnaq resulting in skeletal defects that are imperceptible with our current analysis. An alternative explanation is that gnaq is simply not expressed in the lower jaw progenitor cells and therefore does not contribute to lower jaw patterning in zebrafish. The latter explanation is consistent with our gene expression analyses using scRNA-seq (Fig. S6), qPCR (Fig. S7) and in situ hybridization (Fig. S8), which suggest that gnaq expression is absent or expressed at lower levels relative to gna11a and gna11b in cranial NCCs. This difference may reflect divergent regulation of Gq/11 genes between teleosts and mammals. The teleost genome, which has undergone an additional round of duplication compared with mammals, exhibits significant divergence in the cis-regulatory elements (CREs) of duplicated genes and paralogs that result in different spatial and temporal expression patterns compared with other vertebrate species (Kassahn et al., 2009). The CREs for Gq/11 genes in teleosts may have also undergone changes that result in expression patterns that differ from mammals. A comparative analysis of the CREs and expression profiles for Gq/11 genes in zebrafish and mice may provide insights into evolutionary mechanisms that drive changes to gene regulatory networks for cell signaling proteins, which could have implications for understanding how the form and function of anatomical structures have diverged between teleosts and mammals (Pires-daSilva and Sommer, 2003; de Mendoza et al., 2014).

The second difference was that zebrafish and mice are differentially sensitive to loss of Gq/11 signaling. Gnaq^−/−;Gna11^−/− mouse embryos exhibit mid-gestational lethality due, in part, to myocardial hypoplasia (Offermanns et al., 1998). However, zebrafish embryos that are triple homozygous for Gq/11 knockout alleles did not exhibit evidence of embryonic lethality prior to 6 dpf based on general morphology and normal Mendelian inheritance of alleles (Table S3). Although heart development was not examined in this study, it is possible that our triple homozygous Gq/11 mutant zebrafish have similar heart defects as mice given the conservation of Gq/11 function in craniofacial development. Unlike mouse embryos, zebrafish larvae can tolerate defective heart function up to 7 dpf due to a lack of dependence on oxygen from the circulatory system (Stainier et al., 1996). Therefore, this property of zebrafish physiology could be leveraged to examine the role of Gq/11 in heart development, which is not fully understood (Malbon, 2005; Wettschureck and Offermanns, 2005).

Our study also provides additional evidence for the utility and specificity of YM as a pharmacological tool to interrogate the role of Gq/11 in vivo (Shibata et al., 2022). Embryos treated with 100 μM YM (the highest concentration that we could use; see Materials and Methods) produced phenotype and gene expression changes similar to embryos with partial loss of Edn1/Ednra signaling (Figs 2 and 3) (Miller and Kimmel, 2001; Walker et al., 2006; Nair et al., 2007; Ruest and Clouthier, 2009). The window of sensitivity to YM also agrees with the window of sufficiency for endothelin signaling in lower jaw patterning from previous studies (Miller et al., 2000, 2003; Walker et al., 2006; Kimmel et al., 2007; Ruest and Clouthier, 2009; Vieux-Rochas et al., 2010; Alexander et al., 2011; Zuniga et al., 2011; Barske et al., 2016, 2018; Meinecke et al., 2018). Further, applying YM on edn1 heterozygous mutants produced craniofacial phenotypes that were more severe in the spectrum of Edn1/Ednra-class phenotypes (Fig. 4). Although we acknowledge that the use of YM in zebrafish studies is limited by the concentration we can use, this additive effect of YM on sensitized backgrounds could be advantageous for studying other Gq/11-coupled receptors that might otherwise be embryonic lethal when knocked out. The effect of YM on cranial NCCs in zebrafish embryos is in stark contrast to the potency and efficacy of YM and its sister compound, FR900359, in cell culture assays, which exhibit an IC50 in the nanomolar concentration range (Nishimura et al., 2010; Schrage et al., 2015). The diminished inhibitory effect of YM in zebrafish embryos may be caused by poor perfusion into relevant tissues, inefficient cell permeability, or degradation. FR900359 exhibits slightly different physicochemical properties (Schlegel et al., 2021) and thus may exhibit greater effects on zebrafish embryos than YM.

By establishing Gq/11 as the sole signaling mediator for Ednra in lower jaw patterning, we can begin to elucidate mechanisms of jaw patterning that remain poorly characterized. As demonstrated by our study and others, correct positioning of the upper and lower jaw patterning domains requires a specific balance of signaling inputs by Edn1/Ednra, Jagged/Notch and Nr2f nuclear receptors (Sato et al., 2008b; Zuniga et al., 2010; Alexander et al., 2011; Tavares and Clouthier, 2015; Barske et al., 2016, 2018; Sucharov et al., 2019; Kurihara et al., 2023). It remains unknown how these signaling pathways are integrated to elicit specific transcriptional outputs that refine the positional identities of cranial NCCs residing in and around the upper jaw–lower jaw boundary. Future studies will examine mechanisms of signal integration.

Protein sequences were aligned using Clustal Omega (Madeira et al., 2019). The following are accession numbers for sequences used. For human proteins: Gq (NP_002063.2), G11 (NP_002058.2) and Gi1 (NP_002060.4). For zebrafish proteins: Gq (NP_001138271.1), G11a (NP_001038501.1), G11b (NP_001007774.1), G14 (NP_001003753.1), G14a (XP_683989.2), G15.1 (NP_001003626.2), G15.2 (XP_002667410.2), G15.3 (translated from mRNA sequence for gna15.3, XR_659583.3), G15.4 (NP_001038454.1), Gi1 (NP_957265.1), Gs (XP_001335732.1) and G12a (NP_001013295.1).

All work was approved by the University of Colorado Institutional Animal Care and Use Committee (Protocol No. 00188). Zebrafish (Danio rerio) adults and embryos were raised and staged according to Kimmel et al. (1995) and Westerfield (2007). Mutant and transgenic lines were maintained in the AB strain. The following mutant and transgenic lines were previously described: sucker/edn1^tf216 (Miller et al., 2000), Tg(fli1a:EGFP)^γ1 (Lawson and Weinstein, 2002) and Tg(sox10:mRFP)^vu234 (Kirby et al., 2006). In addition, five new lines were created: mutant alleles for edn1^co3009, gnaq^co3010, gna11a^co3015 and gna11b^co3014, and the transgenic line Tg(hsp70l:Gq-Q209L-IRES-EGFP,cmlc2:EGFP). All lines are available upon request. Construction of strains is described below.

All mutant alleles were generated with the CRISPR/Cas9 gene-editing system as previously described (Hwang et al., 2013; Jao et al., 2013). Gene-specific target sequences were generated with CHOPCHOP v.2 (Labun et al., 2016), and guide RNAs (gRNAs) were generated as previously described (Bassett et al., 2013). To generate edn1^co3009, two target sequences were identified in exon 2. To generate gnaq^co3010 and gna11b^co3014, target sequences were identified in exon 4 (Fig. S9). To generate gna11a^co3015, a target sequence was identified in exon 3 (Fig. S9). The following are the target sequences used for each gene: edn1, 5′-GGAATAAGAGATGCTCCTGC-3′and 5′-GGACATAATATGGGTGAACA-3′; gnaq, 5′-GGCTGGGTGGGAATGTAGGAA-3′; gna11a, 5′-GGCGAGAGGTCGATGTCGAGA-3′; gna11b, 5′GGTGGGAAGGTACGAAGATT-3′. The individual target sequences were then incorporated into an oligo consisting of 5′-[T7 promoter sequence]-[Target sequence]-[Start of scaffold sequence]-3′. The T7 promoter sequence is 5′-AATTAATACGACTCACTATA-3′, and the start of the scaffold sequence is 5′-GTTTTAGAGCTAGAAATAGC-3′. DNA templates for gene-specific gRNAs were then generated by performing PCR with the target sequence-containing oligo and a separate oligo containing the gRNA scaffold sequence: 5′-GATCCGCACCGACTCGGTGCCACTTTTTCAAGTTGATAACGGACTAGCCTTATTTTAACTTGCTATTTCTAGCTCTAAAAC-3′. All oligos were purchased from IDT DNA Technologies. PCR reactions were performed with KAPA HiFi HotStart ReadyMix (Roche, 07958927001) (1× KAPA Taq, 500 nM target sequence oligo, 500 nM gRNA scaffold oligo) using the following cycling parameters: 98°C for 30 s; 40 cycles of 98°C for 10 s, 60°C for 10 s; 72°C for 15 s; 72°C for 10 min. Following purification of the PCR product using the QIAquick PCR Purification Kit (QIAGEN, 28104), the gRNA was synthesized using the MEGAscript T7 Transcription kit (Invitrogen, AM1334). Cas9 mRNA (codon-optimized for zebrafish) was generated as previously described (Jao et al., 2013). Briefly, pT3TS-nCas9n (Addgene plasmid #46757) was linearized with XbaI restriction enzyme and then used as the template for in vitro transcription using the mMESSAGE mMACHINE T3 Transcription Kit (Invitrogen, AM1348).

Embryos were injected at the one-cell stage with 2 nl of solution containing 100 ng/µl Cas9 mRNA, 50 ng/µl gRNA and 0.025% Phenol Red (Sigma-Aldrich, P0290). To determine whether the gRNAs were excising the genomic loci of interest, at least ten embryos per injected clutch were lysed in genomic DNA extraction buffer (10 mM Tris pH 8, 2 mM EDTA, 0.2% Triton X-100, 0.2 mg/ml Proteinase K) at 24 hpf, and the genomic loci of interest were PCR-amplified with genotyping primers. PCR amplicons were then analyzed by gel electrophoresis, with the presence of smears or ladders serving as indicators of potential genomic lesions. Clutches indicating genomic lesions were grown to adults (F0) and then crossed with wild-type ABs to generate F1 embryos. F1 adults were screened for genomic lesions at loci of interest by PCR, and the nature of the genomic lesions were characterized by Sanger sequencing. F1 adults harboring frameshift mutations were crossed with wild-type ABs to identify those that transmit mutant alleles to offspring in the expected Mendelian frequency. Germ line-stable F1 adults were kept and outcrossed to wild-type ABs to generate F2 fish used in this study. The following genotyping primers were used for each allele: edn1^co3009, 5′-TAGGTGCTCCAGCATCTTTG-3′ and 5′-GGAGCGTTTCCAAGTCCATA-3′; gnaq^co3010, 5′-TATGATGCCTTTTTGTCCACAG-3′ and 5′-CATTGTCAGACTCGACGAGAA-3′; gna11a^co3015, 5′-TCAAGTCGTAAAATGGGTTGTG-3′ and 5′-TAAAAGCAGCAAATGACGACAC-3′; gna11b^co3014, 5′-GAAGCATCCTTTACCAAACCAC-3′ and 5′-CGAACGAAGGCAGATGAATAAT-3′. The frameshift mutation in the gna11b^co3014 allele introduces a BceAI restriction site. Thus, PCR reactions from the gna11b^co3014 genotyping assay were digested with BceAI to further distinguish PCR amplicons from wild-type and mutant alleles.

The hsp70l:Gq-Q209L-IRES-EGFP,cmlc2:EGFP plasmid was generated using the Tol2kit as described by Kwan et al. (2007). First, the middle entry vector containing Gq-Q209L, pME-Gq-Q209L, was generated. The Gq-Q209L cDNA was PCR-amplified from pc3.1-Gq-Q209L (Kanai et al., 2022) and then cloned into the pME-MCS vector using T4 ligase (New England Biolabs, M0202). pME-Gq-Q209L was then combined with p5E-hsp70l, p3E-IRES-GFPpA and pDestTol2CG2 using LR Clonase II Plus Enzyme Mix (Invitrogen,12538-120). To generate transposase mRNA, the pCS2FA-transposase plasmid was linearized with NotI restriction enzyme and then used as the template for in vitro transcription using the mMESSAGE mMACHINE SP6 Transcription Kit (Invitrogen, AM1340).

Embryos at the one-cell stage were injected with a 2 nl solution containing 12.5 ng/µl transposase mRNA, 12.5 ng/µl hsp70l:Gq-Q209L-IRES-EGFP,cmlc2:EGFP and 0.025% Phenol Red. Injected embryos were grown to adults (F0) and then crossed to wild-type ABs. Offspring (F1) were screened for the presence of the transgene by examining GFP expression in the heart. Transgene-positive, germ line-stable F1 fish were then crossed into edn1^+/− fish to generate Tg(hsp70l:GqQ209L-IRES-EGFP); edn1^+/− fish. The plasmids pME-MCS, p5E-hsp70l, p3E-IRES-GFPpA, pDestTol2CG2 and pCS2FA-transposase were provided by Kristen Kwan (Kwan et al., 2007).

Embryos generated by intercrossing wild-type ABs were dechorionated at 14 hpf in Petri dishes coated with 0.5% agarose using E2 media (Westerfield, 2007) containing 2 mg/ml Pronase (Roche, 10165921001) at 25°C for 10-15 min. Embryos were subsequently washed in fresh E2 media in agarose-coated Petri dishes three times. Dechorionated embryos in agarose-coated plates were then incubated with DMSO or YM [Adipogen (AG-CN2-0509-M001) or Cayman Chemicals (29735)] dissolved in DMSO, from 16 to 36 hpf. At 36 hpf, embryos were washed with fresh E2 media three times and then processed for skeletal preparations at 6 dpf. Embryos for fluorescence in situ hybridization with immunofluorescence were fixed overnight at 36 hpf in 4% paraformaldehyde (PFA) in PBS. Note that the product datasheet for YM states that the solubility limit in DMSO is 10 mM (9.6 mg/ml), which is the stock concentration we used. However, for zebrafish embryo media, 1% DMSO is the upper limit that does not cause developmental defects (Hoyberghs et al., 2021). To ensure that we did not exceed 1% DMSO in the embryo media, 100 µM YM was the highest working concentration we could use.

Larvae were stained with Alizarin Red and Alcian Blue as previously described (Walker and Kimmel, 2007; Brooks and Nichols, 2017).

Cranial neural crest cells were isolated from zebrafish embryos using a fli1a:EGFP and sox10:mRFP double labeling technique as previously described (Barske et al., 2016; Askary et al., 2017; Mitchell et al., 2021) with some modifications. In brief, 140 double-positive embryos were dechorionated at 14 hpf with 2 mg/ml Pronase. Seventy embryos were then treated with 100 µM YM or DMSO from 16 to 36 hpf. At 36 hpf, ten embryos from each condition were set aside for skeletal preparations at 6 dpf. The heads of the remaining embryos were severed, pooled and de-yolked by gentle pipetting in Ca²⁺-free Ringer's solution (116 mM NaCl, 2.9 mM KCl, 5 mM HEPES, pH 7.0). The heads were then dissociated into single cells with mechanical agitation and enzymatic digestion for 15 min at 4°C with cold-activated protease from Bacillus licheniformis (Sigma-Aldrich, P5380) in dissociation buffer (10 mg/ml B. licheniformis protease, 125 U/ml DNAse, 2.5 mM EDTA, PBS), followed by enzyme neutralization with stop solution (30% fetal bovine serum, 0.8 mM CaCl₂, in PBS). Cells were then centrifuged (400 g for 5 min), resuspended in cell suspension buffer [1% fetal bovine serum, 0.8 mM CaCl_2, Leibovitz's L-15 Medium (Gibco, 21083-027)], filtered through a 70 µm strainer (PluriSelect, 43-10040-40), centrifuged again (400 g for 5 min) and resuspended in sorting buffer (1% fetal bovine serum, 1 mM EDTA, 25 mM HEPES). EGFP and mRFP double-positive cells were then purified by FACS using a MoFlo XDP100 (Beckman Coulter). A total of 72,000 and 86,000 double-positive cells were purified for DMSO and YM-treated samples, respectively. Cells were once again centrifuged (400 g for 5 min) and resuspended in cell suspension buffer. After confirming cell viability with Trypan Blue (Gibco, 15250-061), 8000 cells for each condition were submitted to the University of Colorado Anschutz Medical Campus Genomics and Microarray Core for scRNA-seq. Single-cell capture, barcoding and library generation was performed using 10x Genomics single cell RNA-sequencing technology (Chromium Controller and Next GEM single Cell 3′ Kit v.3.1, Dual Index). Libraries were processed for paired-end sequencing with a read depth of 75,000/cell using NovaSeq6000 (Illumina).

Sequence reads were aligned to the D. rerio reference genome (GRCz11) and converted to count matrices with Cell Ranger v.5.0.1 (10x Genomics). Subsequent analyses of count matrices were performed with Seurat v.4.3.0 in R studio (Hao et al., 2021). Individual Seurat objects were created for DMSO and YM-treated samples, and cells in each dataset were filtered for quality. For DMSO-treated samples, 3700 cells were obtained after filtering for cells with 200-3000 unique feature counts and less than 5% mitochondrial DNA. For YM-treated cells, 4506 cells were obtained after filtering for cells with 200-3500 unique feature counts and less than 5% mitochondrial DNA. Datasets for DMSO or YM-treated samples were then merged, normalized and scaled, with cell cycle phase-associated genes regressed out. A list of S- and G2M-phase genes built into Seurat was converted to zebrafish identifiers using BioMart and then used for cell cycle regression (Nestorowa et al., 2016). Converted gene lists were manually checked and edited for any changes in nomenclature/duplications between the reference and Ensembl database. Linear dimensional reduction was then performed with principal component analysis (PCA), and the top 30 PCAs were used to generate and visualize clusters with the uniform manifold approximation and projection (UMAP) nonlinear dimensional reduction technique. Marker genes for clusters were obtained with a Wilcoxon Rank Sum test, and cluster identities were determined using published, experimentally verified marker genes (Nair et al., 2007; Talbot et al., 2010; Askary et al., 2017; Barske et al., 2018, 2016; Fabian et al., 2022). The approximate population of cells representing the intermediate domain were determined with a feature plot encompassing known intermediate domain marker genes, ccn2b, dlx3b, dlx4a, dlx4b, emx2, fgfbp2a, foxc1b, foxd1, fsta, grem2b, igfbp5b, msx1a, nkx3-2 and shox (Fig. S4). Differential expression analysis was performed using Wilcoxon Rank Sum test between equivalent clusters of cells in DMSO and YM-treated samples.

RNA was collected from 40 heads of 28 hpf zebrafish embryos using Direct-zol RNA MiniPrep Kit (Zymo Research, R2050) following the manufacturer's protocol. cDNA was generated using the QuantiTect Reverse Transcription Kit (QIAGEN, 205311). qPCR assays were prepared using the QuantiTect SYBR Green PCR kit (QIAGEN, 204143) following the manufacturer's protocol and read in a CFX Connect Real-Time PCR Detection System (Bio-Rad). The following QuantiTect Primers were purchased from QIAGEN (249900): gnaq (GeneGlobe ID: QT02091502), gna11a (GeneGlobe ID: QT02065042), gna11b (GeneGlobe ID: QT02062060), actb1 (GeneGlobe ID: QT02174907). The primer pair for dlx2a was: 5′-gaaacgctttcggccccta-3′ and 5′-ccattcggatttcaggttcgc-3′ (Weinschutz Mendes et al., 2020). actb1 was used as the reference gene to calculate relative gene expression.

Colorimetric in situ hybridization was performed as previously described (Thisse and Thisse, 2008). Two-color fluorescence in situ hybridization was performed as previously described (Talbot et al., 2010). Fluorescence in situ hybridization with immunofluorescence was performed as previously described (Mitchell et al., 2021). The following probes have been previously described: dlx2a (Akimenko et al., 1994), dlx4b (Ellies et al., 1997), dlx5a (Walker et al., 2006) and nr2f5 (Barske et al., 2018). The plasmid for nr2f5 was provided by Lindsey Barske (Cincinnati Children's Hospital Medical Center, OH, USA). To generate probes for gnaq, gna11a and gna11b, gene fragments were first cloned into the pCRII plasmid with the TOPO-TA Cloning Kit (Invitrogen, 450640) following the manufacturer's protocol. The following primer pairs were used to amplify gene fragments: gnaq, 5′-cgttaacacgggaggaacac-3′ and 5′-gagctgtcggtcgatctcat-3′ (Takeuchi et al., 2017); gna11a, 5′-tactcgcacttcacctgtgc-3′ and 5′-tgggaaaggcgttttatttg-3′; and gna11b, 5′-caccgacacagagaacatcc-3′ and 5′-acattcatcgatgcgagttg-3′.

Embryos processed for whole-mount fluorescence in situ hybridization and immunofluorescence were mounted in 0.4% agarose in PBS on a glass-bottom dish (MatTek, P35G-1.5-10-C), with images taken using a 20× air objective on a DMi8 microscope (Leica) equipped with an Andor Dragonfly confocal unit (Oxford Instruments). Images were subsequently processed using Imaris image analysis software (Oxford Instruments). Embryos processed for whole-mount colorimetric in situ hybridization and skeletal preparations were imaged on a SZX12 stereo microscope (Olympus) equipped with a SC100 camera (Olympus). Sections for colorimetric in situ hybridization were imaged using a 20× objective on a BX51 microscope (Olympus) under Nomarski optics, with images captured using a DP71 camera (Olympus). Sections processed for immunohistochemistry were imaged with a 63× oil immersion objective on a SP8 confocal microscope (Leica), and image processing was subsequently performed in Fiji (Schindelin et al., 2012).

Embryos generated from crossing Tg(hsp70l:GqQ209L-IRES-GFP,cmlc2:EGFP);edn1^+/− and edn1^+/− adults were heat-shocked at 16 hpf for 10 min in a 38°C water bath and then returned to 28.5°C. Transgenic and non-transgenic embryos were separated based on cardiac EGFP expression and monitored daily until 4 dpf to account for embryos developing cardiac edema, and then processed for skeletal preparations. Embryos processed for in situ hybridization were fixed overnight in 4% PFA in PBS at 28 hpf.

Zebrafish embryos collected from hsp70l:Gq-Q209L and sox10:mRFP mating pairs were heat-shocked at 21 hpf and fixed 3 h later at 24 hpf in 4% PFA overnight at 4°C. Embryos were washed in PBS for 5 min three times, successively incubated in 15% sucrose/PBS and 30% sucrose/PBS, embedded in Tissue-Tek O.C.T. Compound (Sakura, 4583), and flash-frozen in an ethanol/dry ice bath. Embryos were then sectioned at 20 μm on a CM1900 cryostat (Leica). Sections were hydrated in PBS for 5 min, permeabilized in 0.1% Triton X-100/PBS (PBS-Triton) for 5 min, and then incubated in blocking buffer (2% bovine serum albumin, 2% goat serum, 0.1% Triton X-100, PBS) for 2 h at room temperature. Sections were then incubated in blocking buffer containing 1:200 dilution of primary antibody [rabbit P-p44/42MAPK (T202/Y204), Cell Signaling Technology, 4370, lot #17)] overnight at 4°C. Sections were washed for 15 min in PBS-Triton three times and then incubated in blocking buffer containing 1:500 dilution of secondary antibody (goat anti-rabbit Alexa 647, Invitrogen, A21245) for 2 h at room temperature. Sections were then washed in PBS-Triton for 15 min×3, incubated in PBS containing 0.1 µg/ml Hoechst 33342 (Invitrogen, H1399) for 15 min, washed in PBS for 5 min, and mounted in Prolong Gold Antifade Reagent (Invitrogen, P36930).

Statistical analyses for the one-way ANOVA with Holm–Šídák multiple comparisons test, chi-square tests and unpaired t-tests were performed in Prism (GraphPad).

We thank Christine Archer, Ainsley Gilbard and Olivia Gomez for zebrafish care; the Genomics and Flow Cytometry Shared Resources at the University of Colorado Cancer Center (both funded by the Cander Center Support Grant P30CA046934 from the National Cancer Institute) for assistance with sample preparation (single cell capture, library preparation and FACS) for scRNA sequencing; members of the Nichols lab, Dr Joan Hooper and the craniofacial biology interest group (‘11 at 11’ joint lab meetings) for helpful comments and feedback.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.204396.reviewer-comments.pdf