Hirschsprung disease (HSCR) is the most common form of congenital obstruction of the bowel, with an incidence of ~1/5000 live births. The incidence varies significantly between ethnic groups, the highest being in Asia (2.8/10,000 live births) [1]. HSCR results from a failure of the neural crest (NC) cells, which give rise to the enteric nervous system (ENS), to migrate, proliferate, differentiate, or survive in the bowel wall, resulting in aganglionosis of the distal part of the gastrointestinal tract. This results in clinically severe and sometimes life-threatening bowel obstruction. As HSCR is a highly heritable disorder, genetic variation (mutations) in the genomes of these patients must largely explain disease development. So far more than 15 HSCR susceptibility genes, six linkage regions [1], and three associated loci [2, 3] have been found. The genes identified belong to a limited number of pathways relevant to the development of the ENS, among which the RET and the endothelin pathways are the most important. RET (encoding a tyrosine kinase) is the major gene with >80% of all known mutations. These have been mainly identified in ~50% of familial (mostly long-segment HSCR (L-HSCR), total colonic aganglionsis (TCA)) and up to 20% of sporadic (mostly short-segment HSCR (S-HSCR)) cases [4]. However, the identified genes and variants explain no more than 25% of the overall genetic risk of all HSCR cases [2, 3]. These findings indicate that the majority of the disease risk must be due to as yet unidentified rare or common variants in the known HSCR genes or, more likely, variants in yet unknown genes, acting alone or in combination.

The most popular approach to the analysis of whole exome sequencing (WES) data includes selecting genes that can be functionally linked to the pathways already known to be involved in the disease under study. Variants in genes totally unlinked to the known genes or pathways are largely neglected. This study aimed to determine the contribution of rare exonic, non-synonymous DNMs to HSCR without any a priori selection. Therefore, not only did we perform “standard” exome sequencing analyses, followed by burden tests and in silico prediction, but we also carried out an unbiased in vivo analysis of the mutated genes in a zebrafish model.

Over the past years a large number of papers have been published on de novo mutation screening in human diseases. This has resulted in the identification of many new disease-associated genes. Genes are considered as true disease causing when there is a significant excess of de novo mutations among unlinked patients. This works well for diseases that are relatively homogeneous or for which many patients can be investigated. For the more heterogeneous rare diseases for which only small cohorts are available this poses a problem. Often possible disease causing genes are found in a single patient. How to decide whether this finding is of importance? Expression of the gene in the relevant tissues can be considered as additional evidence, as is network analysis. However, making strong statements for private disease genes is, and will be, extremely difficult. It also results in a bias towards genes in the known disease-causing gene networks. Genes not fitting the current knowledge are often discarded as uninteresting. In the current study we wanted to take this all one step further.

Therefore, to assess if the mutated genes play a role in ENS development, we performed two rounds of functional analyses. We opted for an in vivo approach using the zebrafish model system. We knocked down the expression of zebrafish orthologs of 12 of the 13 genes in which loss-of-function or missense DNMs had been identified. Nine genes were successfully knocked down and four of them resulted in loss of neurons in the distal gut. It was noted that in some morphants the proportion of larvae displaying a HSCR-like phenotype was comparatively low. Although in most cases we could improve the efficacy by increasing the amount of morpholino injected, this also led to a higher death rate and more dysmorphic larvae caused by the general toxicity of the morpholino. Therefore, we chose the dose that caused the minimum death rate and dysmorphic rate with slightly lower, but still valid, efficacy. It should also be noted that in some of the morphants the gut development appeared to be affected. Perhaps it was not surprising as from the in situ hybridization results some of these genes seemed to be expressed in the surrounding intestinal tissues as well. Therefore, when hypothesizing the function of these genes in ENS development, in addition to the intrinsic effect on enteric neural crest cells, one should not rule out the possibility of the extrinsic influence via the intestinal tissues. Noteworthy, the SBMOs targeting three of the orthologs (aff3, scube3, and vezf1a) did not knockdown the target transcripts as expected, highlighting the limitation of morpholinos [23]. To strength our morpholino data we opted for CRISPR knockout, which allowed phenotype analysis in F0 larvae. The HSCR-like phenotype was reproduced in the CRISPR knockout of dennd3a and b, ncl1, nup98, and tbata. Although the presence of indel mutation was confirmed by T7E1 assay in all larvae injected with gRNA, the number of larvae displaying the phenotype was relatively low. We suspected this was mainly due to the varying CRISPR efficiency (number of indels present in the target gene over the total number of the target gene copy in a larva) between larvae.

For those four mutations from the newly validated genes (DENND3, NCLN, NUP98, and TBATA), two of them (NCLN:Q166* and DENND3:K640fs) are loss-of-function mutations that disrupt gene translation; the other two (NUP98: N1662S and TBATA:R53C) are predicted as highly deleterious by Polyphen2 and a logistic regression model that combine different predicted scores [24]. This in silico evidence (Table S8) strongly supports the potential impact of the mutations on their genes [15]; hence, we do not report any functional assay to validate them in this study.

The finding of these four genes recapitulating HSCR in zebrafish clearly demonstrates that genes that would have never been followed up based on the usual gene selection criteria should not be ignored. Using bioinformatics prediction and statistics, we would have focused on RET and CKAP2L as they were the only genes significantly enriched for rare variants in patients.

