During craniofacial development, different populations of cartilage- and bone-forming cells develop in precise locations in the head. Most of these cells are derived from pluripotent cranial neural crest cells and differentiate with distinct developmental timing and cellular morphologies. The mechanisms that divide neural crest cells into discrete populations are not fully understood. Here, we use single-cell RNA sequencing to transcriptomically define different populations of cranial neural crest cells. We discovered that the gene family encoding the Alx transcription factors is enriched in the frontonasal population of neural crest cells. Genetic mutant analyses indicate that alx3 functions to regulate the distinct differentiation timing and cellular morphologies among frontonasal neural crest cell subpopulations. This study furthers our understanding of how genes controlling developmental timing shape craniofacial skeletal elements.
The vertebrate head skeleton has many cartilaginous and bony elements. Craniofacial skeletal patterning, or how skeletal cells reproducibly form these different structures with their appropriate shapes and positioning, has been a subject of intense study for many years (Creuzet et al., 2005; Olsson et al., 2005; Graham, 2003). The mechanisms patterning the lower jaw are the most well understood (Medeiros and Crump, 2012; Schilling and Kimmel, 1997). In contrast, relatively little is known about the mechanisms patterning the anterior neurocranium, which serves as the skull base. The anterior neurocranium consists of the skeletal elements anterior to the junction of the trabeculae and the basal plate (‘polar cartilage’ region) (Eberhart et al., 2006; Cubbage and Mabee, 1996; Schilling and Kimmel, 1997). In larval zebrafish, the individual neurocranium skeletal elements are the edge and medial portions of the ethmoid plate, the paired trabecular cartilages, and the dermal parasphenoid bone.
Craniofacial skeletal development occurs in several phases in all vertebrates. In the first phase, the neural crest cells (NCCs) that form most of the facial skeleton migrate away from the dorsal neural tube into three locations in the embryo: the frontonasal region, the anterior pharyngeal arches (arches one and two) and the posterior pharyngeal arches (arches three to seven in zebrafish). In zebrafish, these three groups of post-migratory NCCs undergo complex morphogenetic movements before giving rise to the skeleton of the anterior neurocranium, the jaw, and the gill supports, respectively (Crump et al., 2006; Schilling and Kimmel, 1994, 1997; Birkholz et al., 2009; Swartz et al., 2011; Wada et al., 2005; Szabo-Rogers et al., 2010; Sasaki et al., 2013; Talbot et al., 2016). During patterning, post-migratory NCCs in each region are subdivided into distinct transcriptional programs and adopt different cellular identities dependent upon their position. These specific transcriptional programs determine the future cellular fate of NCCs. As a result, NCCs differentiate into chondrocytes and osteoblasts, the cells that produce cartilage and bone, with region-specific developmental timing. Wild-type development reproducibly results in correctly positioned cell populations with different transcriptional identities giving rise to the individual skeletal elements. Studying genetic mutants in which specific skeletal elements are affected deepens our understanding of the mechanisms of craniofacial development (Mork and Crump, 2015).
Transcription factor-encoding genes that are enriched in certain NCC populations have been shown to be important for patterning. For example, the distal-less (Dlx) gene family is enriched in pharyngeal arch NCCs (Talbot et al., 2010). Specific Dlx gene expression combinations specify an identity code that divides the pharyngeal arch NCCs into subpopulations (Depew et al., 2005). Accordingly, mutations in vertebrate Dlx genes can result in homeosis, or the transformation of one structure into another. For example, ventrally residing cells in mouse Dlx5/6 mutants are homeotically transformed to reflect a dorsal identity. The altered identity of these cells results in skeletal element shape changes including duplicated dorsal skeletal elements (Depew et al., 2002; Beverdam et al., 2002). Although the Dlx code is well described, an identity code regionalizing frontonasal NCCs remains undefined.
Recent work using RNA sequencing has advanced our understanding of the regionalization of the zebrafish facial skeleton (Askary et al., 2017; Barske et al., 2016), uncovering pharyngeal arch NCC subpopulations with different transcriptional identities. However, frontonasal NCC specific transcriptional programs have not yet been examined with transcriptomic analyses. Single-cell RNA sequencing (scRNA-seq) technology is ideal for understanding the precise regionalization of the heterogeneous NCC populations, including the anterior neurocranium progenitors. This technique has not yet been reported for zebrafish NCCs.
Several lines of evidence suggest that the frontonasal NCCs are regionalized into distinct subpopulations. Zebrafish frontonasal NCCs give rise to differently shaped anterior neurocranium skeletal elements, including the edge and medial ethmoid plate, the trabecular cartilages and the parasphenoid bone. Additionally, different populations of skeletal cells can have distinct cellular morphologies, even within a given skeletal element. For example, cartilage cells at the lateral edges of the ethmoid plate are more columnar, whereas the chondrocytes of the medial ethmoid plate are more cuboidal. Interestingly, Alcian Blue staining of cartilage reveals that the ethmoid plate edge chondrocytes have comparably stronger Alcian Blue staining in the surrounding extracellular matrix than the medial ethmoid plate chondrocytes. This difference in staining may be due to distinctive chondrocyte differentiation timing between different populations of frontonasal NCCs. Consistent with this idea, the chondrocytes at the edges of the ethmoid plate differentiate early in craniofacial development, whereas the NCCs that become medial ethmoid plate chondrocytes differentiate relatively late (Eberhart et al., 2006; Swartz et al., 2011; Wada et al., 2005; Dougherty et al., 2013). These findings indicate that the frontonasal NCCs are divided into subpopulations before giving rise to discrete skeletal elements. However, the gene code subdividing the frontonasal NCCs into different identities is unknown.
The aristaless-related transcription factor-encoding genes belong to the paired class of homeodomain protein-encoding genes. These proteins are characterized by a DNA-binding homeodomain, an aristaless (OAR) domain at the carboxyl terminus, and a mostly divergent N terminus (Galliot et al., 1999; Pérez-Villamil et al., 2004). The Alx genes are aristaless-related genes that are important for multiple developmental processes in diverse vertebrates. Human mutations in ALX genes are associated with craniofacial disorders that are broadly classified as frontonasal dysplasia (FND). There are three classes of FND: FND1, FND2 and FND3, which are associated with mutations in ALX3, ALX4 and ALX1, respectively (Uz et al., 2010; Kayserili et al., 2009; Twigg et al., 2009; Ullah et al., 2018). Although most ALX-associated human craniofacial defects are recessive conditions, some abnormalities have been reported in ALX4 heterozygous humans (Farlie et al., 2016). All FND1-associated ALX3 sequence variations are located in the homeobox. Skeletal dysmorphologies in these individuals include short nasal bridges, bifid nasal tips, broad columella, widely separated nares and dysmorphic ethmoid bone (Twigg et al., 2009; Ullah et al., 2018). The etiology of these ALX3-associated disease phenotypes is not fully understood.
Genetic studies in mice demonstrate that the Alx genes are important for mouse skeletal development, including formation of limb and craniofacial structures. Strong's Luxoid mouse is a spontaneously occurring semidominant mutation in Alx4 resulting in polydactyly. These mice show evidence of limb mispatterning indicative of identity transformation (Takahashi et al., 1998; Qu et al., 1997). Engineered mutations in mouse Alx1, Alx3 and Alx4 all exhibit craniofacial defects (Qu et al., 1998; Beverdam et al., 2001; Kayserili et al., 2009; Lakhwani et al., 2010), and mouse Alx3;Alx4 double mutants have severe cleft face phenotypes attributed to increased apoptosis (Beverdam et al., 2001). Yet, a role for these genes in NCC patterning has not been proposed, and the precise function of Alx genes during craniofacial development remains vague. In particular, Alx3 mutant mice were initially reported to have no craniofacial phenotype (Beverdam et al., 2001), but later reports with the same allele indicate incompletely penetrant neural tube closure defects and clefting (Lakhwani et al., 2010). These severe defects are not observed in human FND1 patients with ALX3 mutations, nor have the human disease phenotypes characteristic of FND1 been found in Alx3 mutant mice (Farlie et al., 2016). Thus, further study in additional animal models is required to resolve the function of Alx3 in vertebrate craniofacial development.
Zebrafish are a powerful genetic model for understanding craniofacial development, allowing cellular analyses of NCC migration, morphogenetic movement, patterning and differentiation. Although there are four known zebrafish Alx genes: alx1, alx3, alx4a and alx4b – genetic mutant analysis has not been reported for any of them. In a published series of experiments with morpholino knockdown of the Alx family, only morpholinos targeted to alx1 produced craniofacial dysmorphologies in zebrafish. This work reported that alx1 morpholinos disrupt early NCC migration (Dee et al., 2013). However, morpholino injection phenotypes are not always concordant with germline genetic loss of function phenotypes (Kok et al., 2015). Sometimes this discordance is due to genetic compensation and transcriptional adaptation (El-Brolosy et al., 2019; Rossi et al., 2015; Ma et al., 2019). But in other cases, morpholino injection can produce non-specific defects unrelated to the morpholino-targeted gene (Tessadori et al., 2020). Thus, stable germline mutations remain important for understanding gene function in zebrafish. Other expression studies indicate that the Alx genes may pattern identity in the zebrafish pectoral fin (Nachtrab et al., 2013). However, a role for this gene family in patterning skeletal identity in the head has not yet been proposed.
Here, we present scRNA-seq data defining discrete NCC populations at 24 h post-fertilization (hpf). The frontonasal, anterior pharyngeal arch and posterior pharyngeal arch NCC populations each have unique transcriptomic identities, revealing gene expression programs that are specific to each population. We find that the Alx gene family is strongly and specifically expressed in the frontonasal NCCs at this stage; the anterior and posterior pharyngeal arches are essentially devoid of Alx gene family expression. We confirm the transcriptomic expression data with in situ gene expression studies at cellular resolution in intact embryos. We used CRISPR/Cas9 to develop stable germline mutants for all members of the zebrafish Alx family. Surprisingly, of the mutations we generated, only those in alx3 produce craniofacial phenotypes. Careful mutant analyses reveal that zebrafish alx3 mutants develop craniofacial phenotypes restricted to the anterior neurocranium. Quantitative cell tracking demonstrates that alx3 function is dispensable for NCC migration and morphogenetic movement. We also demonstrate that alx3 function is not required for NCC survival or proliferation. Rather, our data demonstrate that alx3 is required for discrete differentiation timing among frontonasal NCC cell populations, as well as different cell morphology between populations of skeletal cells. These findings motivate a model in which alx3 functions to pattern or subdivide the zebrafish anterior neurocranium into distinct populations of cells, conferring specific shapes to neurocranium skeletal elements.
Single-cell RNA sequencing identifies four subpopulations of post-migratory neural crest cells
To profile NCCs transcriptomically, we isolated cells double positive for the transgenes sox10:mRFP and fli1a:EGFP by fluorescence-activated cell sorting (FACS) from 24 hpf wild-type embryos (Fig. 1A,B). Sorting cells that are double positive for these transgenes allows exquisitely specific purification of NCCs for downstream transcriptomic profiling (Askary et al., 2017; Barske et al., 2016). By 24 hpf, NCCs have migrated to their destination and region-specific transcriptomic programs have begun to emerge. Isolated NCCs were subjected to scRNA-seq analyses and bioinformatic clustering in Seurat. This processing resolved individual NCCs into four distinct populations based on their transcriptomic profiles (Fig. 1C). Our analyses produced lists of marker genes that are enriched in each population (Fig. 1D-F). We used previously published in situ hybridization gene expression and mutant phenotype data to infer the identity of the different clusters. For example, the Dlx genes are strongly represented in the largest cluster (46% of the cells) with three of the top 14 marker genes in the Dlx family (Fig. 1D). Many studies report that Dlx genes are strongly expressed in zebrafish pharyngeal arches one and two as well as more weakly expressed in posterior arch NCCs at this stage (Akimenko et al., 1994). Furthermore, loss of dlx5a affects skeletal structures derived from pharyngeal arches one and two (Talbot et al., 2010; Sucharov et al., 2019). This same cluster was also enriched for barx1, a gene with expression and mutant phenotypes associated with arches one and two (Nichols et al., 2013). These findings suggest that the largest cluster in our analyses contains cells from pharyngeal arches one and two. We identified a cluster with 13% of the cells to be the posterior arches; four of the 14 top marker genes are Hox genes (Fig. 1E). The Hox gene family codes for identity in the anteroposterior axis and these members of the Hox code are all restricted to the posterior arches (Prince et al., 1998; Vaccari et al., 2010; Laue et al., 2008a; Chen et al., 2010; van der Velden et al., 2013). This same cluster was also enriched for prdm1a, a gene with both expression and mutant phenotypes associated with the posterior arches and their skeletal derivatives (Ding et al., 2013; Birkholz et al., 2009). Additionally, the posterior arches three to seven are the last to develop and the least mature, which may explain why this population is relatively small at this stage compared with other populations of NCCs. The third major cluster contains 40% of the cells and was more difficult to characterize; many of the markers were not previously associated with a specific subset of NCCs (Fig. 1F). cyp26b1 was the top marker gene for this population and expression and function of this gene is associated with midline anterior neurocranium cartilages derived from frontonasal NCCs (McGurk et al., 2017; Laue et al., 2008b; Piotrowski et al., 1996). We propose that this population represents frontonasal NCCs, the progenitors of the anterior neurocranium skeleton. We also detected a small (1%) population of cells that was made up of three separate ‘satellite’ clouds that were more similar to each other than to any of the major clusters. One of these clouds was enriched for markers associated with the melanocyte lineage, including pcdh10a (Williams et al., 2018) and mitfa (Lister et al., 1999) (Fig. S1). Classical studies demonstrate that cranial NCCs that become melanocytes are lineage restricted early in development (Schilling and Kimmel, 1994). Therefore, we propose that some of the NCCs in this small population are already committed to become melanoblasts at this early stage (Table S1, Fig. S1).
Zebrafish NCCs are specified and initiate migration in a stereotypical anterior to posterior fashion (Raible et al., 1992; Schilling and Kimmel, 1994). We hypothesized that the most anterior, earliest specified NCCs might be more developmentally advanced than the populations that form later. To test this hypothesis, we performed pseudotemporal analyses using Monocle 3 (Trapnell et al., 2014). We found that the predicted developmental age from these analyses matched the stereotypical anterior-to-posterior development of NCCs (Fig. S2); the frontonasal cluster is the most developmentally advanced population followed by the anterior arches, then the late-arriving posterior arches are the most immature. These analyses confirm both our assignment of cluster identity, and the anterior-to-posterior timing of NCC development and indicate that pseudotemporal analysis might be useful for understanding cell state in our system. All individual gene expression feature plots from our Seurat and Monocle 3 processing are available in two easily searchable interactive web interfaces (https://apulvino.shinyapps.io/sAlxMapp/, https://apulvino.shinyapps.io/mAlxMapp/).
Alx gene family expression is highly enriched in zebrafish frontonasal NCCs
The putative frontonasal population was the most difficult to identify due to many of the marker genes not being well studied with regards to zebrafish craniofacial development. To better understand the genes that might specifically function in these cells, we closely examined the genes enriched in these cells compared with all the other cells in our dataset. Analysis of gene markers for this cluster reflected a disproportionately higher enrichment of Alx gene family transcripts in the frontonasal cluster. alx3 and alx4a are particularly enriched in the frontonasal cluster; both genes were represented in the top five frontonasal cluster markers. Further examination of this gene family in our data set revealed that all four zebrafish Alx genes were specifically expressed in this population and largely excluded from the other clusters (Fig. 2A). In the frontonasal NCCs, alx4a and alx3 were the most strongly expressed Alx genes followed by alx1 then alx4b (Fig. 2B). It is worth highlighting that mis-annotation of the alx3 locus in the zebrafish genome (GRCz11) presented a challenge in our analyses, see Materials and Methods and Fig. S3 for details.
To confirm that the alx3 gene is strongly and specifically expressed in NCCs residing in the frontonasal domain, we performed whole-mount fluorescent in situ hybridization analysis on 24 hpf embryos. We used the fli1a:GFP transgene employed in our FACS (Fig. 1) to label post-migratory NCCs. The strongest alx3 expression was in fli1a:EGFP-positive NCCs in the frontonasal domain (Fig. 2C-E). There was no detectable signal in anterior or posterior pharyngeal arch cells (Fig. S4). To further confirm that cells in this putative frontonasal NCC cluster reside in the frontonasal region of intact embryos, we examined alx4a gene expression at this stage with in situ hybridization. Similar to alx3 expression at 24 hpf, we observed that alx4a expression is strongest in NCCs residing between the eye and the nasal placode, medial to the nasal placodes, and in cells ventral to the nasal placodes (Fig. S5). Fluorescent in situ hybridization demonstrates that alx4a expression is minimal in the anterior and posterior pharyngeal arches, as predicted from our sequencing analyses. These results confirm our interpretation of the scRNA-seq analyses, that the cells from the cluster enriched for Alx genes reside in the frontonasal domain in intact embryos.
We next examined the temporal expression dynamics of alx3 by performing in situ hybridization with transgenes marking either early migrating NCCs (sox10:EGFP) or post-migratory NCCs (fli1a:EGFP). Early in craniofacial development, alx3 is expressed in an anterior stream of migrating NCCs (Fig. 3A). During morphogenic movement, alx3 is expressed in a subset of anterior NCCs (Fig. 3B). Finally, when NCCs are differentiating into chondrocytes alx3 is expressed in cells near the midline of the roof of the mouth where the ethmoid plate will form (Fig. 3C). These results indicate that alx3 expression is highly enriched in a subset of frontonasal NCCs and continues to be expressed in the anterior neurocranium skeletal cells that they become.
alx3 function is specifically required for anterior neurocranium skeletal development, including establishing distinct edge versus medial ethmoid plate cell morphology
To understand Alx gene function during craniofacial development, we mutagenized the entire Alx gene family using CRISPR/Cas9 (Table S2). The alx1, alx4a and alx4b mutant alleles we generated do not show overt phenotypes as homozygous mutants. Future work outside the scope of this report will use full deletion, or ‘transcriptionless’ alleles (El-Brolosy et al., 2019) to determine definitively any requirement for these genes during zebrafish craniofacial development. Here, we focus on alx3 mutant phenotypes. Two independent predicted loss-of-function alx3 alleles produced robust phenotypes in craniofacial structures derived from frontonasal NCCs (Fig. 4A-C″). Both alleles encode premature termination codons upstream of the homeobox and are predicted to encode truncated proteins lacking the DNA-binding homeodomain (Fig. 4D). Human ALX3 variants associated with FND1 are recessive and are also predicted to disrupt the homeodomain. Larvae from intercrosses between adults heterozygous for alx3 were stained for cartilage and bone at 6 days post-fertilization (dpf) with Alcian Blue and Alizarin Red. For both alleles, loss of one copy of alx3 resulted in skeletal phenotypes at both the anterior and posterior ends of the ethmoid plate. Furthermore, loss of both copies considerably increased penetrance of the defects (Fig. 4E, Fig. S6A).
At the anterior end of the ethmoid plate, chondrocyte cellular morphology was disrupted. In wild types, both edge and medial chondrocyte populations were observed (Fig. 4A-A″), as described above. However, in alx3 mutants all the ethmoid plate cells appeared to have uniform morphology and stained strongly for Alcian Blue (Fig. 4B-C″). To quantify these differences in cell morphology, we measured the length and width of individual chondrocytes positioned at the edges as well as the medial ethmoid plate. We found that whereas wild-type edge and medial cells have significantly different morphologies, edge and medial chondrocyte morphology was indistinguishable in mutants (Fig. S7).
We discovered another neurocranium phenotype near the posterior ethmoid plate. Here, we observed ectopic cartilage and an anteriorly truncated parasphenoid bone (Fig. 4C-C″). To quantify both the anterior and posterior ethmoid plate alx3 mutant phenotypes, we measured the parasphenoid bone area, and the ethmoid plate width and length (Fig. 4F-H). We found that both the parasphenoid bone area and ethmoid plate width were significantly reduced in alx3 mutants compared with wild types. By contrast, ethmoid plate length was significantly increased in alx3 mutants compared with wild-type skeletons. Unlike other zebrafish craniofacial mutants (Kimmel et al., 2020; DeLaurier et al., 2014; Nichols et al., 2016), our anterior neurocranium measurements do not indicate significantly increased phenotypic variation among alx3 mutants compared with wild-type siblings (Table S2). We were able to first observe these cartilage and bone skeletal phenotypes at 3 dpf (Fig. S6B). Both alx3 mutant alleles showed indistinguishable phenotypes with similar penetrance (Fig. 4E), and failure to complement, therefore we focused only on alx3dco3003 for further study.
The anterior neurocranium skeletal phenotypes we discovered motivated us to examine other regions of the craniofacial skeleton for defects associated with the alx3 mutation. We found that skeletal structures of the viscerocranium were indistinguishable between wild types and alx3 mutants. Careful quantification confirmed that the alx3 mutant craniofacial phenotypes are highly specific to skeletal structures of the anterior neurocranium (Fig. S6C). The viscerocranium appeared grossly unaffected by the alx3 mutation, demonstrating that the ethmoid plate and parasphenoid phenotypes are not due to general developmental delay, which would affect the entire skeleton. We also determined that some alx3 homozygous mutants survive to adulthood allowing us to examine maternal-zygotic mutants. We found that maternally deposited alx3 is not a major contributor to craniofacial development (Fig. S6D).
NCC migration and morphogenetic cell movements do not require alx3 function
Morpholino-based studies injecting reagents directed toward the alx3 paralog alx1 showed global defects, including heart edema. The craniofacial phenotypes in these published morpholino experiments include absent or hypoplastic skeletal elements from both the viscerocranium and the neurocranium owing to defective NCC migration (Dee et al., 2013). To test whether early NCC migration or later tissue morphogenesis is defective in the alx3 germline mutants, we analyzed time-lapse recordings of transgenically labeled NCCs using sox10:RFP and fli1a:EGFP. Post-NCC migration, at 24 hpf, the alx3-expressing frontonasal NCCs are medial, lateral and ventral to the nasal placode (Fig. 2C-E). In our time-lapse recordings, we observed similar numbers of NCCs medial, lateral and ventral to the nasal placode between wild-type and alx3 mutant embryos (Fig. 5A). Automated counting of six wild-type and six mutant 24 hpf animals demonstrated no significant difference in the number of frontonasal NCCs by t-test (Table S2), indicating that frontonasal NCC migration does not require alx3 function. We next determined whether post-migratory morphogenic cell movement is affected in alx3 mutants, as has been observed in other zebrafish craniofacial mutants (Sasaki et al., 2013). To this end, we quantified multiple cell movement parameters with cell-tracking software. Every aspect of cell movement we analyzed was unaffected in alx3 homozygous mutants compared with wild types (Fig. 5A,B). We next examined whether NCC proliferation or survival require alx3 function. Cell proliferation and apoptosis assays revealed that the number of proliferating cells and cells undergoing apoptosis were indistinguishable between wild types and mutants (Figs S8 and S9). Thus, we conclude that alx3 function is not required for NCC movement, proliferation or survival.
Chondrocyte differentiation timing requires alx3 function
The ethmoid plate consists of two distinct chondrocyte populations with distinct differentiation timing. Specifically, edge cells differentiate earlier than the medial cells in the developing ethmoid plate (Dougherty et al., 2013). In alx3 mutants, we observed an excess of cells with edge-like strong Alcian Blue staining in the ethmoid plate (Fig. 4B-C″). Therefore, we predicted an excess of cells with edge-like differentiation timing in alx3 mutants. To test this, we analyzed time-lapse recordings of progenitor cells differentiating into chondrocytes during ethmoid plate development. The sox10:mRFP transgene marks the NCCs (Figs 1, 3 and 5), whereas the sox9:EGFP transgene marks cells with chondrocyte identity. We can use these markers to quantify the number of progenitors versus chondrocytes in the ethmoid plate. Time-lapse recordings with these transgenes indicated that there are significantly fewer progenitors, and more differentiated chondrocytes in alx3 mutants compared with wild types early in development (Fig. 6A,B). However, at later time points these differences were no longer detected, suggesting that progenitor cells may be differentiating prematurely in alx3 mutants. We next quantified the ratio of differentiated chondrocytes versus progenitor cells in these live-animal recordings. We found that the ratio of differentiated cartilage cells to progenitor cells was significantly increased in alx3 mutants early in development, confirming that there is altered chondrification timing in mutants. However, the ratios were not significantly different at later stages (Fig. 6C). Thus, early in NCC differentiation, there are alx3 mutant cells that precociously differentiate into chondrocytes of the ethmoid plate. These differences are not observable later when wild types have ‘caught up’ to the early differentiating cells in mutants. At these early stages, we observed accelerated chondrocyte differentiation in both edge and medial regions, suggesting that both the trabecular and ethmoid plate chondrocytes appear to differentiate precociously.
We next examined later stages for regional differences in chondrogenesis timing. To this end, we quantified sox9:EGFP fluorescence along a vector spanning three regions of the mediolateral axis of the ethmoid plate: left edge, medial domain and right edge at 3 dpf (Fig. 6D,E). These experiments revealed that in wild-type larvae, the amount of sox9:EGFP signal is significantly higher in the edge regions compared with the medial region. Thus, in wild types sox9-expressing chondrocytes are primarily detected in the edge regions and largely absent from the medial region at this stage. In contrast, differentiating chondrocytes in alx3 mutants are similarly detected in all three regions; there is no significant difference in maximum fluorescence between the edges and the medial zones (Fig. 6D,E). Moreover, there are significantly more chondrocytes in the medial region of alx3 mutants compared with wild types. These data suggest that in wild types early chondrocyte differentiation is largely restricted to the edges, whereas in alx3 mutants, the chondrocytes differentiate as a uniform front at both the edges and at the midline.
alx3 may function as part of a mediolateral patterning code for the anterior neurocranium
The expression of alx3 is highly specific to the frontonasal population of NCCs (Fig. 7A). We find that alx3 regulates the regional timing of frontonasal NCC differentiation into chondrocytes. Furthermore, the morphology of the differentiated chondrocytes in alx3 mutants is more uniform across the mediolateral axis compared with wild types. These discoveries motivate the hypothesis that alx3 function is required in NCCs for specifying different cell population identities within the anterior neurocranium (Fig. 7B). In a future study, molecular markers of cell identity would formally test this model. Profiling wild-type and alx3 mutant cells with scRNA-seq might reveal transcriptomic signatures of identity transformation. The wild-type scRNA-seq dataset we present here may also support this model. We employed pseudotemporal and trajectory analyses to generate a predicted lineage hierarchy for the frontonasal population (Fig. 7C). These analyses predict which frontonasal cells are more progenitor-like and those that are more mature or differentiated. The lineage tree for this population includes a single ‘branch node’ inferred at the most immature cell state, representing a cell fate decision. Strikingly, one lineage radiating from the immature cell state is enriched for alx3 expression. These data motivate the hypothesis that alx3 functions to promote differentiation down one trajectory. Future studies will directly test a model in which alx3 functions in the frontonasal NCC lineage decision between medial, late-differentiating, cuboidal chondrocytes versus edge, earlier-differentiating, columnar chondrocytes (Fig. 7D).
Within this context, it is tempting to speculate that the entire Alx gene family functions as a mediolateral identity code in frontonasal NCCs. This hypothetical code is reminiscent of the Dlx gene dorsoventral code or the Hox gene anteroposterior code in the pharyngeal arches. Although we find that disrupting some individual Alx members, such as alx1, alx4a and alx4b, do not overtly affect the craniofacial skeleton, future analysis of animals in which multiple Alx genes are disabled in conjunction might reveal interesting phenotypes and further examine the possibility of an Alx code. Nevertheless, our data demonstrate that alx3 function is required during craniofacial development for discrete developmental timing and morphology among different populations of anterior neurocranium skeletal cells.
Anterior neurocranium development is sensitive to alx3 gene dosage
To the best of our knowledge, all disease-associated human ALX3 alleles are recessive. Therefore, it was surprising that both of the alx3 mutant alleles we recovered are partially dominant. Although we have not formally ruled out antimorphic or neomorphic activity of these alleles, we interpret that alx3 is haploinsufficient in zebrafish, or that loss of a single copy affects skeletal development. We found that alx3 is the only member of this family with a single-mutant phenotype in zebrafish, further emphasizing the extreme sensitivity of craniofacial development to alx3 dosage. Given this sensitivity, it is surprising that this gene has been repeatedly lost in evolution. At least three vertebrate genomes, including the amphibian Xenopus tropicalis, the squamate reptile Anolis carolinensis and the bird Gallus, independently lost alx3 (McGonnell et al., 2011). Indeed, alx3 is not essential for zebrafish survival in our laboratory strain. In light of these results, we propose that alx3 gene loss might serve as fuel for craniofacial shape evolution. Loss of alx3 function produces changes in the upper face for natural selection to act upon without producing organismal lethality. Support for a model in which Alx genes function in the evolution of craniofacial shape comes from work demonstrating that the ALX1 locus is associated with divergent beak shape among Darwin's finches (Lamichhaney et al., 2015).
Alx gene functional conservation
Here, we capitalize on the strengths of the zebrafish system to demonstrate that frontonasal cells differentiate more uniformly across the mediolateral axis of the ethmoid plate in alx3 mutants. These observations beg the question: do NCCs in mammalian Alx3 mutants also undergo aberrant cell differentiation? Thus far, this mechanism has not been proposed. However, the phenotypes in both human and mouse largely reflect changes in the skeleton about the mediolateral axis in the frontonasal region, which might result from altered differentiation timing. Future studies are required to determine whether the function of zebrafish alx3 we discovered here is conserved across vertebrates.
MATERIALS AND METHODS
Zebrafish strains and husbandry
All zebrafish experiments used the AB strain. Animals were maintained and staged according to established protocols (Westerfield, 1993; Kimmel et al., 1995). All our work with zebrafish has been approved by the University of Colorado Institutional Animal Care and Use Committee (IACUC). Protocol number for animal research: 00188. The following transgenic lines have been previously reported: Tg(fli1a:EGFP)y1 (Lawson and Weinstein, 2002), Tg(sox10:mRFP)vu234 (Kirby et al., 2006), sox9a:EGFP (sox9azc81Tg) (Eames et al., 2013).
Thirty-two double transgenic flia:EGFP;sox10:mRFP embryos were dissociated into single-cell suspension using cold-active protease from Bacillus licheniformis, DNase, EDTA and trituration (Adam et al., 2017). Approximately 19,000 live, double-positive cells were sorted into tubes pre-coated with fetal bovine serum with a MoFlo XDP100 and subjected to scRNA-seq.
Sorted cells (at 2.94×105 cells/ml with 85% viability) were loaded into the 10x Chromium Controller aiming to capture 4000 cells. Using the NovaSEQ6000, 2×150 paired-end sequencing Chromium 10x Genomics libraries were sequenced. Recovered sequencing reads were mapped to a custom Cell Ranger reference based on the GRCz11/danRer11 genome assembly using the Cell Ranger pipeline (10x Genomics). This reference contained the alx3 transcript from Ensembl release 91 (GRCz10) mapped to the GCRz11 genome in place of the release 101 (GRCZ11) version, due to inaccurate truncation of the 3′ UTR in the most recent version (Fig. S3). The Cell Ranger pipeline removed empty droplets from the dataset based on the distribution of cell barcodes. Additional analysis was performed in R using the Seurat package (Butler et al., 2018; Stuart et al., 2019). Initial quality filtering of the data excluded cells with fewer than 100 genes [>0 unique molecular identifiers (UMIs) per gene]. Each gene in the dataset was required to be present in a minimum of five cells. Regarding the capture of dead or dying cells, a minimum presence of 250 genes (at <2.5% of mitochondrial UMIs detected) were required. After quality filtering, 3329 cells with an average of 2775.822 genes detected per cell were used for analyses. Cell-type clusters were projected and dimensionally reduced via principal component analysis (PCA) and uniform manifold approximation and projection (UMAP). Using the ‘FindMarkers’ function, default, non-parametric, Wilcoxon rank sum tests were used to determine a list of genes differentially expressed per cluster compared with the other three. All of the top candidate genes per cluster had P<0.001. The differentially expressed gene lists generated by this method are presented in Fig. 1 and Table S1. In Fig. 1, we generated the sorted list so a minimum of 10% of cells express genes with positive average log fold change, corresponding to upregulated expression levels. In Table S1, a minimum of 25% of cells must be expressing a given gene with only positive average log fold changes. We found few changes to the lists in using either of the cutoff thresholds. These more stringent lists are presented as a series of plots using Seurat's ‘FeaturePlot’ function (Fig. S1). Visual confirmation of cell type clusters and Alx gene localization was confirmed using Seurat's ‘FeaturePlot’ function. Using Seurat's ‘VlnPlot’ function, violin plots for the four Alx genes were generated to indicate the distribution of cells expressing these genes across clusters. We performed clustering with Monocle 3 software using the Leiden algorithm and dimensions were reduced with PCA and UMAP. The raw, feature-barcode matrix can be accessed from the GEO database (accession number GSE163826).
In situ hybridization, immunostaining, and skeletal preparation
We designed the alx3 probe based on our transcriptomic data. scRNA-seq is biased toward the 3′ UTR of transcripts; therefore, we reasoned that the location where most of the reads aligned would be a suitable region of the transcript to design an in situ probe. Moreover, 3′ UTRs provide the highest degree of sequence specificity. We used 5′-GGCTTCACATGTTTAAGATACAGC-3′ and 5′-TAATACGACTCACTATAGGGATGGGGCTAAAACACAGCAT-3′ to amplify the template for RNA probe synthesis. For alx4a we used a previously published probe (Dee et al., 2013). We performed in situ hybridization as described (Talbot et al., 2010) with minor modifications. After in situ hybridization, animals were immunostained to reveal EGFP protein. Briefly, animals were fixed overnight in BT fix [4% paraformaldehyde (PFA), 4% sucrose, 0.12 mM CaCl2 in PBS] at 4°C. Embryos were permeabilized in 150 mM Tris-HCl pH 9.0 at 70°C for 15 min and acetone for 20 min at −20°C, then washed six times in 0.1% Tween 20 in 1× PBS (1× PBST). Blocking solution (1% bovine serum albumin, 2% goat serum, 1% DMSO, 1× PBST) was applied for 2 h at room temperature followed by 1:2000 mouse monoclonal anti-GFP primary antibody (Clontech, JL-8, 632381) overnight in blocking solution. Embryos were rinsed in 1× PBST five times, then stained with 1:2000 Alexa Fluor 488 donkey anti-mouse secondary antibody (Thermo Fisher, A-21202, lot 536050) diluted in blocking solution for 1 h at room temperature. Embryos were washed five times in 1× PBST before imaging. Skeletal preparations were performed as previously described (Walker and Kimmel, 2007; Brooks and Nichols, 2017) and larval skeletons were scored for phenotype penetrance before genotyping.
The Alx gene family was mutagenized as previously described (Jao et al., 2013; Barske et al., 2016). For alx1, the guide RNA 5′-GGAGAGCAGCCTGCACGCGA-3′ produced a 2 bp deletion in exon 2; the alx1co3002 allele is predicted to encode an early stop codon at amino acid 172. For alx4a, the guide RNA 5′-ATCTCAGCAACCCTTACCC-3′ produced an 8 bp deletion in exon 2; the alx4aco3000 allele is predicted to encode an early stop at amino acid 197. For alx4b, the guide RNA 5′-GGACAGACCTCACTGAAGCA-3′ produced a 7 bp deletion in exon 2; the alx4bco3001 allele is predicted to encode an early stop codon at amino acid 262. All mutant sequences are presented in Table S2. For alx1, alx4a and alx4b, we analyzed skeletal preparations from crosses between germline-stable heterozygotes at 6 dpf. Skeletal phenotypes appeared wild type for all animals in these crosses and genotyping recovered homozygous mutants for alx1, alx4a and alx4b, which were indistinguishable from wild types, so we did not further consider these mutants here. Future work outside the scope of this report will characterize any phenotypes that might arise when these mutant alleles are combined to generate double and triple mutants. For the alx3 mutant alleles we used in this study, exon 2 of alx3 was targeted using the following RNAs: 5′-GGCGTCCAACGGCTTGACTG-3′ and 5′-GATAGGTCAATAGAGTCCGC-3′. Both guide RNAs were injected together along with Cas9 mRNA. G0 injected fish and their F1 offspring were screened for mutagenesis with PCR amplification of a portion of exon 2 with forward primer 5′-CTATCCCGCTCTGGACTCAG-3′ and reverse primer 5′-CAGTGCCAGCTGTTCTCTTG-3′, which flank the targeted sites. Sanger sequencing revealed the sequence changes in the two germline stable alx3 alleles. The first, alx3dco3003, is a 37 bp deletion and a 2 bp insertion resulting in a premature stop codon and a predicted 136 amino acid protein compared with 363 amino acids in the wild type. The second, alx3dco3004, consists of two separate 4 bp deletions resulting in a missense sequence followed by a premature stop codon and a predicted 191 amino acid protein. These predicted truncated proteins lack the DNA-binding and C-terminal domains. For genotyping, heterozygotes and homozygous mutants were discriminated from wild types based on PCR amplification using these same PCR primers. The alx3dco3003 allele produces phenotypes that correlate with genotype in the F3 generation.
Nomarski imaging, fluorescent microscopy, time-lapse recordings and cell tracking
Alcian Blue/Alizarin Red-stained skeletons were dissected and flat mounted for imaging on a Leica DMi8 inverted microscope equipped with a Leica DMC2900 camera. Fluorescent images and time-lapse recordings were captured using a Leica DMi8 microscope equipped with an Andor Dragonfly 301 spinning disk confocal system. Images were analyzed using Imaris and FIJI image analysis tools. Acquisition parameters and fluorescence adjustments were applied linearly and equally to all samples.
Cell proliferation and survival analyses
For the proliferation assay, dechorionated 24 and 48 hpf Tg(sox10:mRFP)vu234 animals were incubated with 400 µM 5-thynyl-2′-deoxyuridine (EdU) in DMSO or DMSO vehicle control in embryo media at 28.5°C. After a 24 h incubation for both time points, larvae were rinsed, and placed back in embryo media for 30 min at 28.5°C. Larvae were anesthetized with tricaine and fixed in 4% PFA for 1 h at room temperature, then rinsed with 1× PBS/1% Triton X-100 and processed according to the manufacturer's protocol for the CLICK-iT-EdU-Alexa 488 kit (Invitrogen). Immunostaining was performed as described above with 1:1000 rat anti-RFP (Chromotek, 5F8, Lot 120110) and 1:2000 goat anti-rat Alexa 568 (Thermo Fisher, A11077, Lot 1853640). For the cell survival assay, dechorionated 24 and 48 hpf Tg(fli1a:EGFP)y1 animals were fixed with 4% PFA for 1 h at room temperature and stained with rabbit anti-cleaved caspase-3 (Asp175) antibody (Cell Signaling Technology, 9661L, Lot 26) and goat anti-rabbit Alexa 568 (Thermo Fisher, A11011, Lot 5436A) as above. For both assays, larvae were imaged as described above.
For penetrance scores, we performed the Fisher's exact test using GraphPad. For morphometric measurements, NCC migration parameters and fluorescence quantifications, the means of all wild types were compared with mutants by unpaired t-test using GraphPad. For vector fluorescence intensity analyses, maximum fluorescence data points were averaged within each genotype and processed with a Kruskal–Wallis test reporting Dunn P-values with and without Benjamini–Hochberg false discovery rate adjustments. All sample sizes and exact P-values are included in Table S2.
We are grateful to Michi Kanai, and Raisa Bailon for careful reading of this manuscript. The rest of the Department of Craniofacial Biology, especially the Artinger laboratory, participated in helpful discussions.
Conceptualization: J.M.M., J.T.N.; Methodology: J.M.M., J.T.N.; Validation: J.M.M., J.T.N.; Formal analysis: J.M.M., J.S., A.T.P., A.E.G., J.T.N.; Investigation: J.M.M., J.T.N.; Resources: E.P.B., J.T.N.; Data curation: J.S., A.T.P., A.E.G.; Writing - original draft: J.M.M., J.T.N.; Writing - review & editing: J.M.M., J.S., A.T.P., E.P.B., J.T.N.; Visualization: J.M.M., J.T.N.; Supervision: J.T.N.; Project administration: J.T.N.; Funding acquisition: J.M.M., J.T.N.
This work was supported by the National Institutes of Health [R56 DE027955, R00 DE024190, R00 DE024190 S1 to J.T.N. and F32 DE029995 to J.M.M.]; and the single-cell RNAseq grant program from the RNA bioscience initiative at Anschutz Medical Campus, University of Colorado to J.T.N. Deposited in PMC for release after 12 months.
All sequencing data are available through Gene Expression Omnibus accession number GSE163826.
Peer review history
The peer review history is available online at https://journals.biologists.com/dev/article-lookup/148/7/dev197483/
The authors declare no competing or financial interests.