(A) For this study we included 326 controls (group 4) and 58 HSCR patients. We determined the RET and / or known disease gene coding mutations of HSCR patients (n = 58). The HSCR patients with coding mutations were included in group 2 (n = 15). We determined the presence of associated anomalies in the rest of the HSCR patients, including them in either the group containing associated anomalies (HSCR-AAM; group 1, n = 23) or in the group without associated anomalies (group 3, n = 30). (B) For all subgroups of HSCR patients and the controls Copy Number profiles were determined. To select candidate genes, we ranked CNVs according to their frequency in the unaffected controls from the Deciphering Developmental Disorders project (for illustrative purposes a screenshot of their UCSC genome browser track; https://genome-euro.ucsc.edu/ is depicted). We determined if genes included in the rare CNVs were “ENS genes”: genes with increased expression in isolated ENS (with and without the addition of GDNF) cells compared to the whole intestine (for illustrative purposes a screenshot of a heatmap of differential gene expression of known disease genes derived from brb-array tools (https://brb.nci.nih.gov/BRB-ArrayTools/) is depicted. Additionally, we determined if a gene was known as a Constrained Coding Region (CCR)(For illustrative purposes a screenshot of gnomAD browser website where such scores can be found; https://gnomad.broadinstitute.org/). Next, we determined whether disruption of the main candidate genes resulted in a reduction of enteric neurons in zebrafish. (C) In parallel, we evaluated the contribution of predisposing haplotypes across groups. Combined, we suggest that the RSnc (predisposing haplotypes), RSrcv (deleterious rare variant burden) and RScnv (deleterious Copy Number Variation) result in the genetic risk for Hirschsprung disease. Abbreviations: CNV = Copy Number Variation, ENS = Enteric Nervous System, RSnc = Risk Score non-coding variants, RSrcv = Risk Score rare coding variants, RScnv = Risk Score CNV, CCR = Constrained Coding Region.

(A) Graph showing the number of rare CNVs found in the four groups. No statistical differences were found. (B) Graph showing the size in base pairs (bp) of the rare CNVs found in the four groups. Rare CNV size in group 1 (HSCR-AAM) was significantly larger than those found in group 3 (HSCR without a causal variant identified) and group 4 (controls). (C) Graph showing the size (bp) of the rare CN losses found in the four groups. Rare CN loss size in group 1 (HSCR-AAM) was significantly larger than those found in group 4 (controls). (D) Graph showing the number of “ENS genes” included in the CNVs found in the four groups. The number of “ENS genes” included in CNVs of group 1 (HSCR-AAM) was significantly higher than those found in group 4 (controls). (E) Graph showing the number of “ENS genes” included in the CN losses found in the four groups. The number of “ENS genes” included in CN losses of group 1 (HSCR-AAM) was significantly higher than those found in group 4 (controls). (F) Graph showing the number of CCR genes included in the CNVs found in the four groups. The number of CCR genes included in CNVs of group 1 (HSCR-AAM) was significantly higher than those found in group 4 (controls). (G) Graph showing the number of CCR genes included in the CN losses found in the four groups. The number of CCR genes included in CN losses of group 1 (HSCR-AAM) was significantly higher than those found in group 4 (controls). (H) Graph showing the number of “ENS genes” that are also a CCR, included in the CN losses found in the four groups. These genes were considered as the HSCR candidate genes. The number of CCR “ENS genes” included in CNVs of group 1 (HSCR-AAM) was significantly higher than those found in group 4 (controls). Each dot represents one patient. Error bars represent standard deviation. Statistical analysis used: one-way ANOVA followed by students t-test. ** p < 0.01 *** p < 0.001; Exact p-values and statistics can be found in S2 Table.

(A) Visual representation of the CRISPR/Cas9 complex injections using the Alt-R CRISPR-Cas9 System from Integrated DNA Technologies (IDT) [34]. (B) Legend showing the various phenotypes observed in zebrafish. Larvae are characterized as category I (e.g. ENS not affected) when there is a full colonization of the intestine by phox2bb-GFP+ cells which are dispersed at a density comparable to that observed in untreated animals. Larvae are characterized as category II when there is an overall reduction observed in the presence of phox2bb-GFP+ cells along the total length of the intestine. Larvae are characterized as category III when an absence of phox2bb-GFP+ cells is observed in the most distal part of the gut, depicted by the white arrowhead and magenta line. Category IV classification is used when the phox2bb-GFP+ have migrated into the intestine but only reach in the mid-gut as in the example depicted by the white arrowhead. The green arrowheads as well as the magenta line indicate the range in which fish are categorized as category IV. The most severe alteration in the ENS is categorized as V, where the mid-gut and distal-gut are deprived of neurons, as depicted by the magenta line and the ENS has migrated no further than the white arrowhead. (C) Disruption of four genes induced ENS phenotypes in the zebrafish larvae, of which two were statistically significant (accumulated data from multiple experiments). (D) In line with the observed ENS phenotypes, quantification of the number of neurons per 100μm showed significantly reduced numbers of neurons upon disruption of three genes (accumulated data from multiple experiments). (E) Schematic representation of sequential injection with CRISPR/Cas9 complexes and ATG-blocking ret morpholino. Graph showing accumulated data of the percentage of fish with ENS phenotypes upon injection of a morpholino targeting Ret translation at a concentration that induced ENS phenotypes in approximately 50% of the fish. Ret morpholino injections in combination with disruption of mapk8a/b and tubb5 shows epistasis. Disruption of gnl1 showed a trend towards higher penetrance of ENS phenotypes. Statistical tests used: students t-test (D), one-way ANOVA followed by Dunnett’s multiple comparisons test (D), and "N-1" Chi-squared test (C and E); * p < 0.05 ** p < 0.01 *** p < 0.001; Exact p-values, number of fish per group and statistics can be found in S11 Table.

(A) Figures and graphs showing that disruption of tbx2a/b in zebrafish does not result in overall significantly reduced numbers of phox2bb-GFP+ cells, however it does induce ENS phenotypes in category II and III in a significant number of zebrafish larvae. (B) Figures and graphs showing that disruption of ufd1l in zebrafish results in significantly reduced numbers of phox2bb-GFP+ cells and increase in the percentage of fish presenting with an ENS phenotype. (C) Figures and graphs showing that disruption of mapk8a/b in zebrafish does not result in overall significantly reduced numbers of phox2bb-GFP+ cells, however it does induce ENS phenotypes in a significant number of zebrafish larvae. Statistical tests used: students t-test and "N-1" Chi-squared test; * p < 0.05 ** p < 0.01 *** p < 0.001; Exact p-values and statistics can be found in S11 Table.

(A) A second injection round targeting gnl1, using a lower dose of ret morpholino, confirmed gnl1 epistasis with Ret. Statistical tests used: "N-1" Chi-squared test; * p < 0.05 (B) Visual representation of the distribution of genetic predispositions over HSCR patient groups. In total 197 patients born between 1973 and 2018 were evaluated by a clinical geneticist in the department of Clinical Genetics, Erasmus Medical Center, Rotterdam. Of these, 114 did not have associated anomalies nor a known syndrome. 29 patients had a known HSCR related genetic syndrome, including Down syndrome (n = 18). 153 out of 197 patients were genetically evaluated for RET gene involvement and 21 had a pathogenic RET variant. (C) Pie charts showing the incidence of rare CNVs containing “ENS genes”, CCR “ENS genes” and coding variants in HSCR patients. (D) Graphical representation of a hypothetical model explaining the relative contributions of the risk scores in our 3 patient groups. Error bars represent standard deviation.