We wondered whether any of these four genes could be linked to the ENS or whether they play relevant roles in neuronal development or NC-derived cell types in general. In fact, by studying these genes in more depth we noticed that all four, despite a lack of obvious connection to the known ENS pathways, are involved in the development of the CNS or the NC, making these not as random as they might first appear.

DENN/MADD domain containing 3 (DENND3) is a guanine nucleotide exchange factor (GEF) that is involved in intracellular trafficking by activation of the small GTPase RAB12 [25]. In zebrafish, Rab12 and other Rab GTPases are highly expressed by pre-migratory NC cells and their expression is dysregulated in Ovo1 morphant zebrafish that display altered migration of NC cells [26]. Independently of RAB12, DENND3 also regulates Akt activity, which is involved in the proliferation and survival of enteric NC cells [25, 27].

Nicalin (NCLN) is a key component of a protein complex that antagonizes Nodal signaling [28], which in vertebrates is involved in induction of the mesoderm and endoderm [29]. In contrast, inhibition of Nodal signaling is required for the specification of human embryonic stem cells into neuroectoderm, including the NC [30, 31]. The antagonizing function of Nicalin on Nodal signaling is therefore consistent with the NC specification required for ENS development.

NUP98 encodes a precursor protein that is autoproteolytically cleaved to produce two proteins: NUP98 from the N-terminus and NUP96 from the C-terminus [32]. A missense DNM was identified in the last exon of the NUP98 gene and therefore affects the NUP96 protein. As in humans, zebrafish Nup96 is produced by cleavage of the Nup98 precursor protein. Since morpholinos act on the mRNA level, both nup98 and nup96 were targeted in our zebrafish experiments. It is therefore unclear whether the observed aganglionosis is caused by loss of Nup98 or Nup96. NUP96 is one of approximately 30 proteins in the nuclear pore complex (NPC) [33] and its expression level regulates the rate of proliferation [34]. Two other members of the NPC (Nup133 and Nup210) are involved in neural differentiation in mice [35, 36]. Moreover, NUP96 interacts with NUP98 and NUP98 is involved in transcriptional regulation of the HSCR genes SEMA3A, DSCAM, NRG1, and the NRG1 receptor ERBB4 in human neural progenitor cells [37]. Therefore, it is likely that loss of both NUP proteins (NUP96 or NUP98) could contribute to HSCR development.

The mouse ortholog of TBATA (Thymus, brain and testes associated) is called Spatial and is highly expressed during differentiation of several tissues [38]. These include the cerebellum, hippocampus, and Purkinje cells in the brain, where TBATA/Spatial is expressed in early differentiating neurons [39]. In mouse hippocampal neurons, TBATA/Spatial is required for neurite outgrowth and dendrite patterning [40].

The four newly identified HSCR candidate genes seem to play a role in neuronal development and could potentially be involved in HSCR (Fig. 3). This also suggests a clear link between CNS and ENS development. This does not come as a total surprise, as several known HSCR genes (e.g., KBP, SOX10, NRG1, IKBKAP, ZEB2, PHOX2B) are involved in both CNS and ENS pathologies [2, 41–43] and the fact that HSCR is strongly associated with Down syndrome.

Besides the fact that several HSCR/neuromuscular genes are known to be associated with CNS defects, the opposite is also described. Many neurological and psychiatric disorders are associated with constipation, and sometimes defects in the ENS are reported [44]. For instance, it has recently been described that mutations in CDH8 result in a specific subtype of autism in combination with gastrointestinal problems. A cdh8 ^−/− zebrafish recapitulates the human phenotype, including increased head size, impairment of gastrointestinal motility, and reduction in post-mitotic enteric neurons [45]. Besides, searching “CNS” and “autism” in Phenolyzer (http://phenolyzer.wglab.org/) returns two ENS genes, APP [46] and MECP2 [47].

Given the fact that HSCR occurs together with neurological disorders more often than expected, it is not surprising that dysfunction of these newly identified neurological-related genes results in dysregulation of the NC-derived cells that form the ENS, and hence in HSCR. These data are further corroborated by the expression patterns we observed for the orthologs of these four genes in zebrafish embryos (Fig. 2), with all four being expressed in both brain and gut.

Using zebrafish as a model we have experimentally shown that some of the genes with de novo mutations appear to be functionally or biologically linked with the ENS. As the effect of de novo variants on the phenotype may depend on the genetic background, it is tempting to speculate that those genes that failed to reproduce the HSCR phenotype in zebrafish could, in fact, contribute to the disease in humans.

While statistical evidence of involvement in HSCR has only been attained for RET, the functional studies presented here support the possible contribution of DENND3, NCLN, NUP98, and TBATA to the disorder. Finding a role for these genes in ENS development will open new research avenues and enhance our knowledge about ENS development and HSCR disease mechanisms. Until now, we believed that the number of cellular processes involved in the development of HSCR was limited. Clearly this idea needs to be revisited as the novel genes we identified are not directly linked to any of the currently known HSCR gene networks.

Our data open new fields of investigation into HSCR pathology and provide insight into the development of both the ENS and CNS. Moreover, as a lesson for all those who work on rare heterogeneous diseases, the study demonstrates that functional analyses of genes carrying DNMs is warranted to delineate the full genetic architecture of rare complex diseases.