Deletion of sf3b4 causes splicing defects and gene dysregulation that disrupt craniofacial development and survival

Nager syndrome (OMIM #154400) is a form of acrofacial dysostosis (AFD), a rare type of disorder characterized by malformations of the craniofacial skeleton and limbs (Halal et al., 1983). Individuals with Nager syndrome specifically present with downslanting palpebral fissures, midface retrusion, micrognathia, defective middle ear ossicles, and hypoplastic or absent thumbs (Trainor and Andrews, 2013). Nager syndrome is a rare disorder, with approximately 100 reported cases ever, worldwide. The craniofacial skeletal structures affected in Nager syndrome are neural crest (NC)-derived.

The NC is an embryonic cell type unique to vertebrates. These cells arise from the neural plate border region of developing embryos, undergo an epithelial-to-mesenchymal transition, migrate throughout the embryo, and eventually give rise to numerous cell types and structures, including melanocytes, the peripheral nervous system in the trunk and most of the craniofacial skeleton in the head. A subdomain of the NC, known as the cranial NC, originating from the mesencephalon and rhombencephalon, migrates as streams to populate the branchial arches. The branchial arches are composed of NC cells surrounding a core of mesoderm and give rise to the majority of the differentiated tissues of the head and neck. The mesoderm-derived cells form the musculature of the head and the endothelial cells of the arch arteries, whereas NC cells develop into craniofacial skeletal elements specific for each branchial arch (Trainor and Krumlauf, 2000). The craniofacial structures affected in Nager syndrome are primarily derived from the first and second branchial arches (Passos-Bueno et al., 2009).

The major cause of Nager syndrome, in approximately 60% of cases, is haploinsufficiency of the splicing factor 3b, subunit 4 (SF3B4) gene (Bernier et al., 2012; Czeschik et al., 2013; Petit et al., 2014). The remaining 40% of cases are of unknown cause. SF3B4 pathogenic variants are predicted to encode truncated proteins with or without altered amino acid sequences at the C-terminus end (Bernier et al., 2012), suggesting that Nager syndrome results from haploinsufficiency. Most cases are sporadic, and both autosomal-dominant and autosomal-recessive inheritances have been reported (Chemke et al., 1988; Hall, 1989). Rodriguez syndrome (OMIM #201170) is another condition due to mutations in SF3B4 (Drivas et al., 2019). Although this syndrome is characterized by similar craniofacial abnormalities as Nager syndrome, the defects are usually more severe and involve lower limb anomalies, cardiac defects and arhinencephaly (Rodríguez et al., 1990). SF3B4 encodes SAP49, a component of the U2 subunit of the major spliceosome (Will and Luhrmann, 2011). The spliceosome is a complex made up of RNAs and proteins that functions to identify non-coding introns in precursor messenger-RNA (pre-mRNA) and promote accurate splicing at the surrounding splice sites. Recognition of the 5′ and 3′ splice sites, as well as proper binding to the pre-mRNA and other parts of the spliceosome complex are required for the splicing process to occur. At each step of the splicing process, different sets of small nuclear ribonucleoproteins (snRNPs) are recruited, and the combination of RNA, small nuclear RNA (snRNA), snRNP and non-snRNP protein interactions allow for splicing to occur at the proper location (Will and Luhrmann, 2011). SF3B4 contains two RNA recognition motifs at its N-terminal end through which it binds upstream to the branch site in the intronic region of the pre-mRNA, where it helps tether the U2 complex to the branch site.

Previous work in mouse models has shown that a loss of SF3B4 affects the axial skeleton and the forebrain, with no apparent craniofacial defects (Yamada et al., 2020; Kumar et al., 2023). This phenotype was associated with a disruption in Hox gene expression and aberrant splicing of several chromatin remodelers (Kumar et al., 2023). In Xenopus laevis, morpholino-mediated knockdown of sf3b4 results in a loss of NC gene expression and causes a reduction or loss of craniofacial cartilages at tadpole stages through a mechanism that involves increased apoptosis (Devotta et al., 2016). A mutation in sf3b4 in zebrafish has been reported to result in a phenotype reminiscent of retinitis pigmentosa, a spliceosomopathy affecting the photoreceptors in the retina (Ulhaq et al., 2023). Despite these efforts, for the most part the precise mechanisms underlying the pathogenesis of Nager syndrome remain largely unknown. It is especially puzzling that variants of a core component of this largely ubiquitous cellular machinery result in such an exquisitely cell type- and lineage-specific defect.

Here, we report the characterization of a Xenopus tropicalis sf3b4 mutant line generated using CRISPR/Cas9 gene-editing technology. We show that in sf3b4 mutants NC induction is unaffected. However, at the tailbud stage these animals show NC migration defects and increased apoptosis in the head region. Although homozygous animals fail to survive beyond the tadpole stage, heterozygous animals are phenotypically largely indistinguishable from wild-type tadpoles. RNA-sequencing (RNA-seq) and Gene Ontology (GO) analyses of embryos at different developmental stages reveal an increase in aberrant splicing events, followed by a massive dysregulation of gene expression in homozygous animals, with an enrichment for terms consistent with defects in NC cell migration and survival as the possible underlying causes of Nager syndrome.

We first analyzed the developmental expression of sf3b4 in Xenopus tropicalis embryos using whole-mount in situ hybridization (WMISH) with a short chromogenic reaction to highlight areas of sf3b4 enrichment. sf3b4 transcripts were first detected at the end of gastrulation [Nieuwkoop and Faber (NF) stage 12.5] in the dorsal anterior ectoderm (Fig. 1A,B). At the early neurula stage (NF stage 14), the sf3b4 expression domain encompassed the neural plate and the neural plate border (Fig. 1C). As neurulation proceeded (NF stage 17), sf3b4 transcripts remained broadly enriched dorsally, including the prospective brain and spinal cord (Fig. 1E), and in a domain that overlapped with the sox10 expression domain in the NC territory (Fig. 1F). As development continued, sf3b4 was detected in the head region, including the developing brain, spinal cord, eyes and migrating NC cells (Fig. 1G-I). At tailbud stages (NF stages 25-31), sf3b4 expression was largely confined to the brain, eyes, otic vesicles, branchial arches and the tailbud (Fig. 1K-M). Control sense probe at stage 14 (Fig. 1D) and stage 25 (Fig. 1J) did not show any signal. Two-color WMISH for sox10 and sf3b4, revealed substantial overlap between these two genes in the migrating NC cells populating the posterior branchial arches (Fig. 1N-R). The expression pattern of X. tropicalis sf3b4 was similar to that of X. laevis sf3b4, with an enrichment dorsally and anteriorly at neurula stage, and in the branchial arches at the tailbud stage (Devotta et al., 2016). Although the expression pattern of SF3B4 in human embryos remains uncharacterized, in the mouse Sf3b4 expression is enriched in the cranial region, including the brain and the branchial arches, as well as the limb bud mesenchyme and the heart at embryonic day 10.5 (Yamada et al., 2020). A longer chromogenic reaction revealed expression of sf3b4 throughout the embryo, including those regions reported in the mouse (C.G., unpublished observations). The sf3b4 expression domain over time encompasses regions of the frog embryo that are broader than the NC territory, therefore this expression pattern cannot account for the cell type-specific effect of sf3b4 mutation in cranial NC and its derivatives.

To generate this custom mutant line, five sgRNAs were synthesized targeting the first three exons of X. tropicalis sf3b4 (Fig. 2A). A mutation in this region has been described in individuals with Nager syndrome from three unrelated families; it is predicted to abolish the methionine initiation codon resulting in a complete lack of protein product (Bernier et al., 2012). Of the 50 embryos injected at the one-cell stage, 20 reached adulthood. After reaching sexual maturity, one female was outcrossed to a wild-type (WT) male and tested for germline transmission. Genotyping of the corresponding F1 embryos showed both −5 bp and −31 bp mutations. The remaining F1 embryos were reared and genotyped as adults by hindlimb web punch sampling. Out of 50 F1 adults, 15 had a confirmed −31 bp mutation. The −31 bp heterozygous F1 sf3b4 adults were intercrossed to generate F2 mutants used for this study. The −31 bp mutation disrupts the start codon of Sf3b4, resulting in a frameshift at the first amino acid, and introducing an early stop codon 29 amino acids downstream (Fig. 2B).

Western blot analysis of stage 25 WT, heterozygous (Het), and homozygous (Null) mutant embryos showed that Sf3b4 protein expression levels are consistent with each genotype and the number of intact copies of the gene (Fig. 2C,D; Fig. S1). By gross morphology, mutant embryos at the tailbud stage (NF stage 25) displayed reduced head length and depth, a phenotype that was more pronounced in Null mutants (Fig. 2E,F).

Because NC-derived craniofacial skeletal elements are primarily affected in Nager syndrome, we analyzed by WMISH the expression of genes to visualize the NC at pre-migratory (snai2, sox10, tfap2e) and migratory (sox9, sox10) stages. Across all three genotypes, we found no difference in the expression of snai2, sox10 and tfap2e at early neurula stage (NF stage 14/15; Fig. 3A,B; Table S1). Because sf3b4 was also expressed in the neural plate (Fig. 1), we examined the expression of sox2 at early neurula stage. This gene was unaffected in mutant (Het and Null) embryos compared to WT sibling controls (Fig. 3A; Table S1). Upon closure of the neural plate (NF stage 20), sox10 expression pattern in the NC territory remained unchanged across genotypes (Fig. 3C; Table S1). At the migratory stage (NF stage 25) there was a notable phenotype in Het and Null mutant embryos characterized by a decrease in the length of a subset of cranial NC streams, as evidenced by sox9 and sox10 expression (Fig. 4A,B; Table S1). Quantification of the phenotype showed a statistically significant reduction of the NC stream length (streams 2 and 3 for sox10, and streams 3 and 4 for sox9) in Het and Null embryos compared to WT siblings (Fig. 4C). The differential effect on sox9 and sox10 may reflect gene-specific sensitivity to the loss of sf3b4. To understand the mechanism driving this phenotype, we performed terminal deoxynucleotidyl transferase dUTP nick end labeling (TUNEL) staining at different developmental stages. We found no change in apoptosis in sf3b4 mutants at NF stage 20 (Fig. S2), which correlates with robust expression of sox10 in the pre-migratory NC territory at this stage (Fig. 3C). However, at the migratory stage (NF stage 25), there was a significant increase in apoptosis in the head region of Het and Null embryos, with a greater number of TUNEL-positive cells in the Null compared to the Het animals (Fig. 4D,E). These results suggest that although NC cell formation is not affected at the neural plate border in sf3b4 mutants, their migration in the branchial arches is impaired.

We next analyzed mutant tadpoles at post-migratory stage (NF stage 40), when NC cells coalesce in the branchial arches into foci of cartilage precursors expressing sox9 and runx2. We noticed that the number of sf3b4 Null animals recovered at this stage was lower than the expected Mendelian ratio (Fig. S3), suggesting that the phenotype described at the migratory stage may not be fully compatible with survival. Our analysis of the sf3b4 mutant tadpoles indicates that sox9 expression is severely decreased in Null compared to Het and WT embryos (Fig. 5A; Table S1). sox9 expression was, however, not completely lost in the mutants, suggesting that the subset of NC cells that migrates into the branchial arches can initiate a cartilage differentiation program. We also analyzed runx2 expression at this stage, a marker for cartilage progenitors. We found that runx2 expression was completely lost in the Null compared to Het and WT embryos (Fig. 5A; Table S1), again indicating that the formation of NC-derived cartilages is severely impaired in Null mutants. We also noticed by gross morphology that the few Null tadpoles that survived up to that stage appeared to present heart defects that were not observed in WT and Het animals (Fig. 5A). To investigate further a possible heart phenotype, we analyzed the expression of several genes that are crucial for cardiac development, namely tbx5, tbx20 and nkx2-5, at multiple stages of development (stages 25, 30 and 35). Whereas tbx5 and nkx2-5 expression was largely unperturbed in Het and Null animals at stage 25 and stage 30, respectively, we observed a marked decreased expression of nkx2-5 and tbx20 in the Null embryos at stage 35 compared to Het and WT embryos (Fig. S4A). Furthermore, transcriptomic analysis of Het versus Null at stage 35 (see below) revealed an enrichment for terms such as ‘cardiac septum morphogenesis’ and ‘heart morphogenesis’ among the downregulated genes (Fig. S4B; Table S2), consistent with a role of Sf3b4 in cardiac development.

Finally, we analyzed the long-term consequences of sf3b4 deletion on craniofacial cartilage formation at NF stage 45 by performing Alcian Blue staining (Fig. 5E). Interestingly, we did not recover any Null tadpoles at this stage; all genotyped animals were either WT or Het (Fig. S3), suggesting that the defects observed at earlier developmental stages are too severe for survival. Head cartilage staining of WT and Het tadpoles revealed that these structures are largely unaffected in Het tadpoles, with a slight decrease in the overall size of the head, which was not statistically significant (Fig. 5B,C). Dissected craniofacial cartilages of sf3b4 Het tadpoles were largely indistinguishable from that of their WT siblings (Fig. 5D).

To identify gene networks disrupted in sf3b4 mutants that may underly the Null phenotype observed, we performed bulk RNA-seq on pools of whole embryos from NF stages 15, 25 and 35 comparing WT, Het and Null at each stage. Principal component analysis (PCA) at all three stages indicated that WT and Het animals are more similar to each other than to the Null (Fig. S5). Differential gene expression analyses indicated very minimal changes at NF stage 15 across genotypes (Fig. 6A, top row). Furthermore, WT and Het samples analyses at all three stages showed little to no difference (Fig. 6A, left column; Table S3), suggesting that loss of one copy of sf3b4 is largely inconsequential at these stages. By contrast, pairwise comparison of Null versus WT and Null versus Het at NF stage 25 and stage 35 showed a significant number of differentially expressed genes (Fig. 6A, middle and right columns). We next used Venn diagrams to show the extent to which gene expression changes in Null versus WT and Null versus Het samples overlap at these two stages (Fig. 6B,C). Interestingly, the majority of differentially expressed genes in Null versus Het largely overlapped with that of Null versus WT at stage 25 (Fig. 6B) and stage 35 (Fig. 6C), although the overall number of differentially expressed genes was much greater at stage 35 (2823 genes) than at stage 25 (352 genes) (Fig. 6C). It is important to point out that the vast majority of these genes were altered in the same direction (up or down) in the two comparisons, as illustrated for a subset of genes in Table S4.

We used the list of overlapping genes at stage 25 (Fig. 6B; 93 genes) and stage 35 (Fig. 6C; 1471 genes) to perform GO analyses for Biological Processes. This revealed an enrichment for terms such as ‘Wnt signaling pathway’ and ‘mRNA splicing, via spliceosome’ for downregulated genes, and terms such as ‘cell division’ and ‘DNA repair’ for upregulated genes at NF stage 25 (Fig. 6D; Table S5). Whereas at NF stage 35, ‘neural crest cell migration’, ‘extracellular matrix organization’, ‘negative regulation of extrinsic apoptotic signaling’ and ‘mRNA splicing, via spliceosome’ terms were enriched for downregulated genes, and ‘cell cycle’, ‘exocytosis’ and ‘visual perception’ terms enriched for upregulated genes (Fig. 6E; Table S6).

We next performed qRT-PCR to validate a subset of downregulated genes in these categories using RNA isolated from new batches of Null and WT stage 35 embryos. Our results revealed a significant decrease in expression levels of tbx1, ednrb and sema4c as representative genes for the ‘neural crest cell migration’ GO term (Fig. 6F), adamts1 and mmp7 as representative genes for the ‘extracellular matrix organization’ GO term (Fig. 6G), and phf5a, eftud2 and sf3b3 as representative genes for the ‘mRNA splicing, via spliceosome’ GO term (Fig. 6H) in Null compared to WT embryos. Taken together, these results demonstrate that gene expression changes in WT versus Het are very minimal at all stages examined, suggesting a limited impact associated with the loss of one copy of sf3b4, whereas when both genotypes were compared to their Null counterpart a similar dysregulation of gene expression was observed, which was significantly more pronounced at stage 35.

Because Sf3b4 is an active component of the splicing machinery, we next examined the number and type of splicing events occurring at each stage (NF stages 15, 25 and 35), comparing each genotype, and focusing on four main events: skipped exons, retained introns, and 3′ and 5′ alternative splice sites (Fig. 7A). As observed for differentially expressed genes (Fig. 6A), we found very little difference between WT and Het when compared to Null at all three stages examined (Fig. 7B; Fig. S6). However, in all cases there was a marked increase in the number of genes with abnormal skipped exon events, especially at stage 25 (Fig. 7B). Venn diagram representation highlights the considerable overlap in genes with aberrant skipped exon between WT and Het when compared to Null at stage 25 (Fig. 7C) and stage 35 (Fig. 7D), with a greater overall number of genes affected at stage 25 (970 genes) than at stage 35 (180 genes) or stage 15 (140 genes) (Fig. 7C,D; Fig. S6).

Interestingly, two-thirds of the genes with aberrant skipped exon events at stage 35 also had skipped exon events at stage 25 (106 out of 163 in WT compared to Null and 106 out of 159 in Het compared to Null), indicating a consistent requirement for Sf3b4 splicing activity (Fig. S7). Genes that show aberrant skipped exon events at stage 25 but not stage 35 may be expressed at low level or not at all at later stages; this is, for example, the case for chd2. This gene has an aberrant skipped exon event at stage 25, but not stage 35 when chd2 expression levels are much lower (Fisher et al., 2023).

We used the list of overlapping genes at stage 25 (Fig. 7C; 837 genes) and stage 35 (Fig. 7D; 142 genes) to perform GO analyses for Biological Processes. The analysis at stage 25 revealed an enrichment for terms such as ‘RNA splicing’, ‘DNA repair’, ‘regulation of embryonic development’ and ‘regulation of apoptotic process’, among others (Fig. 7E; Table S7), and at stage 35 these terms included ‘mRNA splicing, via spliceosome’, ‘clathrin-dependent endocytosis’, ‘mRNA export from nucleus’ and ‘ribosome assembly’ (Fig. 7F; Table S8). We next performed qRT-PCR to validate a subset of genes displaying abnormal skipped-exon events from these categories using RNA isolated from new batches of Null and WT stage 25 embryos. Our results confirmed that eftud2 (RNA splicing), ino80 (regulation of embryonic development), brd8 (regulation of apoptotic process) and ccna2 (cell division) show aberrant skipped exons in Null compared to WT (Fig. 7G). These results indicate that the complete loss of Sf3b4 function affects the spliceosome activity by preferentially promoting atypical skipped-exon events more widely at stage 25 than at stage 35, whereas the loss of one copy of sf3b4 had limited impact on spliceosome activity at both stages.

Here, we report the generation and molecular characterization of a Xenopus tropicalis CRISPR/Cas9 sf3b4 mutant line, as a novel tool to investigate the pathogenesis of Nager syndrome, a condition that affects NC-derived craniofacial structures. We show that in the absence of Sf3b4 function NC induction is not affected. However, NC cell migration is disrupted at the tailbud stage, and coupled with an increase in apoptosis in the head region, a phenotype that is much more pronounced in Null than in Het animals. sf3b4 Null animals failed to survive beyond the tadpole stage, but Het tadpoles were indistinguishable from their WT siblings, forming largely normal craniofacial cartilages. Temporal RNA-seq analysis confirmed at the transcriptome level the limited difference between WT and sf3b4 Het embryos, and revealed aberrant pre-mRNA splicing events in stage 25 Null embryos, most notably characterized by an increase in the number of genes with atypical skipped exons, followed by a massive dysregulation of gene expression at stage 35. GO analysis of differentially expressed genes in sf3b4 Null embryos revealed an enrichment in terms pertaining to mRNA splicing, apoptosis, cell cycle and NC cell migration. Altogether, our data indicate that sf3b4 haploinsufficiency is compatible with normal development in X. tropicalis, whereas homozygous deletion of sf3b4 is detrimental to survival and the formation of NC-derived cartilage progenitors, consistent with an autosomal recessive mode of inheritance. We also found that the sf3b4 Null embryos exhibit heart defects as evidenced by the reduced expression of several cardiac genes at stage 35, and an enrichment for GO terms related to cardiac morphogenesis in the transcriptomic analysis of Het versus Null embryos.

Splicing factor variants have been implicated in numerous diseases, referred to as spliceosomopathies. They form a singular group of diseases due to the nature of their phenotype affecting often a single cell or tissue type, despite the predicted ubiquitous activity of the spliceosome in all cell types of the body (Griffin and Saint-Jeannet, 2020; Beauchamp et al., 2020). The tissues primarily affected in these pathologies include the retina (retinitis pigmentosa), spinal cord (amyotrophic lateral sclerosis), bone marrow (myelodysplastic syndromes), limb and the craniofacial skeleton (Griffin and Saint-Jeannet, 2020). Understanding the cell/tissue specificity of these phenotypes has been challenging.

Mouse embryos with heterozygous Sf3b4 mutation show axial skeleton and forebrain defects with no apparent craniofacial abnormalities, whereas Sf3b4 homozygous deletion is embryonic lethal (Yamada et al., 2020; Kumar et al., 2023). We also report that X. tropicalis sf3b4 homozygous deletion is not compatible with survival, although these animals survive early on, presumably due to a maternal supply of Sf3b4 mRNA or protein (Fisher et al., 2023). Similarly to the mouse phenotype, X. tropicalis sf3b4 heterozygous animals do not show any obvious craniofacial anomalies. Whereas in humans SF3B4 haploinsufficiency causes severe craniofacial abnormalities, the lack of an overt craniofacial phenotype in heterozygous mouse and frog embryos could be explained by the activity of other genes that compensate for the loss of sf3b4.

Interestingly, we also detected a heart defect in null embryos at later stages of development, which we posit contributes to the poor survival of the nulls at tadpole stages. The mechanisms underlying this phenotype remain unclear, and will require further investigation. Although cardiac neural crest cells participate minimally to the Xenopus heart, primarily contributing to the aortic sac and arch arteries (Lee and Saint-Jeannet, 2011), we speculate that cardiac NC cell migration and survival is also affected in the null mutants. Individuals with Nager syndrome do not typically present with cardiac defects; however, individuals with Rodríguez syndrome have cardiac abnormalities, including atrial and ventricular septal defects (Rodríguez et al., 1990).

RNA-seq analyses of mouse Sf3b4 heterozygous mutants reveal a disruption in Hox gene expression and aberrant splicing of several chromatin remodelers known to regulate Hox genes (Kumar et al., 2023). We also found downregulation of the Hox genes hoxa1 and hoxb6 in our RNA-seq data, as well as aberrant splicing of several chromatin remodelers, including kmt2e, kmt2a, phf20, phf10, scmh1, chd7, chd6 and chd2. The same group recently generated a conditional deletion of Sf3b4 in the NC lineage using the Wnt1-Cre2 transgenic line (Kumar et al., 2024). These animals exhibit craniofacial and cardiac defects with variable penetrance and expressivity. The authors propose that the levels of Sf3b4 in non-NC derived neighboring tissues directly influence the severity of the phenotype. At later stages, both NC- and mesoderm-derived craniofacial bones and cartilages were hypoplastic in the mutants, together with abnormal cranial ganglia and cardiac outflow tract defects, presumably due to changes in gene expression related to NC development, cell survival and proliferation (Kumar et al., 2024).

In X. laevis, morpholino-mediated knockdown of sf3b4 results in a reduction of NC gene expression at neurula and tailbud stages, causing a reduction of NC-derived craniofacial cartilages at tadpole stages through a mechanism that involves apoptosis (Devotta et al., 2016). We did not observe a similar reduction in NC gene expression in the sf3b4 X. tropicalis mutant embryos. Although the X. laevis morphant phenotype was partially rescued by expression of human SF3B4 (Devotta et al., 2016), we cannot exclude the possibility that the morpholino may still also have off-target effects that could explain the observed reduced expression of premigratory NC genes. There is also the possibility that compensatory mechanisms may be at play in the CRISPR/Cas9-induced sf3b4 mutation, with a subset of genes compensating at early stages for the loss of sf3b4, a phenomenon that is not seen to the same extent in the context of morpholino knockdowns (Rossi et al., 2015). Regardless, in both models, cell death appears to be a common root cause in the presentation of the phenotype. Although X. tropicalis sf3b4 Het embryos exhibit a mild NC phenotype at tailbud stages, they appear to regulate/compensate this defect as they show virtually no craniofacial differences with WT animals at tadpole stage. Our data indicate that loss of Sf3b4 in X. tropicalis manifests differently than in humans. The sf3b4 Het embryos do not show an overall robust disease phenotype, and in many respects are much more similar to WT than Null embryos, suggesting that redundant mechanisms are at play in Xenopus. This is further supported by our RNA-seq data (see below). Although this tool does not fully recapitulate the human disease condition, it is an excellent system to investigate the molecular processes regulated by Sf3b4 during NC and craniofacial development.

Our RNA-seq analyses at three different stages of NC development – pre-migratory (NF stage 15), migratory (NF stage 25) and post-migratory (NF stage 35) – indicate that there is very little change in gene expression between WT and sf3b4 Het animals (Fig. 6A), suggesting that loss of one copy of sf3b4 is largely inconsequential in X. tropicalis, which is consistent with the phenotypes of these embryos (Fig. 5). When both genotypes were individually compared to sf3b4 Null, we observed a similar pattern of gene dysregulation, with a much greater number of genes affected at stage 35 than at stage 25 or stage 15 (Fig. 6B,C; Fig. S6). GO analyses of these gene sets revealed terms such as ‘mRNA splicing, via spliceosome’, ‘apoptosis’ and ‘neural crest cell migration’ enriched in genes downregulated at stage 35, whereas the most statistically significant GO term enriched among upregulated genes was ‘visual perception’. Genes downregulated related to ‘mRNA splicing, via spliceosome’ include several genes that have been linked to other craniofacial spliceosomopathies, such as eftud2 (Lines et al., 2012), txnl4 (Wieczorek et al., 2014) and phf5a (Harms et al., 2023). We also found multiple genes causative of retinitis pigmentosa (a retina-specific spliceosomopathy, characterized by degeneration of photoreceptors), including prpf3, prpf8 and prpf31 (Tanackovic et al., 2011), and snrnp200 (Zhao et al., 2009). This is consistent with a recent study describing a photoreceptor phenotype in sf3b4 mutant zebrafish larvae (Ulhaq et al., 2023). Also downregulated were genes related to ‘apoptosis’, including cdkn1a (cyclin dependent kinase inhibitor 1A), a negative regulator of p53 (El-Deiry et al., 1994) and tp53 (tumor protein 53) itself (Chipuk et al., 2003). This is significant, as previous work has implicated apoptosis in craniofacial malformations underlying loss of Eftud2, Snrpb and Txnl4a in X. laevis (Park et al., 2022), as well as the specific involvement of the p53 pathway in Eftud2 mutant mice (Beauchamp et al., 2021), indicating a common mechanism underlying craniofacial spliceosomopathies. Under the ‘neural crest cell migration’ GO term, we recovered several members of the semaphorin family (sema3b, sema4c, sema6a and sema6c), a semaphorin receptor (nrp1; neuropilin 1), ednrb (endothelin receptor type B) and two transcription factors, hif1a and tbx1. Most of these factors have been linked to NC development in several species, and pathogenetic variants of TBX1 and EDNRB cause two neurocristopathies, DiGeorge syndrome (Merscher et al., 2001) and Hirschsprung disease (Amiel et al., 1996), respectively. Future studies will determine whether the downregulation of these genes in sf3b4 Null mutants underlies aspects of the NC-specific phenotype of this craniofacial spliceosomopathy.

The RNA-seq analyses combined with rMATS, a computational tool to detect splicing events, revealed very limited changes in the number and type of splicing events occurring between WT and sf3b4 Het animals (Fig. 7B), again suggesting that one copy of sf3b4 in frogs is sufficient for spliceosome activity. When both genotypes are individually compared to sf3b4 Null samples, a similar pattern of disrupted splicing events emerges, with a strong bias toward skipped exons, and a much greater number of genes affected at stage 25 than at stage 35 (Fig. 7C,D). This bias toward skipped exons could be due to the role of Sf3b4 in tethering the U2 complex to the branch site during splicing (Champion-Arnaud et al., 1994). In the absence of Sf3b4, the branch site stability may decrease, leading to the spliceosome falling off and exons skipping as alternate branch points are used instead. Interestingly, however, an enrichment in skipped-exon events has been reported for mutation in other components of the spliceosome, including Eftud2 (Beauchamp et al., 2021) and Snrpb (Alam et al., 2022), which carry very different functions in the spliceosome.

GO analyses identified several terms pertaining to ‘mRNA splicing’, ‘regulation of apoptotic process’, ‘DNA repair’, ‘cell division’ and ‘regulation of embryonic development’ were enriched at stage 25, and at stage 35 the most statistically significant GO term enriched was ‘mRNA splicing, via spliceosome’, followed by ‘clathrin dependent endocytosis’. How these factors relate to NC and craniofacial development remain to be investigated. There is very little overlap between the genes that are dysregulated and the genes that are aberrantly spliced in our dataset. This could be explained by the fact that genes that are aberrantly spliced are still expressed at WT levels, but the expression of their downstream effectors are up- or downregulated due to these changes in splicing events. This question remains open as these mechanistic details require further elucidation.

The temporal transcriptomic analysis described here offers a unique perspective on the sequence of events and the mechanisms of Sf3b4 pathogenesis. sf3b4 mutants show a massive disruption in spliceosome activity occurring primarily between NF stage 15 and stage 25, followed by a significant dysregulation of gene expression occurring around NF stage 35. This suggests that changes in gene expression are the likely consequence of spliceosome activity disruption and therefore may not be the primary cause of the phenotype, but rather the contributing factor. Future studies will focus on the plethora of genes misregulated in Null mutants at stage 35 (Table S6) to identify the gene networks and pathways underlying this craniofacial condition.

X. tropicalis were obtained from The National Xenopus Resource (NXR; Woods Hole, MA, USA) housed as described (McNamara et al., 2018; Shaidani et al., 2020). Adult sf3b4 −31/+ mutants were intercrossed by in vitro fertilization to produce embryos for this study. Females were given 20 U of pregnant mare serum gonadotropin (BioVender, RP17827210000) and 200 U of human chorionic gonadotropin (BioVender, RP17825010) to induce egg laying (Wlizla et al., 2018).

Antisense digoxygenin (DIG)-labeled probes (Genius kit, Roche) were synthesized using template cDNA encoding X. tropicalis sf3b4 (Horizon, MXT1765-202789577) and X. laevis sox10 (Aoki et al., 2003), snail2 (Mayor et al., 1995), tfap2e (Hong et al., 2014), twist1 (Hopwood et al., 1989), sox9 (Spokony et al., 2002), runx2 (Kerney et al., 2007), tbx5 (Brown et al., 2005), tbx20 (Brown et al., 2002) and nkx2-5 (Cleaver et al., 1996). WMISH was performed as previously described (Harland, 1991). For double in situ hybridization, sf3b4-DIG- and sox10-FITC-labeled RNA probes were sequentially detected using FITC- and DIG alkaline phosphatase-conjugated antibodies. sox10 was visualized using Magenta Phosphate (Biosynth) and the color reaction to detect sf3b4 was performed using 5-bromo-4-chloro-3-indolyl-phosphate (Roche), as described (Saint-Jeannet, 2017).

Guide RNAs (sgRNAs) were designed utilizing CRISPRScan (https://www.crisprscan.org/) targeting the first three exons of sf3b4 (Moreno-Mateos et al., 2015): T1: GGTGCCACGGTGTATGTCGG, T2: GGAACATGATAAAGCTCTAT, T3: GGGGTCTCTCATGATCTTGG, T4: GGCTTCGGACGCAGCCATTG and T5: GGGGCGAGAAAATGGCGGCT. 5′ dinucleotides were converted to GG for increased mutagenic activity (Gagnon et al., 2014). The SP6 MEGAscript kit (Ambion, AM1330) protocol was followed for synthesis of all sgRNAs. Injections using the X. tropicalis Nigerian line (RRID: NXR_1018) at the one-cell stage consisted of 500 pg of guide RNA and 1000 pg of Cas9 protein into each embryo. One F0 female was outcrossed by in vitro fertilization to a WT male, resulting in F1 embryos with −5 bp or −31 bp mutations. The −31/+ (Xtr.sf3b4^emNXR; RRID: NXR_3056) sf3b4 mutant line is available through the NXR (https://www.mbl.edu/xenopus).

To genotype adult frogs, tissue samples were taken from webbing on hindlimbs using single-use biopsy punches (VWR, 21909-140). Heterozygous mating pairs were intercrossed to produce embryos for in situ hybridization and RNA-seq analyses. Embryos collected for in situ hybridization were fixed in MEMFA (10 ml 10× MEMFA salts, 10 ml 37% formaldehyde, 80 ml nanofiltrated H₂O) overnight at 4°C. Embryos were stored long term at −20°C in 100% ethanol. Ventral tissue (stages 15-20) and tail clips (stages 25-45) for genomic DNA (gDNA) extraction were taken from fixed embryos after rehydration in PBS. For embryos used in RNA-seq analysis, the ventral (stage 15) and posterior (stage 25 and stage 35) parts of the embryos were dissected and collected for genotyping, and the remaining tissues were immediately preserved on dry ice and stored at −80°C.

gDNA extractions were performed using the QIAGEN DNeasy Blood & Tissue Kit (69506). PCR amplification was carried out using the following primers for the targeted region: forward primer 5′-AATGAAACACCCTCTATGCGC-3′ and reverse primer 5′-AGAGATGGAGCCTGCACC-3′. The PCR product was then purified by following the NucleoSpin PCR Clean-up procedure (Macherey-Nagel; 740609.250) and sent to Genewiz, South Plainfield, NJ, USA and The Keck Facility, Woods Hole, MA, USA for Sanger sequencing to confirm genotypes.

Pools of ten embryos were homogenized in lysis buffer and concentrated, and western blot analysis was performed as previously described (Devotta et al., 2016). Primary antibodies were: anti-Sf3b4 polyclonal antibody (Proteintech, 10482-1-AP; 1:2000) and anti-α-tubulin antibody (Sigma-Aldrich, T9026; 1:500). Secondary antibodies were: donkey anti-rabbit (EMD Millipore, MAB201P) and donkey anti-mouse IgG (Abcam, ab6820) coupled to horseradish peroxidase (1:10,000).

Using ImageJ software, a vertical line was drawn behind the eye to define head depth, and a horizontal line extending from the back of the eye to the front of the face from the eye to define head length (Fig. 2E, top embryos).

Total RNA was extracted from embryos using the RNeasy microRNA isolation kit (QIAGEN) and the RNA samples were digested with RNase-free DNase I to eliminate genomic DNA. RT-qPCR analysis was performed using the Power SYBR™ Green RNA to C_T™ 1 step RT-PCR kit (Applied Biosystems, 4389986) on a QuantStudio 3 Real-Time PCR System (Applied Biosystems) using gene-specific primer sets (Table S9).

The in situ hybridization signal was measured using ImageJ software. The image was reverted to an 8-bit gray image and inverted. The area of staining was selected, and the mean gray scale was calculated for each area. The left and right sides of each embryo were calculated separately and then averaged together for the final value. Two-tailed Welch's unpaired t-test was performed using Prism GraphPad (v.10.01) to determine statistical significance.

The length of NC streams was measured using ImageJ software by drawing a vertical line between the dorsal midline and the ventral-most aspect of the sox10 and sox9 expression domains in each stream (streams 2, 3 and 4). The ratio of stream length to total embryo depth (dorsal-ventral axis) was then calculated. Two-tailed Welch's unpaired t-test was performed using Prism GraphPad (v.10.0.1) to determine statistical significance.

TUNEL staining was carried out as previously described (Hensey and Gautier, 1998; Devotta et al., 2016). To quantify cell death, embryos were imaged individually, and a set region of the head encompassing the branchial arches and excluding the eye was chosen in which the number of TUNEL-positive cells were manually counted for each embryo. Two-tailed Welch's unpaired t-test was performed using Prism GraphPad (v.10.0.1) to determine statistical significance.

Tadpoles at NF stage 45 were fixed in MEMFA overnight at room temperature, rinsed in tap water, skinned and eviscerated. Dissected heads were then stained in Alcian Blue solution (0.06% in 30% acetic acid and 70% ethanol) for 24 h, rinsed in 95% ethanol for 2 h, rehydrated and macerated in 2% potassium hydroxide overnight at room temperature. Specimens were then cleared in increasing concentrations of glycerol (20%, 40%, 60% and 80%) in 2% potassium hydroxide and imaged on a Leica M165 stereomicroscope (Leica Microsystems Inc.).

Post-genotyping, pools of five embryos were used for RNA extraction using the QIAGEN RNeasy Micro Kit (74004). RNA was eluted in RNase-free water and sent to the New York University Genome Technology Center for RNA sequencing. RNA QC/QA was performed on a bioanalyzer before performing automated stranded RNA-seq library preparation with polyA selection. Samples were run on the Illumina NovaSeq 6000 system with SP 100 cycle flow cells.

RNA-seq data were analyzed using the sns rna-star pipeline (https:igordot.github.io/sns/routes/rna-star.html). Adapters and low-quality bases were trimmed using Trimmomatic (v.0.36) (Bolger et al., 2014). Sequencing reads were mapped to the reference genome (XENTR_10.0; https://download.xenbase.org/xenbase/Genomics/JGI/Xentr10.0) using the STAR aligner (v.2.7.3) (Dobin et al., 2013). Alignments were guided by a Gene Transfer Format (GTF) file, which was converted from a General Feature Format (GFF3) file. The mean read insert sizes and their standard deviations (+s.d.) were calculated using Picard tools (v.2.18.20) (http://broadinstitute.github.io/picard). The genes-samples counts matrix was generated using featureCounts (v.1.6.3) (Liao et al., 2014), normalized based on their library size factors using DEseq2 (v.1.30.1) (Love et al., 2014), and differential expression analysis was performed. The read per million normalized BigWig files were generated using deepTools (v.3.1.0) (Ramirez et al., 2016). To compare the level of similarity among the samples and their replicates, we used two methods: PCA and Euclidean distance-based sample clustering. All the downstream statistical analyses and generating plots were performed in R environment (v.4.0.3) (https://www.r-project.org/). Differentially expressed genes were determined with a cutoff of less than −0.4 or more than 0.4 log2FC.

rMATS (v.4.0.2) (Shen et al., 2014) was used for detecting differential alternative splicing events from RNA-seq data. Events with a mean of inclusion junction counts or a mean of skipped junction counts <10 were removed. Only events with an inclusion level difference of more than 0.1 or less than −0.1, with P<0.05 were included.

The results of GO analysis were generated by DAVID 2021 (Huang et al., 2007) using gene lists from differential expressed genes and splicing events with conversion to human terms. Non-conserved terms were removed from analysis. Venn diagrams were generated using DeepVenn (Hulsen et al., 2022 preprint).

We thank members of the J.-P.S.-J. laboratory past and present for their support and helpful discussions, and Benjamin McKenzie (Eugene Bell Center for Regenerative Biology and Tissue Engineering, and National Xenopus Resource, Marine Biological Laboratory, Woods Hole, MA) for technical help. We are grateful to the NYU Langone Genome Technology Center for expert library preparation and RNA sequencing. The work benefited from the support of Xenbase (http://www.xenbase.org/; RRID:SCR_003280).