Combinatorial multiomic analysis from a pedigree of Sox10 Dom Hirschsprung mice implicates Dach1 as a modifier of Enteric Nervous System development

Combinatorial multiomic analysis from a pedigree of Sox10Dom Hirschsprung mice implicates Dach1 as a modifier of Enteric Nervous System development | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Combinatorial multiomic analysis from a pedigree of Sox10 Dom Hirschsprung mice implicates Dach1 as a modifier of Enteric Nervous System development View ORCID Profile Joseph T. Benthal , View ORCID Profile Justin A. Avila , View ORCID Profile Jeffrey R. Smith , View ORCID Profile E. Michelle Southard-Smith doi: https://doi.org/10.1101/2025.11.02.686123 Joseph T. Benthal 1 Program in Human Genetics, Vanderbilt University , Nashville, TN, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Joseph T. Benthal Justin A. Avila 2 Vanderbilt Brain Institute, Vanderbilt University , Nashville, TN, USA 3 Stowers Institute for Medical Research , Kansas City, MO, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Justin A. Avila Jeffrey R. Smith 4 Genetic Medicine, Vanderbilt University School of Medicine , Nashville, TN, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Jeffrey R. Smith E. Michelle Southard-Smith 1 Program in Human Genetics, Vanderbilt University , Nashville, TN, USA 2 Vanderbilt Brain Institute, Vanderbilt University , Nashville, TN, USA 4 Genetic Medicine, Vanderbilt University School of Medicine , Nashville, TN, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for E. Michelle Southard-Smith For correspondence: michelle.southard-smith{at}vanderbilt.edu Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF ABSTRACT Background Hirschsprung disease (HSCR) is characterized by absence of enteric ganglia (aganglionosis) along variable lengths of the distal intestine. This disorder results from deficient colonization of fetal intestine by enteric neural crest-derived cells (ENCDCs). HSCR exhibits complex, multifactorial inheritance with penetrance and severity varying widely even within a family. SOX10 is among causal genes that predispose to aganglionosis. Yet, how gene interactions influence severity of HSCR aganglionosis is not understood. We previously used an F 1 -intercross strategy to map genetic modifiers of aganglionosis in the Sox10 Dom HSCR mouse. Here we employ an extended pedigree mapping approach and complementary omics analyses of the developing Enteric Nervous System (ENS) to identify modifier loci that affect severity of HSCR aganglionosis. Results To identify loci that modify aganglionosis extent we undertook genome-wide association study (GWAS) of an extended pedigree of Sox10 Dom mice on a mixed C57BL/6J x C3HeBFeJLe-a/a background. GWAS uncovered genetic modifiers of aganglionosis severity in this cohort. Within each modifier interval, we prioritize candidate genes based on gene expression in the developing ENS, proximity to open chromatin regions in ENCDCs, and presence of conserved SOX10 binding motifs. This strategy identified known genes in ENS development as well as multiple novel genes. Dach1 emerged as a top priority gene for modifying ENCDCs migration and thus influencing aganglionosis severity. Conclusions This study identifies genome intervals and intrinsic genes that modify Sox10 Dom aganglionosis severity and that are candidate modifiers of human HSCR severity. AUTHOR SUMMARY Hirschsprung disease is a complex genetic neurodevelopmental disorder that causes loss of neurons in the distal bowel. The length of gut lacking neurons, called “aganglionosis”, in HSCR patients can vary widely even between affected siblings. Multiple genes are Mendelian causative for HSCR, but little is known about the gene interactions responsible for the notable variation in aganglionosis severity. In this study, we use a mouse model of HSCR to identify genomic regions, “modifiers”, associated with length of aganglionosis. Genes active within these regions are then identified in RNA-sequencing and open chromatin data from progenitor cells, which form the enteric nervous system. These omics approaches identify both known and novel genes that can affect enteric neuron development and may underly HSCR severity in patients. Dach1, already known for effects on neuronal progenitor proliferation and migration in other aspects of the nervous system, emerged as the top priority gene. The findings greatly expand the gene network that influences HSCR aganglionosis. INTRODUCTION The ENS is essential for normal gastrointestinal (GI) motility and analyses of mouse models have identified genes that are essential for ENS development. ENS neurons and glia that make up the myenteric and submucosal ganglia along the entire length of the intestine are formed by ENCDCs that colonize the fetal gut during development ( Lake and Heuckeroth 2013 ). Gut colonization begins at 9.5 days post coitus (dpc) in mice as neural crest cells invade the foregut and then migrate along the full length of the gut by 14.5dpc ( Lake and Heuckeroth 2013 ). The wavefront of migrating enteric neural crest-derived cells (ENCDCs) leaves behind cells that differentiate into neurons and glia. Disruption of initial ENCDC migration can produce GI motility disorders such as HSCR or Waardenburg-Shah syndrome while disrupted differentiation of ENCDCs is thought to contribute to chronic intestinal pseudo-obstruction ( Mallory et al. 1986 ; Pingault et al. 1998 ). The characteristic phenotype for HSCR is aganglionosis—lack of enteric ganglia—of the distal colon at varying lengths, caused by deficient colonization of ENCDCs ( Lake and Heuckeroth 2013 ). Multiple mouse models mimic human HSCR, including Sox10 Dom ( Lane and Liu 1984 ). The Sox10 Dom HSCR model also recapitulates the variability of aganglionosis length that remains a poorly understood characteristic of human HSCR ( Southard-Smith et al. 1999 ). Mouse models on controlled genetic backgrounds offer opportunity to identify contributing genes to the variability of HSCR aganglionosis. Genetic mapping studies in HSCR mouse models to find genes contributing to aganglionosis severity have had limited success due to large genomic intervals from standard crosses and lack of ENCDC expression data for candidate gene prioritization. Initial work to define how genetic background affects aganglionosis was performed by Cantrell et al. who compared lengths of intestinal aganglionosis in congenic C3HeB/FeJ and C57BL/6J. Sox10 Dom strains (2004). C57BL/6J. Sox10 Dom mice exhibit notably more severe aganglionosis than C3HeB/FeJ. Sox10 Dom mice ( Cantrell et al. 2004 ). Owens et al. mapped five broad genomic modifier intervals (ranging from 8 to 30 cM) associated with aganglionosis via genome-wide linkage scan of Sox10 Dom/+ F 2 progeny derived from C3HeB/FeJ and C57BL/6J congenic lines (2005). Further evidence that gene interactions can alter extent of aganglionosis was reported by Maka et al., who found that loss of Sox8 increases extent of aganglionosis in crosses with Sox10 lacZ knockout mice (2005). Multiple gene defects contributing to the severity of HSCR aganglionosis were also demonstrated for Ret +/- ; Ednrb S/S mutants (McCallion et a., 2003). Previously, the ability to prioritize genes that may influence severity of aganglionosis was limited by lack of gene expression profiles for migrating ENCDCs. Recent success in bulk and single cell sequencing of ENS progenitor populations offers new avenues to distill candidate genes in modifier intervals. Stavely and colleagues recently utilized bulk RNA-seq of ENCDCs at the migrating wavefront and the cells behind the wavefront to identify differentially regulated genes at 11.5dpc ( Stavely et al. 2023 ). In addition, Zhao and colleagues performed scRNA-seq on entire mouse gut over the time course of ENCDC migration (2022). This dataset allows detection of genes that are expressed either within the ENCDCs or in the surrounding gut environment, an important aspect, since interactions between progenitors and the fetal gut mesenchyme are known to influence extent of migration. Such datasets offer unique opportunities to prioritize candidate genes that might be functional among cells that contribute to aganglionosis phenotype. Prior genetics analyses of human HSCR have focused primarily on Mendelian causes. Relatively limited efforts have pursued the characteristic variability of HSCR disease penetrance and severity. Multiple case-control GWAS have identified variants associated with human HSCR ( Garcia-Barcelo et al. 2009 ; Kim et al. 2014 ; Tang et al. 2016 ; Fadista et al. 2018 ). These studies identified common genetic variation at RET , NRG1 , and SEMA3C/D loci that is significantly associated with HSCR. A study investigating additional phenotypic anomalies that can accompany HSCR implicated copy number variation at the SOX2 , MAPK10 , ZFHX1B , PHOX2B , and SEMA3A loci; severity of aganglionosis in these patients was not considered ( Jiang et al. 2011 ). To date, six human studies identified variants or genes that influence the penetrance of HSCR ( Bolk et al. 2000 ; Parisi et al. 2002 ; de Pontual et al. 2007 ; Garcia-Barcelo et al. 2009 ; Jiang et al. 2011 ; Li et al. 2023 ). To our knowledge, no human study has undertaken quantitative analysis to identify genetic variation modifying the extent of HSCR aganglionosis length, which is responsible for clinical severity. In contrast, mouse models of HSCR enable precise measure of aganglionosis length on controlled genetic backgrounds to dissect genetics underlying severity. Given the ability to control genetic background in mice and combine transcriptomic resources to prioritize candidate genes, we sought to refine Sox10 Dom/+ modifier intervals and to identify candidate genes that influence length of aganglonosis. We investigated a 10-generation pedigree of 830 heterozygous affected Sox10 Dom mice that yielded improved resolution, capitalizing upon the increased number of meiotic recombination events relative to standard F 2 intercrosses. The pedigree analysis accomplished this by reintroducing the greater-effect alleles of the B6 genetic background at each generation ( Cantrell et al. 2004 ). The Sox10 Dom line was maintained by crosses of affected heterozygous males with wild type (WT) F 1 B6C3Fe- a / a females. Pedigree individuals were genotyped with a linkage mapping SNP set appropriate for its genetic resolution. Genome-wide association in this pedigree improved the resolution of known loci, and identified novel loci that modify aganglionosis length, including some with sex-biased effects. Candidate genes from each modifier interval were prioritized based on gene expression among ENCDCs and fetal gut during ENCDC migration, chromatin accessibility in fetal enteric neuronal progenitors, and conserved SOX10 transcription factor binding motifs. From the hundreds of genes within modifier intervals, 19 met multiple genetic and omics criteria. Among these we observed both novel and known genes involved in ENS development including Ednrb , Nrg1 , and Phox2b . Phox2b and Mcm3 appeared to be sex-biased modifiers. A top candidate gene identified through the omics analyses was Dach1 , which has not previously been associated with aganglionosis, yet is expressed in migrating ENCDC, is flanked by SOX10 conserved binding motifs, and resides within accessible chromatin in ENS neuronal progenitors. The modifier genes identified from this analysis are likely to influence ENCDC migration or differentiation in the developing gut. RESULTS Genome-wide SNP analysis identifies Sox10 Dom aganglionosis modifiers While maintaining the B6C3Fe-a/a. Sox10 Dom/+ strain, we observed notable phenotype variability among Sox10 Dom/+ pups produced by iterative outcrosses to B6C3Fe- a / a wildtype females with some pups dying from severe aganglionosis at postnatal day (P)7 while others survived beyond a year of age. To understand genetic variation associated with severity of aganglionosis, we quantitatively assessed aganglionosis for 830 P7-10 B6C3Fe-a/a. Sox10 Dom/+ pups. Sox10 Dom/+ pups were distinguished from their WT littermates via characteristic white ventral spotting, white feet, and confirmed by direct genotyping of the mutation ( Cantrell et al. 2004 ; Owens et al. 2005 ). Extent of aganglionosis was assessed by whole-mount acetylcholinesterase staining collected over multiple pedigree generations ( Fig 1A ; Supp. Table 1 ) using established methods ( Cantrell et al. 2004 ; Owens et al. 2005 ). Because Sox10 Dom pedigree maintenance relied on breeding Sox10 Dom males to WT F 1 B6C3Fe- a / a female mice, alleles for greater severity/aganglionosis were reintroduced at each generation. Reduced severity/protective alleles were selected for by mating Sox10 Dom males that survive to sexual maturity. A standard mouse linkage mapping panel including 1449 SNPs of which 876 were informative in the pedigree was used to genotype each Sox10 Dom/+ pup ( Supp. Table 2 ). No significant difference between males and females were detected in length of intestine or length of aganglionosis (length of intestine t-test p=0.084; length of aganglionosis Wilcoxon p=0.64; Fig 1B,C ,D,E ). The aganglionosis phenotype in Sox10 Dom/+ pups was zero-inflated rather than normally-distributed ( Fig 1D ), replicating prior studies ( Cantrell et al. 2004 ; Owens et al. 2005 ). Download figure Open in new tab Fig 1. Phenotype distribution of Sox10 Dom pedigree mice exhibits no significant separation of phenotype by sex. A The breeding strategy used for the Sox10 Dom mouse extended pedigree used in this study. B,C Length of intestine split by sex displayed on histograms and a box plot shows normal distribution and no significant difference between males and females. D,E Length of intestinal aganglionosis split by sex on histograms and a box plot shows a non-normal, zero-skewed distribution with no significant difference between males and females. To identify genomic intervals associated with aganglionic length in the B6C3Fea. Sox10 Dom pedigree, we tested association between SNPs distinguishing genomic intervals originating from C3HeB/FeJLe-a/a (C3Fea) or C57BL/6J (B6) with the variable aganglionosis phenotype. This employed Genome-wide Efficient Mixed Model Analysis (GEMMA), which mitigates false positives that could otherwise arise due to relatedness of inbred mouse strains (Zhou et al 2012). Measured agangionic length was evaluated as a proportion of total intestinal length. Sex was included as a covariate (443 females, 387 males). We prioritized consideration of associated loci that A) were significant after multiple testing correction, B) had a logarithm of the odds (LOD) score of ≥3, or C) replicated observations of an independent F 1 -intercross study ( Owens et al. 2005 ). Even so, we also discuss additional nominally significant loci due the pleiotropic nature of the trait as shown by Owens et al., (2005) . In the main GEMMA analysis, an interval on chromosome 5 (LOD score = 3.9) was significantly associated with aganglionic length after false discovery (FDR) rate correction ( Fig 2A , Supp. Fig 1A , Supp. Table 3 ). The greatest significance was observed at rs13478309 (chr5:6714120 of mm10, p-Wald = 8.310292e-05); the B6 allele G was associated with increased aganglionosis. This SNP is ∼48 kb from the transcription start site of Phox2b (mm10; Perez et al. 2024 ). This finding is consistent with the well-known role of Phox2b in ENS development ( Pattyn et al. 1999 ). Additional nominally significant loci were present on multiple chromosomes illustrated in Fig 2A . Four of these are novel by comparison to the prior F 1 -intercross analysis ( Owens et al. 2005 ). Sox10 did not overlap with the loci on chromosome 15 ( Fig 2A , Supp. Table 3 ). Download figure Open in new tab Fig 2. GEMMA in the Sox10 Dom pedigree identifies significant regions associated with aganglionosis. A) Manhattan plot visualizing association analysis of quantitative aganglionic length, excluding unaffected Sox10 Dom mutation carriers. B) Analysis of quantitative aganglionic length as in A, with exclusion of chromosomes 5 (harboring Phox2B ), 15 (harboring Sox10 ), and X (with potential to accentuate alternative loci). C) Alternative analysis comparing Sox10 Dom mutation carriers that are unaffected to carriers with any affected aganglionic length (binary analysis). Blue lines on Manhattan plots indicate marginal significance (unadjusted P-Wald<0.05), dotted red lines indicate Bonferroni-adjusted significance (adjusted P-Wald<0.05), and the solid red line indicates false discovery rate significance (adjusted P-Wald<0.05). To uncover additive effects of variants on other chromosomes and to account for the zero-inflated aganglionosis phenotype distribution, we evaluated two alternative models. First, given Phox2b ’s known role in ENS development, the known issues that can arise when performing analyses on the X chromosome, and the previously established “leave one chromosome out” method, we opted to perform our analysis excluding the SNPs on chromosomes 5, X, and 15, where Phox2B and Sox10 reside, respectively ( Pattyn et al. 1999 ; Broman et al. 2006 ; Lippert et al. 2011 ; Yang et al. 2014 ). This approach can augment the statistical signal of loci on alternative chromosomes. Using this approach, we observed additional nominally significant modifier loci on chromosomes 7 and 10 ( Fig 2B , Supp. Fig 1B , Supp. Table 4 ). Second, many Sox10 Dom mice in our pedigree do not have detectable aganglionosis, which could obscure SNPs associated with length versus the presence of aganglionosis ( Fig 1C ). To account for this, we converted aganglionosis measurements to a pseudo case-control binary phenotype. We then compared Sox10 Dom mice who exhibited aganglionosis to mutation carriers without any detectable aganglionosis. Nominally significant genetic association with absence or presence of aganglionosis was observed at the previously detected loci as well as additional regions (chromosomes 1, 3, 5, 7, 9, 13, and 15) ( Fig 2C , Supp. Fig 1C , Supp. Table 5 ). These pedigree-based association scan results expanded upon the previously published F 1 -intercross loci modifying Sox10 Dom aganglionosis ( Owens et al. 2005 ). The previously observed modifiers on chromosomes 5, 8, 11, and 14 were again detected ( Table 1 ; Owens et al. 2005 ). Additionally, the pedigree analysis detected multiple novel modifiers on chromosomes 1, 2, 3, 13, and 19 that were nominally significant ( Table 1 ). The chromosome 3 loci observed in this pedigree-based analysis were distinct from that of the intercross ( Table 1 ; Owens et al. 2005 ). Directions of effect (DOE) for loci replicating in the prior F 2 -intercross study were concordant ( Supp. Fig 2A,B,F-H , Supp. Tables 1,2 , see Owens et al. 2005 Fig 4 ). However, several of the loci appeared to have sex-biased effects on aganglionosis in the pedigree, which had not previously been observed in the F 1 -intercross ( Supp. Fig 2A-E,I-N , Supp. Tables 1,2 ). Given this sex-biased observation, sex-stratified analysis was conducted. View this table: View inline View popup Table 1. Sox10 Dom aganglionosis modifiers compared to Owens et al. View this table: View inline View popup Table 2. Sex-biased Sox10 Dom aganglionosis modifier locations. View this table: View inline View popup Table 3. Top 2 GEMMA peak SNPs per Sox10 Dom aganglionosis modifier interval definition. Assessing potential sex-bias among modifiers of Sox10 Dom aganglionosis Because some of the SNP variants in our pedigree associations exhibited marginally significant sex-biased DOE ( Supp. Fig 2A-E,I-N , see methods for shorthand terms used for each genome-wide analysis, Supp. Tables 1,2 ), we conducted sex-stratified genome-wide analyses for both quantitative and binary aganglionosis phenotypes using GEMMA (Zhou et al. 2012). The male-specific genome-wide scan replicated the Phox2b on chromosome 5, which was again genome-wide significant ( Fig 3A , Supp. Fig 1D , Supp. Table 6 ). This locus was notably attenuated in the female-specific scan. Several other loci were nominally significant in the male-specific scan, including on chromosome 14 near Ednrb , a critical gene for ENS development ( Fig 3A , Supp. Table 6 ).The female-specific quantitative scan only yielded nominally significant loci that did not remain significant upon multiple testing correction ( Fig 3B , Supp. Fig 1E, Supp. Table 7). Download figure Open in new tab Fig 3. Sex-specific genome-wide association for aganglionosis modifiers in the Sox10 Dom pedigree. Manhattan plots visualizing GEMMA genome-wide scan results of male- (A) and female-specific (B) association analysis of quantitative aganglionosis length, excluding unaffected Sox10 Dom carriers. Manhattan plots visualizing GEMMA genome-wide scan results of male- (C) and female-specific (D) analysis comparing Sox10 Dom mice that are affected by any length of aganglionosis to those mice that are unaffected (binary analysis). Blue lines on Manhattan plots indicate marginal significance (unadjusted P-Wald<0.05) and dotted red lines indicate Bonferroni-adjusted significance (adjusted P-Wald<0.05). Download figure Open in new tab Fig 4. Differential gene expression in the developing gut within aganglionosis modifier intervals defines a first list of candidate genes. A Pipeline for prioritizing Sox10 Dom aganglionosis modifier interval genes by differential gene expression. B UMAP of reprocessed scRNA-seq data from Zhao et al. 2022 split by developmental timepoint and colored by their cell type definitions. The red asterisk indicates the location of the enteric neural crest-derived cells (ENCDCs). C Dot plot showing expression of differentially expressed genes near the top two LOD peaks per GEMMA GWAS (out of the nearest 10 genes per top 2 LOD peaks, those that are differentially expressed). Cell types have been consolidated to those that are known to be relevant for migrating enteric neural crest-derived cells and those with shared names (numbered) are combined to one identifier. Dots to the left of each gene correspond to which GEMMA genome-wide scan each gene originates. The sex-specific binary association analysis only identified nominally significant loci. Male-specific binary phenotype analysis loci mirrored those of the quantitative analysis, but were less prominent ( Fig 3C , Supp. Fig 1G and F, Supp. Table 8 ). Female-specific binary phenotype analysis also yielded only nominally significant variants ( Fig 3D , Supp. Table 9 ). Altogether, these findings suggest that some Sox10 Dom modifiers of aganglionosis severity may exert sex-biased effects ( Tables 2,3 ). Prioritization of candidate genes based on expression in the fetal mouse intestine To identify relevant candidate genes potentially underlying the variability of Sox10 Dom aganglionosis, we first assessed which genes within modifier intervals are expressed in ENS precursor cells and the fetal gut environment during ENCDC migratory stages. Genes were detected within modifier intervals as defined in Methods. We initially leveraged scRNA-seq data from total mouse fetal gut collected over 9.5 to 15.5dpc that profiled ENCDCs and surrounding gut mesenchyme ( Fig 4A,B ; Zhao et al. 2022 ). After reprocessing this scRNA-seq data, we filtered for genes within modifier intervals that have log-normalized pseudobulk expression greater than 1 (see Methods; Table 4 , Supp. Table 10 ; Karolchik et al. 2004 ; Zhao et al. 2022 ; Hao et al. 2024 ). Differential gene expression analysis for the expressed genes by age and cell type determined which genes are expressed in specific cell types and timepoints in the developing mouse gut ( Table 4 , Supp. Table 11 ). We detected differential expression of Phox2b and Ednrb , genes that are highly expressed in ENCDCs at all timepoints and which are well known genes that influence ENS development ( Fig 4C , Supp. Fig 3C ; Cantrell et al. 2004 ; Elworthy et al. 2005 ; Hakami et al. 2006 ; Stanchina et al. 2006 ; Nagashimada et al. 2012 ). The top female associated SNP had a LOD score of >3, while the second most significant female SNP had a higher LOD score than other modifiers that were replicated from the F 2 study ( Owens et al. 2005 ). Therefore, we further examined differentially expressed genes nearest the top two female associated SNPs including Mcm3 , Il17f , and Isl1 . Mcm3 and Isl1 are expressed in ENCDCs ( Fig 4C , Supp. Fig 3C ). View this table: View inline View popup Table 4. Number of genes within each modifier interval set prioritized to those that are expressed and differentially expressed in scRNA-seq of the fetal gut. Expression of genes in the gut environment influencing migration of ENCDCs is consistent with known roles for environmental factors on migration of neural crest progenitors ( Nagy and Goldstein 2017 ). To evaluate modifier interval candidate gene gut mesenchyme expression and potential to affect communication with ENCDCs during bowel colonization, we used CellChat to estimate ligand-receptor communication in the scRNA-seq dataset between ENCDCs and other cell types ( Supp. Fig 4A, Supp Table 12 ; Jin et al. 2021 ; Zhao et al. 2022 ; Jin et al. 2024 ). We observe five genes within modifier intervals that CellChat identified as part of communication between ENCDCs and other gut cells during migration along the bowel including Ppia , Col1a1 , Lgals9 , Ednrb , and Ngr1 ( Supp. Fig 4A, Supp Table 12 ). While none of these genes are within modifier intervals that are significant after multiple testing correction, they fall within intervals that were detected in the F 1 -intercross study ( Owens et al. 2005 ). Expression of these CellChat-identified modifier interval genes and their binding partners shows at least one of each are expressed in ENCDCs in the scRNA-seq dataset ( Supp. Fig 4B ). Interestingly, CellChat identified predicted enriched signaling of the Nrg1 - Erbb3 signaling pathway at 13.5dpc ( Supp. Fig 4A,B ) when the hindgut is being colonized. Nrg1 expression is present in the developing ENS and is a known risk gene for HSCR ( Maka et al. 2005 ; Jiang et al. 2017 ; Memic et al. 2018 ). Nrg1 ’s signaling partner Erbb3 is highly expressed in ENCDCs ( Supp. Fig 4B ). Taken together, these enriched predicted signaling pathway modifier interval genes add another prioritization layer based on the potential to influence ENCDC-fetal gut environment communication during bowel colonization. Differential expression of modifier interval genes in the migrating wavefront of enteric neural crest-derived migrating cells further prioritizes candidate genes We further examined expression of genes within modifier intervals utilizing bulk RNA-seq from the migrating wavefront of ENCDCs during intestinal colonization at 11.5dpc ( Stavely et al. 2023 ). Fetal intestinal colonization by ENCDCs relies heavily on migratory capacity of the leading edge of the ENCDC wavefront, exhibiting distinct transcriptional features ( Young et al. 2002 ; Young et al. 2004 ; Simpson et al. 2007 ; Young et al. 2014 ; Zhang et al. 2019 ; Stavely et al. 2023 ). At 11.5dpc, ENCDCs are transitioning around and across the midgut fold to colonize the hindgut. Differentially regulated genes at this stage are likely relevant for ENCDC migration deficits in the Sox10 Dom mouse gut. We therefore prioritized genes within Sox10 Dom aganglionosis modifier intervals that were differentially expressed in the ENCDC wavefront versus the lagging cells ( Fig 5A , Table 5 ). Of these, five genes were downregulated ( St18 , Hoxb5 , Prph , Slc10a4 , and Uchl1 ) and 9 genes were upregulated ( Myh10 , Col1a1 , Alkbh5 , Pdgfra , Ptprg , Sh2b3 , Dach1 , Hs3st3b1 , and Tgfbi ) in the migrating wavefront ( Fig 5B , Table 5 , Supp. Table 13 ). Several of these genes including Hoxb5 , Pdgfra , Ptprg , Alkbh5 , Dach1 , and Tgfbi have been implicated in various aspects of neural crest development ( Fu et al. 2003 ; Lui et al. 2008 ; Wang et al. 2013 ; Zhu et al. 2014 ; Kim et al. 2020 ; Mo et al. 2020 ; Wang et al. 2021 ; Liu et al. 2024 ). Download figure Open in new tab Fig 5. Overlap between differentially expressed genes at the enteric neural crest-derived migrating wavefront and genes within modifier intervals identifies migratory candidate genes. A Venn diagram visualizing pipeline for prioritizing Sox10 Dom aganglionosis modifier interval genes by differential gene expression at the migrating wavefront. Red asterisk indicates overlap between modifier intervals and upregulated genes in the migrating wavefront. B Volcano plot of all differentially expressed genes of the migrating wavefront enteric neural crest-derived cells at 11.5dpc versus the lagging/stationary cells. Red and blue dots indicate genes upregulated and downregulated in the wavefront, respectively. Genes labeled are those within modifier intervals. C UMAP of scRNA-seq from Zhao et al. 2022 enteric neural crest-derived cells (ENCDCs) based on source data annotations isolated from the main dataset and split by timepoint. D Expression of Sox10 , Elavl4 , and Phox2a to label progenitor and neuronal cells, respectively. E Feature plots of Zhao et al., 2022 ENCDCs colored by expression of genes differentially downregulated (left) and genes differentially upregulated (right) in the migrating wavefront. UMAPs visualizing expression genes differentially downregulated in the wavefront are not split by time point while UMAPs for genes differentially upregulated are split by time. Dots next to gene names indicate from which modifier interval set each gene falls within. View this table: View inline View popup Table 5. Genes differentially expressed in the enteric neural crest-derived cell wavefront within Sox10 Dom aganglionosis modifier intervals. To validate expression in a separate dataset and determine if these differentially expressed genes at the migrating ENCDC wavefront are expressed in specific groups of cells within the ENCDC population, we subset the whole gut developmental scRNA-seq dataset to focus just on ENCDCs ( Fig 5B,C ; Stavely et al. 2023 ). Expression of Sox10 and Phox2a identify progenitors and neuronal cells, respectively ( Fig 5D ). We expected modifier candidate genes upregulated in the migrating wavefront of the ENCDCs might be expressed in progenitor cells which are highly migratory, while downregulated genes would be expressed in more mature neuronal ENCDC cells ( Fig 5B,D ). Downregulated genes in the ENCDC wavefront are mostly expressed in more mature, neuronal populations ( Fig 5E —Lagging ). However, wavefront upregulated genes are expressed in either a specific intermediate population ( Col1a1 , Tgfbi , Pdgfra ) or increase in expression over time ( Ptprg , Dach1 ; Fig 5E —Migrating Wavefront ). ScRNA-seq during migration validates expression of these candidate genes and provides insight into ENCDC type-and temporally-specific expression. Evolutionarily conserved SOX10 binding sites within Sox10 Dom modifier intervals overlap or are near genes differentially expressed in the migrating wavefront To further distill relevant candidate genes within Sox10 Dom aganglionosis modifier intervals potentially regulated by SOX10 DNA binding, we leveraged a dataset of dimeric SOX10 binding motifs conserved across chick, mouse, and human genomes ( Gopinath et al. 2016 ). We located 24 conserved SOX10 binding motifs within modifier intervals ( Table 6 , Supp. Table 14 ). We located the closest gene to each conserved SOX10 binding motif, yielding 17 unique genes, some of which are nearby multiple SOX10 binding motifs ( Table 6 ). Most of these binding motifs were intronic to the annotated gene with only one motif positioned outside the nearest gene (∼1.3kb away from Sox2 ; Table 6 ). Each 1Mb region surrounding the modifier interval-contained SOX10 binding motifs was manually examined on the mm10 genome to determine whether there were nearby genes that were already identified via other data modalities or already known to affect ENS development ( Perez et al. 2024 ). This added 6 genes to this list ( Table 6 ). Several of these genes have already been identified in this work by alternate means, including Mcm3 , St18 , Pdgfra , Ptprg , and Dach1 ( Table 6 ). Also identified are several genes that have two or more intronic or close conserved SOX10 binding motifs including Dach1 , Bcas3 (close to Tbx2 ), Tenm2 , Ranbp17 (close to Tlx3 ), and Mecom with four intronic motifs ( Table 6 ). Lastly, other genes not previously identified through other omics were near or overlap with conserved SOX10 binding motifs within modifier intervals, including Bai3 ( Adgrb3 ), Cyp7b1 , Nlg1 , Slit2 , Ppargc1a , Prdm8 (close to Antxr2 ), Fam204a , and Sox2 ( Table 6 ). Several of these genes are expressed in the ENCDCs ( Fig 5 , Supp. Fig 5 ; Zhao et al. 2022 ). Altogether, these 24 conserved SOX10 binding motifs within modifier intervals highlight potential targets of SOX10 binding that could influence the aganglionosis phenotype caused by the Sox10 Dom mutation. View this table: View inline View popup Table 6. Conserved SOX10 TF binding sites within Sox10 Dom aganglionosis modifier intervals. Chromatin accessibility within enteric neuronal progenitors highlights putative regulatory regions within aganglionosis modifier intervals Modifiers can consist of variants in cis -regulatory elements that affect gene expression ( Shih and Fay 2021 ; Kalra et al. 2024 ). To assay open chromatin, that might harbor sites of gene regulatory activity within Sox10 Dom modifier intervals, we utilized single nucleus assay for transposase-accessible chromatin with sequencing (snATAC-seq) within mouse fetal ENS cells ( Supp. Fig 6A ). We captured ENCDCs at 16.5dpc based on expression of a Phox2b H2B-CFP transgene shown to mark glial, neuronal progenitor, and neuronal cells ( Corpening et al. 2008 ; May-Zhang et al. 2021 ). By 16.5dpc, the gut has been fully colonized, yet enteric neurogenesis is ongoing. We processed 13431 Phox2b H2B-CFP+ nuclei and performed unsupervised clustering to obtain 17 clusters ( Fig 6A ; Korsunsky et al. 2019 ; Stuart et al. 2021 ; Hao et al. 2024 ). We then utilized scRNA-seq of WT 15.5dpc ENS to annotate snATAC-seq nuclei via estimated gene expression (see Supplementary Methods; Supp. Fig 6B ), assigning six main groups (progenitors, neuroblast1 and 2, and neuronal branches A, B, and C; Fig 6B , Supp. Fig 6C ; Avila et al. 2025 ). Differential chromatin accessibility analysis was used to localize 71521 differentially accessible chromatin (DAC) regions for each of these six main groups ( Supp. Table 15 ). Annotations did not deviate across groups except for Neuroblast2, which we postulate could be due to fewer nuclei assigned as compared to the other groups (n = 102 Neuroblast2 nuclei / 13431 total nuclei; Supp. Fig 6D , top panel). We filtered these DAC to those within Sox10 Dom modifier intervals, yielding 9151 unique regions ( Supp. Table 16 ). Modifier interval DAC regions in the Neuroblast2 nuclei are limited to promoter and distal intergenic categories, and this might be, again, due to fewer nuclei assigned to this group ( Supp. Fig 6D , bottom panel). Download figure Open in new tab Fig 6. Whole gut 16.5dpc Phox2b H2B-CFP+ snATAC-seq-derived transcription factor binding motif enrichment in differentially accessible chromatin identifies modifier interval-contained enriched TF motifs. A UMAP of snATAC-seq nuclei colored by unsupervised clusters. B UMAP displaying the result of label transfer of cell types from scRNA-seq to snATAC-seq. C Scatterplot of enriched TF binding motifs from the modifier interval-contained differentially accessible chromatin from progenitor and neuroblast clusters filtered by expression of their corresponding genes in the fetal gut ENS. The scatterplot points are colored by which modifier interval set each TF motif is enriched, with duplicates resulting from the analysis being performed by specific interval (“assay”). The percentage difference between the number of TF motifs in the DAC and the background open chromatin (“percent.diff”) is indicated by the size of the point. To determine whether conserved SOX10 binding sites from Gopinath and colleagues localized to DAC regions within modifier intervals, we examined the overlap between these datasets ( Gopinath et al. 2016 ). Two of the conserved SOX10 binding motifs overlapped with modifier interval DAC. These were intronic to Ranbp17 (close to Tlx3 ; chr11:33400347-33400362) and Dach1 (chr14:97891443-97891462; Table 6 , Supp. Table 17 ). Transcription factor (TF) binding is important for the transcription of genes, and prominent candidate genes identified in this study are TFs including Phox2b , that have established roles in enteric neurogenesis ( Pattyn et al. 1999 ). We performed TF binding motif enrichment on modifier interval-contained DAC from progenitor and neuroblast groups to evaluate potential regulators of ENS development in Sox10 Dom modifier intervals. We focused on progenitors and neuroblasts because these are most similar to migrating ENCDCs. This approach examined both open chromatin in progenitors and neuroblasts as well as regions that are closed in differentiating enteric neuronal cells. This resulted in 797 enriched TF binding motifs in modifier interval DAC compared to the rest of accessible chromatin ( Supp. Table 18 ). We filtered these TF binding motifs to select for those whose encoding genes are expressed in the ENCDC scRNA-seq data, revealing 367 motifs ( Supp. Table 19 ). This enrichment analysis coupled with expression identifies the following TFs from modifier intervals: Phox2b , Rbpj , Hoxb2/3/4/5/9 , Sox2 , Pou5f1 , Tlx2 , Tbx3 , Pou3f1 , and Mecom ( Fig 6C , Supp. Table 19 ). Several of these TFs participate in ENS development ( Borghini et al. 2007 ; Taylor et al. 2007 ; López et al. 2018 ; Memic et al. 2018 ; Wright et al. 2021 ; Avila et al. 2025 ). Prioritization of top aganglionosis modifier candidate genes across data modalities To rank candidate genes within Sox10 Dom aganglionosis modifier intervals ( Supp. Table 20 ) we assessed both genetic and multiomics data in a prioritization pipeline ( Fig 7 ). In this process, we applied two scoring metrics: a modifier prioritization score and an omics evidence score. For our modifier interval prioritization score, we assigned 0 to 3.5 points to each gene within modifier intervals based on the following criteria to minimize false positives: 1) peak SNPs of modifier intervals significant after multiple testing correction were assigned one point, 2) peak SNPs of modifier intervals with LOD score of ≥3 were assigned one point, and 3) intervals replicated from the prior F 1 -intercross study were assigned 1.5 points ( Table 7 , Supp. Table 21 ; Owens et al. 2005 ). We defined the preponderance of evidence omics score over a range of 0 to 5 corresponding to the number of separate lines of evidence implicating candidate genes. This incorporated single cell differential expression, estimated cell communication (CellChat), differential expression at the ENCDC wavefront, snATAC-seq-derived DAC regions in ENCDCs, ENCDC-expressed genes whose TF binding motifs are enriched within DAC regions in modifier intervals, and proximity to evolutionarily conserved SOX10 binding motifs ( Fig 7 , Supp. Table 22 ). Genes from each modifier interval were assessed independently through our prioritization pipeline to identify relevant candidates with potential to influence aganglionosis ( Fig 7 , Supp. Table 22 ). Thirty genes emerged that appear in three or more of our omics analyses. Of these, 11 with modifier interval prioritization scores of 0, indicating lack of strong genetic evidence, including Mecom and Sox2 , were excluded from further analysis. Nineteen genes exhibited modifier prioritization scores >0 and multiple lines of omics evidence ( Table 7 ). The top candidate modifier gene, Dach1 , exhibits five lives of omics evidence. Three additional genes have four supporting lines of evidence ( Col1a1 , Pdgfra , and Tbx2 ; Table 7 , Supp. Table 22 ). Several genes that emerged from the pipeline are already known to be involved in ENS development and HSCR including Phox2b , Ednrb , and Nrg1 , which effectively serve as internal controls. Download figure Open in new tab Fig 7. Analysis pipeline used to identify candidate genes from Sox10 Dom pedigree-derived aganglionosis modifier intervals. Outline for candidate gene analysis pipeline following GEMMA genome-wide scans to identify genomic regions that modify the aganglionosis phenotype in Sox10 Dom mice, then use bulk and single cell sequencing strategies to identify candidate causal genes. After identifying candidate genes, these are then used to find gene body variants in C3Fe as compared to mm39 (C57BL/6J), then compared to human phenotype-genotype association data. DGE: differential gene expression. View this table: View inline View popup Table 7. Cumulative omics evidence supporting priority genes within Sox10 Dom aganglionosis modifier intervals. After comparison across data modalities to identify candidate genes with the highest relevance scores, we probed each top candidate gene sequence for C3Fe variants, as compared to B6, with predicted high impact on expression ( McLaren et al. 2016 ). Of these candidate genes, Phox2b and Col1a1 contained C3Fe variants with predicted functional impact ( Table 8 , Supp. Table 23 ). Phox2b contained a C3Fe variant (T>TC insertion, chr5:67255011, transcript ENSMUST00000174251.2 UCSC Genome Browser, mm39) located at a splice junction boundary that could affect splicing efficiency and overall gene expression levels ( Perez et al. 2024 ). Col1a1 also contained a C3Fe variant predicted to impact splicing ( Table 8 , Supp. Table 23 ). Identification of sequence variants between the C3Fe and B6 strains for these candidate genes adds further evidence to the likelihood that strain differences may exacerbate severity of Sox10 Dom aganglionosis. View this table: View inline View popup Table 8. Relating candidate genes to genetic variants, related phenotypes in mouse and human. Candidate modifiers of aganglionosis related to human GI phenotypes To relate our Sox10 Dom aganglionosis modifier candidate genes to human GI and HSCR disease loci, we assessed whether these genes’ human orthologs were in LD with SNPs identified as significant from four human HSCR GWAS and GWAS of stool frequency ( Garcia-Barcelo et al. 2009 ; Kim et al. 2014 ; Tang et al. 2016 ; Fadista et al. 2018 ; Bonfiglio et al. 2021 ). Two separate HSCR GWAS summary statistics had significant SNPs that were in LD with genes identified in this study, including Phox2b , Col1a1 , Rbpj , Nrg1 , Ednrb , Hoxb2/5/13 , and Uchl1 ( Table 8 , Supp. Tables 24,25,26 ; Garcia-Barcelo et al. 2009 ; Tang et al. 2016 ). The GWAS of stool frequency had significant SNPs that were in LD for Col1a1 , Rbpj , Nrg1 , and Elf1 ( Supp. Table 27 ; Bonfiglio et al. 2021 ). While not in an LD block with significant SNPs, DACH1 is ∼2.6Mb away from an associated SNP in the stool frequency GWAS ( Bonfiglio et al. 2021 ). To further relate our candidate Sox10 Dom aganglionosis modifier genes to other human GI, neural, and neural crest phenotypes, we probed the publicly available exome-based PheWAS sever Genebass. We evaluated whether variants in human orthologs of our candidate mouse genes might be associated with ENS-related traits. We utilized a text/term-based search among the top 20 phenotypes that were returned for each of the orthologs (Supplemental Methods; Table 8 ; Karczewski et al. 2022 ; Genebass 2024). Of the orthologs examined in PheWAS, each had phenotypes based on partial string searches that captured terms based on neural crest, nervous system, or GI (Supplemental Methods). Among these is DACH1 with four different phenotypes relating to GI surgery ( Table 8 , Supp. Table 28 ). Interestingly, among the terms found for DACH1 is “Anterior resection of rectum and exteriorisation of bowel”, a process often used to treat HSCR. DISCUSSION While HSCR patient studies have largely focused on causal genes, animal models, such as Sox10 Dom , are advantageous for defining the role of genetic background on aganglionosis severity. Prior Sox10 Dom modifier mapping studies identified very large genomic intervals and relatively few individual genes emerged as modifiers of aganglionosis severity and penetrance ( Owens et al. 2005 ). Omics from relevant cell and tissue types can help prioritize large numbers of genes in modifier intervals to select candidates with greatest potential to influence aganglionosis. In this study, we conducted genome-wide analyses on an extended pedigree of the Sox10 Dom mouse model of HSCR to refine aganglionosis modifier intervals and identify candidate modifier genes. We identified multiple genomic intervals that are associated with severity of Sox10 Dom aganglionosis, some of which replicate intervals identified in the prior F 2 study, and others that are novel regions ( Owens et al. 2005 ). 7 Several nominally significant modifier intervals are associated in a sex-bias manner. Utilizing multiple omics datasets including scRNA-seq, bulk RNA-seq, snATAC-seq, and TF binding motifs, we analyzed relevance of modifier interval candidate genes. This strategy revealed multiple highly relevant candidate modifier genes that are expressed in developing gut and which either have known roles in ENS development or exert effects on other aspects of cell migration. Our analysis confirms Phox2b is a candidate modifier of Sox10 Dom aganglionosis and identifies multiple novel candidate genes likely to influence aganglionosis severity. Among these, Dach1 stands out as the highest priority candidate based on omics data and prior evidence of this gene affecting cell migration ( Davis et al. 2001 ; Wu et al. 2008 ; Kim et al. 2020 ). Quantitative trait mapping in mice has blossomed with the availability of mouse genetic resources like the collaborative cross ( Lorè et al. 2020 ; Dorman et al. 2021 ). However, challenges remain for identifying disease modifiers that rely on inclusion of a mutant allele. We took advantage of a multi-generational Sox10 Dom mouse pedigree for GEMMA genome-wide analyses. This approach compliments our previous F 1 -intercross strategy to map Sox10 Dom aganglionosis modifiers by taking advantage of additional recombination in pedigree generations. Aganglionosis modifiers identified here replicate intervals located in the F 1 -intercross including a genome-wide significant modifier on chromosome 5 and marginally significant modifiers on chromosomes 3 (shifted), 8, 11, and 14 ( Owens et al. 2005 ). We also detect novel marginally significant modifiers, including those that are sex-biased ( Tables 1,2 ). Our dual score approach—confidence of modifiers and lines of evidence through omics analysis— allowed us to prioritize candidate genes from large genomic intervals. Among the total 6216 annotated genes in modifier intervals, 683 were expressed in either ENCDCs or the surrounding gut mesenchyme through which these cells migrate ( Zhao et al. 2022 ). Subsequent gene filtering identified 5 cell communication candidate genes, 26 candidates near conserved SOX10 binding motifs, 9151 DAC regions enriched for motifs of 367 ENCDC-expressed TFs from snATAC-seq, and 14 genes differentially expressed at the ENCDC migrating wavefront. Thirty candidate genes with cumulative evidence aggregated from multiple modalities emerged as putative modifiers of Sox10 Dom aganglionosis, however only 19 of these also exhibited significant modifier scores >0 ( Table 7 ). Of the 19 candidate genes, Phox2b , Rbpj , and Uchl1 were in modifiers that 1) had peak SNPs that were significant past the threshold for multiple comparisons, 2) had peak SNPs that had LOD scores ≥3, and 3) overlapped with intervals derived from the prior F 1 -intercross study ( Table 7 , Supp. Tables 21,22 ; Owens et al. 2005 ). Pdgfra and Antxr2 met two of these criteria ( Table 7 , Supp. Tables 21,22 ). Dach1 was within a replicated modifier interval ( Table 7 , Supp. Tables 21,22 ; Owens et al. 2005 ). We then examined these candidate genes for variants of predicted high impact in the C3Fe genome versus C57BL/6J. This analysis revealed potential function-altering C3Fe variants for Phox2b and Col1a1 that could affect splicing ( McLaren et al. 2016 ). The mouse aganglionosis modifier genes identified here may have relevance as modifiers of human GI disease. Therefore, we compared the mouse Sox10 Dom aganglionosis candidate genes to significant SNPs identified from HSCR and stool frequency GWAS ( Garcia-Barcelo et al. 2009 ; Kim et al. 2014 ; Tang et al. 2016 ; Fadista et al. 2018 ; Bonfiglio et al. 2021 ). Of the candidate genes, nine of these genes were in LD blocks with significant HSCR GWAS SNPs and four were in LD blocks with significant stool frequency GWAS SNPs ( Table 8 ). In addition, we use exome-based PheWAS study results from the UKBiobank (Genebass) to find GI, neuronal, and neural crest phenotypes associated with variants in our candidate genes ( Karczewski et al. 2022 ; Genebass 2024). These comparisons suggest that our candidate genes in Sox10 Dom aganglionosis modifier intervals are likely relevant for human HSCR phenotype variation or other GI phenotypes. In this study, we generated and mined a novel snATAC-seq dataset of 16.5dpc fetal ENS collected from Phox2b H2B-CFP mice ( Corpening et al. 2008 ; May-Zhang et al. 2021 ). We used this dataset to prioritize modifier interval candidate genes via differentially accessible chromatin and enriched TF binding motifs. Because this dataset contains chromatin accessibility for populations that include progenitors, neuroblasts, and maturing neuronal linages, it has the potential to be highly informative for researchers interested in regulatory regions that control neuronal diversification and maturation in the fetal ENS. This dataset can be used further for analyses of TF binding motif accessibility changes over neurogenesis for each neuronal trajectory relative to analogous scRNA-seq datasets. There is tremendous temporal and spatial complexity in ENCDC colonization that hinges upon migration, proliferation, and differentiation. Given that HSCR aganglionosis is a consequence of failed colonization, we hypothesize these top modifier candidates regulate different aspects of ENCDC migration ( Fig 8 ). There is substantial prior evidence of Dach1 influencing migration in neural crest, embryonic fibroblasts, and breast cancer ( Davis et al. 2001 ; Wu et al. 2008 ; Kim et al. 2020 ). Specifically, Dach1 was associated with early neural crest migration in Xenopus ( Kim et al. 2020 ). In contrast Col1a1 in zebrafish is implicated downstream of retinoic acid (RA) signaling, which has been shown to be important for ENS development and ENCDC migration ( Sato and Heuckeroth 2008 ; Li et al. 2010 ; Li et al. 2018 ; Uribe et al. 2018 ; Gao et al. 2021 ). Pdgfra affects both migration and survival in oligodendrocytes as regulated by Sox9 and Sox10 while PDGFR signaling has been shown to affect nitrergic enteric neuron specification in vitro ( Finzsch et al. 2008 ; Majd et al. 2025 ). In zebrafish, tbx2a/b CRISPR/Cas9 targeting caused reduced ENS cell density, suggesting function for differentiation or proliferation of ENCDCs ( Kuil et al. 2021 ). Two genes that fell off our omics candidate list due to lack of strong genetic evidence include Sox2 and Mecom . However, these two genes exhibit four lines of omics evidence. Mecom affects neural crest-derived chondrocyte differentiation, orientation, and polarity in mouse and zebrafish models ( Shull et al. 2022 ). Sox2 has been shown to contribute to maintenance of progenitor states of neuronal progenitors and has been identified as a modifier of human HSCR via copy number variation ( Jiang et al. 2011 ; Wakamatsu et al. 2021). These genes may exert small additive effects on Sox10 Dom aganglionosis despite their failure to reach significance in our genetic analysis. The depth of prior analyses assist in understanding how these different genes might contribute to the aganglionosis phenotype in Sox10 Dom mice. Download figure Open in new tab Fig 8. Hypotheses for candidate Sox10 Dom modifier interval genes to alter ENCDC migration. Cartoons displaying migration of enteric neural crest-derived cells (ENCDCs), differentiation of ENCDCs towards more mature ENS fates, and proliferation of ENCDCs at the migrating wavefront. In blue are the typical functions, and in red are cells and functions that are the hypothesized results of perturbation of the candidate genes. Candidate genes hypothesized to be involved in each function are in smaller boxes in the upper right of each function box, with superscript sex symbols for those genes that are within sex-specific modifier intervals. There are several limitations to our study. The omics data analyzed in this study are sourced from the fetal gut and therefore does not include genes whose expression during initial migration of vagal neural crest from neural tube to the foregut, which also could influence Sox10 Dom aganglionosis. Ideally future analysis will capture expression profiles of early migrating vagal neural crest. There is also inherent stochasticity of ENCDC development even in congenic inbred backgrounds which can contribute to variation in aganglionosis ( Kaern et al. 2005 ). Resolution for detecting modifiers on chrX is limited both due to the cross structure and the analysis tools we applied. GEMMA does not include algorithms that can account for copy number or mosaicism that are needed for analysis of chrX SNPs (Zhou et al., 2012). We expect since only male Sox10 Dom mice are bred to propagate the pedigree, and WT B6.C3Fe F 1 female mice are introduced into the pedigree at each generation, there is minimal recombination on chrX. Lastly, future work validating the effects of individual modifier genes identified here will be required. Analysis of individual modifier gene loss on ENCDP migration or the extent of aganglionosis will require complementation test crosses with Sox10 Dom mutants or knockdown studies in homogenous genetic backgrounds. In conclusion, we demonstrate the ability to improve resolution of Sox10 Dom/+ aganglionosis modifier intervals from mapping in an extended pedigree. Combinatorial omics analysis identifies a high priority list of modifier genes that meet multiple genetic and genomic metrics. Detection of Phox2b as a well-known gene in ENS development confirms the pipeline approach. Dach1 emerges as a novel modifier of Sox10 Dom aganglionosis. Cumulative evidence for Dach1 includes genetic interval replication across studies, prior reports of Dach1 mutation on neural crest migration, and omics evidence of SOX10 binding sites in accessible chromatin around this gene. Our analysis enriches the gene network that impacts ENS development and reveals important genes for future variant analyses in human HSCR and other GI motility disorders. METHODS Mouse husbandry All experimental protocols were approved by the Institutional Animal Care and Use Committee (IACUC) at Vanderbilt University Medical Center. All mice were housed in a modified barrier facility on a 14-hour on, 10-hour off-light cycle in high-density caging (Lab Products Inc., #10025) with breeders and all Sox10 Dom mice on (LabDiet 5LJ5) and water ad libitum. The B6C3Fe-a/a. Sox10 Dom line originated at Jackson Lab (Jackson Stock 000290) and was maintained by crosses to B6C3Fe-a/a WT female mice generated on site by crosses of C3FeLe.B6- a /J (Stock # 000198) females crossed to C57BL/6J (Stock # 000664) males. Pedigree offspring over 10 generations were euthanized at postnatal days 7-10 for GI tract collection that was dissected intact from stomach to anus. All offspring were genotyped for the Sox10 Dom allele using established methods ( Cantrell et al. 2004 ). Whole-mount acetylcholinesterase enzyme histochemistry of Sox10 Dom mutants was performed as described by Enomoto and colleagues with subsequent measurement of total gut length from pyloric sphincter to anus and aganglionosis based on acetylcholinesterase stain as previously reported ( Enomoto et al. 1998 ; Cantrell et al. 2004 ). The Tg(Phox2b-HIST2H2BE/Cerulean)1Sout mice (MGI: 5013571), Phox2b H2B-CFP, was maintained by crosses with female C3FeB6F1 ( Corpening et al. 2008 ). Timed matings between female C3FeB6F1 and male Phox2b H2B-CFP mice were conducted to produce fetal intestines at 16.5dpc used for snATAC-seq. Mouse genotyping Genotyping data was generated by the Center for Inherited Disease Research at Johns Hopkins University using the Illumina Mouse Linkage Panel (GoldenGate GS0006826-OPA), on a Sentrix array (1449 SNPs). Array data was analyzed with de novo clustering in GenomeStudio 2.0 with 1369 SNPs meeting quality metrics for analysis of which 876 were informative between C3HeBFeJLe-a/a and C57BL6/J. Analysis of Sox10 Dom extended pedigree using Genome-wide Efficient Mixed Model Analysis (GEMMA) GEMMA was used generally on the population and in sex-specific analyses to generate relatedness matrices (Zhou et al. 2012). We then performed SNP-based aganglionosis association analyses using GEMMA. Analysis code and further details can be found in Supplemental Materials and in supplementary files located on our Zenodo repository. Modifier intervals for downstream candidate gene analysis Sox10 Dom aganglionosis modifier intervals were defined as each interval of adjacent significant SNPs (p-value not adjusted for multiple tests) +/-0.5 megabase (Mb). The +/-0.5Mb window was chosen to account for associations that may not contain the causal SNP(s), gene(s), or locus (loci). We utilize shorthand for each GEMMA association run . Shorthand for each GEMMA analysis and modifier interval set used as input for our analysis pipeline ( Fig 7 ) is as follows: sex-regressed quantitative aganglionosis phenotype with zeroes included, BothSexAllChrom ; sex-regressed quantitative aganglionosis phenotype with zeroes included with removal of SNPs from chromosomes 5, 15, and X, BothSexNo15X5 ; sex-regressed with pseudo case-control “binary” phenotype, BinBothSexAllChrom ; male-and female-specific quantitative aganglionosis phenotype with zeroes included, Male and Female , respectively; male-and female-specific pseudo case-control “binary” phenotype, BinMale and BinFemale , respectively. Genes within modifier intervals were found using the UCSC Table Browser Tool, using the above modifier interval definitions as input ( Karolchik et al. 2004 ). 30 Reprocessing of whole gut 9.5-15.5 days post coitus scRNA-seq data from Zhao et al., 2022 ScRNA-seq data were downloaded from the Gene Expression Omnibus at accession GSE186525 ( Zhao et al. 2022 ). Seurat and SCTransform version 2 were utilized for quality control, clustering, and integration across age and sample ( Hao et al. 2024 , Seurat v5; Choudhary et al. 2022). Seurat’s FindAllMarkers function was used for differential gene expression analysis. Analysis code and further details for scRNA-seq processing can be found in Supplemental materials and on our Zenodo repository. Enteric neural crest-derived cell wavefront differential gene expression (DGE) overlapping with modifier intervals We downloaded differential gene expression data from Stavely et al. on the Gene Expression Omnibus at accession GSE217757 ( Stavely et al. 2023 ). This file was imported into R and genomic loci per significant differentially expressed gene (p-value less than 0.05) were determined using the same strategy as performed in the methods for imprinted genes (see Supplemental Materials). This was done for each GEMMA association-based modifier intervals, both non-sex specific and sex specific. SOX10 conserved binding motifs within modifier intervals SOX10 binding motifs from Gopinath et al. were downloaded and ported into UCSC Genome Browser’s LiftOver tool ( Hinrichs et al. 2006 ; Gopinath et al. 2016 ). Motif regions were processed sequentially from hg18 to mm10 in LiftOver. The binding motifs were then checked to see if they overlap with modifier intervals from all GEMMA associations using R code. Generation and processing of fetal ENS single nucleus assay for transposase-accessible chromatin-sequencing (snATAC-seq) data Cells from pooled 16.5dpc C3Fe.B6 Phox2b H2B-CFP+ fetal guts via cold-active protease dissociation as previously described ( Avila et al. 2025 ). CFP+ high (neuronal) and low (progenitors, glial) intensity cells were isolated by fluorescence-activated cell sorting (FACS) separately and simultaneously ( May-Zhang et al. 2021 ). Nonviable cells were removed via FAC sorting for 7-aminoactinomycin D+ stain. Nuclei were isolated post-FACS using the 10x Genomics protocol CG000209_Rev D. Nuclei counts were achieved via hemocytometer or Countess Automated Cell Counter (Thermofisher). Nuclei were encapsulated using 10x Genomics platform for snATAC-seq on high-and low-intensity CFP+ nuclei separately across two replicates from 6 and 12 pooled fetal guts. Libraries were sequenced on an Illumina NovaSeq 6000 using an S4 flow cell with custom read lengths to support ATAC-seq samples. Paired-end reads targeted greater than 70,000 reads per nucleus. Sequence FASTQs were processed with CellRanger ATAC pipeline version 1.2.0 by the Vanderbilt Technologies for advanced genomics (VANTAGE) shared resource ( Zheng et al. 2017 ). Sequences were aligned to the mm10. Quality control, batch correction, and subsequent analyses such as TF binding motif enrichment via Seurat, Signac, MACS2, and Harmony were performed as described in Supplementary Materials ( Zhang et al. 2008 ; Korsunsky et al. 2019 ; Stuart et al. 2022; Hao et al. 2024 ). Sox10 Dom aganglionosis modifier interval candidate gene prioritization pipeline from mouse datasets Omics analyses used in prioritization were merged into a final table that was filtered to only candidate genes that were in three or more lines of omics evidence ( Supp. Table 21, Table 7 ). A score of 0-3.5 was assigned to each candidate gene based on whether the modifier loci in which the candidate resides were 1) significant after multiple testing correction assigned 1 point, 2) LOD score of ≥3 assigned one point, and 3) overlapped with the modifier intervals from the F 1 -intercross study assigned 1.5 points ( Supp. Table 22 ; Owens et al. 2005 ). 7 We do this to consider whether modifiers in which candidate genes reside could be false-positive associations. Filtering Sox10 Dom aganglionosis modifier interval candidate genes for those with intron or exon variants with predicted high impact C3HeB/FeJ variants (VCF; mm39) for 30 candidate genes were used as input for the online Ensembl Variant Effect Predictor (VEP) for Mus musculus and were filtered for high predicted impact ( Hinrichs et al. 2006 ; Lawrence et al. 2009; McLaren et al. 2016 ). Overlap of Sox10 Dom aganglionosis modifier interval candidate genes with human HSCR and stool frequency GWAS summary statistics Summary statistics or analogous tables were downloaded for stool frequency and four different HSCR GWAS ( Garcia-Barcelo et al. 2009 ; Kim et al. 2014 ; Tang et al. 2016 ; Fadista et al. 2018 ; Bonfiglio et al. 2021 ). If all SNPs were available, false discovery rate (FDR) p-value adjustment was used, and SNPs were filtered based on FDR significance ( Garcia-Barcelo et al. 2009 ; Tang et al. 2016 ; Bonfiglio et al. 2021 ). Each summary statistics table was verified to use SNP coordinates on hg19. SNPs were tested for LD with candidate genes via precalculated LD blocks (Berisa et al. 2016). These genes were then overlapped with their mouse orthologs from Table 8 . Extraction of relevant phenotypes from Genebass-sourced Sox10 Dom aganglionosis modifier interval candidate gene-based PheWAS Exome-based PheWAS from Genebass for variants in each candidate gene were filtered to the top 20 associated phenotypes per candidate ( Karczewski et al. 2022 , Genebass 2024). These phenotypes were then filtered by common GI and nervous system terms as described in Supplementary Materials. AUTHOR CONTRIBUTIONS EMS 2 and JTB conceived the study. JTB performed data and statistical analyses and generated the figures. JAA and JRS performed data collection and analysis. JTB and EMS 2 drafted the manuscript. JRS and EMS 2 edited the manuscript. JTB generated shared data resource links. EMS 2 obtained funding. EMS 2 supervised the study. All authors revised and approved the manuscript. Supplementary Figure Legends Download figure Open in new tab Supplementary Fig 1. QQ Plots of Non-Sex-specific and Sex-specific GEMMA genome-wide scans on the Sox10 Dom mice population. QQ plots visualizing GEMMA genome-wide scan results with inclusion of all chromosomes (A) and exclusion of chromosomes 5, 15, and X (B) for the total quantitative aganglionosis percentage phenotype. C QQ plot visualizing results from a GEMMA genome-wide scan in which a binary phenotype—either the individual mouse has or does not have aganglionosis measured—was used. QQ plots are also shown visualizing GEMMA genome-wide scan results of female- (D) and male-specific (E) runs for the total quantitative percentage aganglionosis phenotype. QQ plots visualizing GEMMA genome-wide scan results of female- (F) and male-specific (G) runs for the binary phenotype—either the individual mouse has or does not have aganglionosis measured. Download figure Open in new tab Supplementary Fig 2. Stratification of allele-specific directions of effect by sex reveal sex effects. A-E,I-N The top 2 most significant SNPs per sex-specific genome-wide scans comparing the percent aganglionosis across individuals by genotype split by sex. F Chromosome 3’s top hit in the first GEMMA run shows a larger effect on males than females with rs13477043, and larger effects in females for genotype combinations for rs3667738, chromosome 8 (G). Similar differences in males and females are observed for rs13481145, the top hit for chromosome 11 (H). Each box plot shows distribution of percentage of aganglionosis and comparative statistics for each allele combination split by sex. Color of dots indicate the GEMMA runs with which each SNP is the peak associated SNP. See methods for the shorthand key for GEMMA association runs, which are the labels used here. Overall p: Kruskal-Wallis; internal p: Wilcoxon test. *, p<0.05; **, p<0.005; ***, p<0.0005; ****, p<0.00005; ns, not significant. A ‡ beside a SNP indicates significance of a SNP past multiple testing correction. Download figure Open in new tab Supplementary Fig 3. Sex-biased modifiers contain the highest number of differentially expressed genes in developing fetal gut cell types. Feature plots display expression via presence and intensity of purple split by developmental timepoint in the reprocessed Zhao et al. 2022 scRNA-seq dataset for candidate genes nearest to the top 2 associated SNPs or the most likely candidate gene based on known ENS developmental biology ( Ednrb ). Download figure Open in new tab Supplementary Fig 4. CellChat identifies modifier interval candidates via estimated cell signaling with neural crest cells. A Dot plot split by time point visualizing probability of activity of significant signaling pathways (y-axis) by cell type (x-axis). Shape of the dot indicates which modifier interval each signaling gene is within. Signaling genes within modifier intervals are underlined. Significance is represented by the size of the dot. Cell types on the left (orange bracket) represent signaling into neural crest cells, while types on the right (blue bracket) represent signaling out of neural crest cells to those cell types. B Dot plot showing expression of genes within aganglionosis modifier intervals and their ligand-receptor partners. Cell types have been consolidated to those that are in Fig 4 , and each cell type has split expression for developmental time in order, excluding 15.5dpc. Download figure Open in new tab Supplementary Fig 5. Expression of Sox10 Dom aganglionosis modifier interval candidate genes from wavefront ENCDC bulk RNA-seq and near conserved SOX10 binding motifs. A UMAP of the neural crest cells from Zhao et al., 2022 highlighting neuronal and progenitor cells. UMAPs of the neural crest cells from Zhao et al., 2022 showing expression in purple of candidate genes are shown that are upregulated in the migrating wavefront of ENCDCs ( B ) or are near or overlapping with conserved SOX10 binding motifs grouped by prior evidence ( C , other data modalities supporting gene as a candidate) and number of binding motifs ( D , two binding motifs; E, one binding motif). Download figure Open in new tab Supplementary Fig 6. Whole gut 16.5dpc Phox2b H2B-CFP+ snATAC-seq label transfer and derived putative cis-regulatory elements within modifier intervals. A Flow chart of analysis pipeline for differentially accessible chromatin contained within modifier intervals. B UMAP of scRNA-seq of 15.5dpc whole gut enteric nervous system cells annotated by supervised clustering used as a template for estimation of cell types in the snATAC-seq. C Estimation of cell types via Seurat and Signac’s LabelTransfer function used to annotate Fig 6B . Clusters from B are on the x-axis and clusters from Fig 6B are on the y-axis. D Annotations of differentially accessible (DA) peaks from all DA peaks (top) and those within modifier intervals split by cluster (bottom). ACKNOWLEDGEMENTS We are grateful to Ashley Cantrell, Elizabeth Stein, and Emily Ferguson for maintenance of the C3FeB6. Sox10 Dom pedigree and assisting with wholemount AChE gut staining; to Joan Breyer for helping grid DNA samples; and to Kevin M. Bradley for assistance with managing genotyping data. Genotyping was performed at the Center for Inherited Disease Research (Johns Hopkins University) via the NIH CIDR Program. We thank the Vanderbilt Technologies for Advanced Genomics (VANTAGE) shared resource for sequencing and CellRanger analysis for the snATAC-seq samples. We appreciate the suggestions made by Dr. Karl Broman on using a split aganglionosis phenotype approach to account for the zero-inflated data for the pedigree analysis. We thank Dr. Laura Reinholdt and colleagues for sharing the C3HeB/FeJ whole genome sequencing dataset that was vital for identifying risk alleles in our study. We are grateful to Drs. Rebecca Ihrie and Mary Chalkley for assisting with use of the Countess Automated cell counter. This work was supported in part by a Core Scholarship funded by P30DK058404 via the Vanderbilt University Medical Center’s Digestive Disease Research Center. HuC/D antibody was a kind gift from Vanda Lennon (Mayo Clinic). Funder Information Declared National Institutes of Health, https://ror.org/01cwqze88 , R01 DK60047 , R01 DK127178 , T32-HD007502 , 1F31DK137637-01 , R01-DK127178-03S1 Footnotes Grant Support: This work was funded by NIH award R01 DK60047 and R01 DK127178 to EMS 2 , with support for JTB on T32-HD007502 and NIDDK 1F31DK137637-01. JAA was supported by T32-HD007502 and R01-DK127178-03S1. Disclosures: Authors do not have any relevant conflicts of interest. Data Access: All Supplementary Tables, Code, fragment files, and snATAC-seq Seurat Object are shared at Zenodo: https://doi.org/10.5281/zenodo.17503749 https://doi.org/10.5281/zenodo.17503749 REFERENCES ↵ Avila JA , Benthal JT , Schafer JC , Flaherty DK , Southard-Smith EM . Single Cell Profiling in the Sox10 Dom Hirschsprung Mouse Implicates Hox genes in Enteric Neuron Trajectory Allocation . Cell Mol Gastroenterol Hepatol. Published online July 21 , 2025 . doi: 10.1016/j.jcmgh.2025.101590 OpenUrl CrossRef Berisa T , Pickrell JK . Approximately independent linkage disequilibrium blocks in human populations . Bioinformatics . 2016 ; 32 ( 2 ): 283 – 285 . doi: 10.1093/bioinformatics/btv546 OpenUrl CrossRef PubMed ↵ Bolk S , Pelet A , Hofstra RM , Angrist M , Salomon R , Croaker D , et al. A human model for multigenic inheritance: phenotypic expression in Hirschsprung disease requires both the RET gene and a new 9q31 locus . Proc Natl Acad Sci U S A . 2000 ; 97 ( 1 ): 268 – 273 . doi: 10.1073/pnas.97.1.268 OpenUrl Abstract / FREE Full Text ↵ Bonfiglio F , Liu X , Smillie C , Pandit A , Kurilshikov A , Bacigalupe R , et al. GWAS of stool frequency provides insights into gastrointestinal motility and irritable bowel syndrome . Cell Genom . 2021 ; 1 ( 3 ): 100069 . Published 2021 Dec 8. doi: 10.1016/j.xgen.2021.100069 OpenUrl CrossRef ↵ Borghini S , Di Duca M , Santamaria G , Vargiolu M , Bachetti T , Cargnin F , et al. Transcriptional regulation of TLX2 and impaired intestinal innervation: possible role of the PHOX2A and PHOX2B genes . Eur J Hum Genet . 2007 ; 15 ( 8 ): 848 – 855 . doi: 10.1038/sj.ejhg.5201852 OpenUrl CrossRef PubMed ↵ Broman KW , Sen S , Owens SE , Manichaikul A , Southard-Smith EM , Churchill GA . The X chromosome in quantitative trait locus mapping . Genetics . 2006 ; 174 ( 4 ): 2151 – 2158 . doi: 10.1534/genetics.106.061176 OpenUrl Abstract / FREE Full Text ↵ Cantrell VA , Owens SE , Chandler RL , Airey DC , Bradley KM , Smith JR , et al. Interactions between Sox10 and EdnrB modulate penetrance and severity of aganglionosis in the Sox10Dom mouse model of Hirschsprung disease [published correction appears in Hum Mol Genet. 2004 Dec 15;13(24):3241] . Hum Mol Genet . 2004 ; 13 ( 19 ): 2289 – 2301 . doi: 10.1093/hmg/ddh243 OpenUrl CrossRef PubMed Web of Science Choudhary S , Satija R . Comparison and evaluation of statistical error models for scRNA-seq . Genome Biol . 2022 ; 23 ( 1 ): 27 . Published 2022 Jan 18. doi: 10.1186/s13059-021-02584-9 OpenUrl CrossRef PubMed ↵ Corpening JC , Cantrell VA , Deal KK , Southard-Smith EM . A Histone2BCerulean BAC transgene identifies differential expression of Phox2b in migrating enteric neural crest derivatives and enteric glia . Dev Dyn . 2008 ; 237 ( 4 ): 1119 – 1132 . doi: 10.1002/dvdy.21498 OpenUrl CrossRef PubMed Web of Science ↵ Davis RJ , Shen W , Sandler YI , Amoui M , Purcell P , Maas R , et al. Dach1 mutant mice bear no gross abnormalities in eye, limb, and brain development and exhibit postnatal lethality . Mol Cell Biol . 2001 ; 21 ( 5 ): 1484 – 1490 . doi: 10.1128/MCB.21.5.1484-1490.2001 OpenUrl Abstract / FREE Full Text ↵ de Pontual L , Pelet A , Clement-Ziza M , Trochet D , Antonarakis SE , Attie-Bitach T , et al. Epistatic interactions with a common hypomorphic RET allele in syndromic Hirschsprung disease . Hum Mutat . 2007 ; 28 ( 8 ): 790 – 796 . doi: 10.1002/humu.20517 OpenUrl CrossRef PubMed Web of Science ↵ Dorman A , Binenbaum I , Abu-Toamih Atamni HJ , Chatziioannou A , Tomlinson I , Mott R , et al. Genetic mapping of novel modifiers for Apc Min induced intestinal polyps’ development using the genetic architecture power of the collaborative cross mice . BMC Genomics . 2021 ; 22 ( 1 ): 566 . Published 2021 Jul 22. doi: 10.1186/s12864-021-07890-x OpenUrl CrossRef PubMed ↵ Elworthy S , Pinto JP , Pettifer A , Cancela ML , Kelsh RN . Phox2b function in the enteric nervous system is conserved in zebrafish and is sox10-dependent . Mech Dev . 2005 ; 122 ( 5 ): 659 – 669 . doi: 10.1016/j.mod.2004.12.008 OpenUrl CrossRef PubMed Web of Science ↵ Enomoto H , Araki T , Jackman A , Heuckeroth RO , Snider WD , Johnson Jr. EM , et al. GFR alpha1-deficient mice have deficits in the enteric nervous system and kidneys . Neuron . 1998 ; 21 ( 2 ): 317 – 324 . doi: 10.1016/s0896-6273(00)80541-3 OpenUrl CrossRef PubMed Web of Science ↵ Fadista J , Lund M , Skotte L , Geller F , Nandakumar P , Chatterjee S , et al. Genome-wide association study of Hirschsprung disease detects a novel low-frequency variant at the RET locus . Eur J Hum Genet . 2018 ; 26 ( 4 ): 561 – 569 . doi: 10.1038/s41431-017-0053-7 OpenUrl CrossRef PubMed ↵ Finzsch M , Stolt CC , Lommes P , Wegner M . Sox9 and Sox10 influence survival and migration of oligodendrocyte precursors in the spinal cord by regulating PDGF receptor alpha expression . Development . 2008 ; 135 ( 4 ): 637 – 646 . doi: 10.1242/dev.010454 OpenUrl Abstract / FREE Full Text ↵ Fu M , Lui VC , Sham MH , Cheung AN , Tam PK . HOXB5 expression is spatially and temporarily regulated in human embryonic gut during neural crest cell colonization and differentiation of enteric neuroblasts . Dev Dyn . 2003 ; 228 ( 1 ): 1 – 10 . doi: 10.1002/dvdy.10350 OpenUrl CrossRef PubMed Web of Science ↵ Gao T , Wright-Jin EC , Sengupta R , Anderson JB , Heuckeroth RO . Cell-autonomous retinoic acid receptor signaling has stage-specific effects on mouse enteric nervous system . JCI Insight . 2021 ; 6 ( 10 ): e145854 . Published 2021 May 24. doi: 10.1172/jci.insight.145854 OpenUrl CrossRef ↵ Garcia-Barcelo MM , Tang CS , Ngan ES , Lui CV , Chen Y , So M , et al. Genome-wide association study identifies NRG1 as a susceptibility locus for Hirschsprung’s disease . Proc Natl Acad Sci U S A . 2009 ; 106 ( 8 ): 2694 – 2699 . doi: 10.1073/pnas.0809630105 OpenUrl Abstract / FREE Full Text Genebass: Gene-based association summary statistics. https://app.genebass.org/ . Accessed 14 November 2024 . ↵ Gopinath C , Law WD , Rodríguez-Molina JF , Prasad AB , Song L , Crawford GE , et al. Stringent comparative sequence analysis reveals SOX10 as a putative inhibitor of glial cell differentiation . BMC Genomics . 2016 ; 17 ( 1 ): 887 . Published 2016 Nov 7. doi: 10.1186/s12864-016-3167-3 OpenUrl CrossRef PubMed ↵ Hakami RM , Hou L , Baxter LL , Loftus SK , Southard-Smith EM , Incao A , et al. Genetic evidence does not support direct regulation of EDNRB by SOX10 in migratory neural crest and the melanocyte lineage . Mech Dev . 2006 ; 123 ( 2 ): 124 – 134 . doi: 10.1016/j.mod.2005.11.004 OpenUrl CrossRef PubMed Web of Science ↵ Hao Y , Stuart T , Kowalski MH , Choudhary S , Hoffman P , Hartman A , et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis . Nat Biotechnol . 2024 ; 42 ( 2 ): 293 – 304 . doi: 10.1038/s41587-023-01767-y OpenUrl CrossRef PubMed ↵ Hinrichs AS , Karolchik D , Baertsch R , Barber GP , Bejerano G , Clawson H , et al. The UCSC Genome Browser Database: update 2006 . Nucleic Acids Res . 2006 ; 34 (Database issue):D590-D598. doi: 10.1093/nar/gkj144 OpenUrl CrossRef PubMed Web of Science Integrative analysis in Seurat v5. https://satijalab.org/seurat/articles/seurat5_integration . Accessed 14 November 2024 . ↵ Jiang M , Li C , Cao G , Yang D , Zhang X , Yang L , et al. Effects of NRG1 Polymorphisms on Hirschsprung’s Disease Susceptibility: A Meta-analysis . Sci Rep . 2017 ; 7 ( 1 ): 9913 . Published 2017 Aug 30. doi: 10.1038/s41598-017-10477-w OpenUrl CrossRef PubMed ↵ Jiang Q , Ho YY , Hao L , Nichols Berrios C , Chakravarti A . Copy number variants in candidate genes are genetic modifiers of Hirschsprung disease . PLoS One . 2011 ; 6 ( 6 ): e21219 . doi: 10.1371/journal.pone.0021219 OpenUrl CrossRef PubMed ↵ Jin S , Guerrero-Juarez CF , Zhang L , Chang I , Ramos R , Kuan CH , et al. Inference and analysis of cell-cell communication using CellChat . Nat Commun . 2021 ; 12 ( 1 ): 1088 . Published 2021 Feb 17. doi: 10.1038/s41467-021-21246-9 OpenUrl CrossRef PubMed ↵ Jin S , Plikus MV , Nie Q . CellChat for systematic analysis of cell-cell communication from single-cell transcriptomics . Nat Protoc. Published online September 16 , 2024 . doi: 10.1038/s41596-024-01045-4 OpenUrl CrossRef PubMed ↵ Kaern M , Elston TC , Blake WJ , Collins JJ . Stochasticity in gene expression: from theories to phenotypes . Nat Rev Genet . 2005 ; 6 ( 6 ): 451 – 464 . doi: 10.1038/nrg1615 OpenUrl CrossRef PubMed Web of Science ↵ Kalra S , Lanno S , Sanchez G , Coolon JD. cis-and trans-regulatory contributions to a hierarchy of factors influencing gene expression variation . Commun Biol . 2024 ; 7 ( 1 ): 1563 . Published 2024 Nov 25. doi: 10.1038/s42003-024-07255-6 OpenUrl CrossRef PubMed ↵ Karczewski KJ , Solomonson M , Chao KR , Goodrich JK , Tiao G , Lu W , et al. Systematic single-variant and gene-based association testing of thousands of phenotypes in 394,841 UK Biobank exomes . Cell Genom . 2022 ; 2 ( 9 ): 100168 . Published 2022 Aug 15. doi: 10.1016/j.xgen.2022.100168 OpenUrl CrossRef PubMed ↵ Karolchik D , Hinrichs AS , Furey TS , Roskin KM , Sugnet CW , Haussler D , et al. The UCSC Table Browser data retrieval tool . Nucleic Acids Res . 2004 ; 32 (Database issue):D493-D496. doi: 10.1093/nar/gkh103 OpenUrl CrossRef PubMed Web of Science ↵ Kim JH , Cheong HS , Sul JH , Seo JM , Kim DY , Oh JT , et al. A genome-wide association study identifies potential susceptibility loci for Hirschsprung disease . PLoS One . 2014 ; 9 ( 10 ): e110292 . Published 2014 Oct 13. doi: 10.1371/journal.pone.0110292 OpenUrl CrossRef PubMed ↵ Kim YK , Lee H , Ismail T , Kim Y , Lee HS . Dach1 regulates neural crest migration during embryonic development . Biochem Biophys Res Commun . 2020 ; 527 ( 4 ): 896 – 901 . doi: 10.1016/j.bbrc.2020.05.009 OpenUrl CrossRef PubMed ↵ Korsunsky I , Millard N , Fan J , Slowikowski K , Zhang F , Wei K , et al. Fast, sensitive and accurate integration of single-cell data with Harmony . Nat Methods . 2019 ; 16 ( 12 ): 1289 – 1296 . doi: 10.1038/s41592-019-0619-0 OpenUrl CrossRef PubMed ↵ Kuil LE , MacKenzie KC , Tang CS , Windster JD , Le TL , Karim A , et al. Size matters: Large copy number losses in Hirschsprung disease patients reveal genes involved in enteric nervous system development . PLoS Genet . 2021 ; 17 ( 8 ): e1009698 . Published 2021 Aug 6. doi: 10.1371/journal.pgen.1009698 OpenUrl CrossRef PubMed ↵ Lake JI , Heuckeroth RO . Enteric nervous system development: migration, differentiation, and disease . Am J Physiol Gastrointest Liver Physiol . 2013 ; 305 ( 1 ): G1 – G24 . doi: 10.1152/ajpgi.00452.2012 OpenUrl CrossRef PubMed Web of Science ↵ Lane PW , Liu HM . Association of megacolon with a new dominant spotting gene (Dom) in the mouse . J Hered . 1984 ; 75 ( 6 ): 435 – 439 . doi: 10.1093/oxfordjournals.jhered.a109980 OpenUrl CrossRef PubMed Web of Science Lawrence M , Gentleman R , Carey V. rtracklayer: an R package for interfacing with genome browsers . Bioinformatics . 2009 ; 25 ( 14 ): 1841 – 1842 . doi: 10.1093/bioinformatics/btp328 OpenUrl CrossRef PubMed Web of Science ↵ Li N , Kelsh RN , Croucher P , Roehl HH . Regulation of neural crest cell fate by the retinoic acid and Pparg signalling pathways . Development . 2010 ; 137 ( 3 ): 389 – 394 . doi: 10.1242/dev.044164 OpenUrl Abstract / FREE Full Text ↵ Li C , Hu R , Hou N , Wang Y , Wang Z , Yang T , et al. Alteration of the Retinoid Acid-CBP Signaling Pathway in Neural Crest Induction Contributes to Enteric Nervous System Disorder . Front Pediatr . 2018 ; 6 : 382 . Published 2018 Dec 3. doi: 10.3389/fped.2018.00382 OpenUrl CrossRef PubMed ↵ Li Z , Lui KN , Lau ST , Lai FPL , Li P , Chung PHY et al. Transcriptomics of Hirschsprung disease patient-derived enteric neural crest cells reveals a role for oxidative phosphorylation . Nat Commun . 2023 ; 14 ( 1 ): 2157 . Published 2023 Apr 15. doi: 10.1038/s41467-023-37928-5 OpenUrl CrossRef PubMed ↵ Lippert C , Listgarten J , Liu Y , Kadie CM , Davidson RI , Heckerman D . FaST linear mixed models for genome-wide association studies . Nat Methods . 2011 ; 8 ( 10 ): 833 – 835 . Published 2011 Sep 4. doi: 10.1038/nmeth.1681 OpenUrl CrossRef PubMed Web of Science ↵ Liu SJ , Casey-Clyde T , Cho NW , Swinderman J , Pekmezci M , Dougherty MC , et al. Epigenetic reprogramming shapes the cellular landscape of schwannoma . Nat Commun . 2024 ; 15 ( 1 ): 476 . Published 2024 Jan 12. doi: 10.1038/s41467-023-40408-5 OpenUrl CrossRef PubMed ↵ López SH , Avetisyan M , Wright CM , Mesbah K , Kelly RG , Moon AM , et al. Loss of Tbx3 in murine neural crest reduces enteric glia and causes cleft palate, but does not influence heart development or bowel transit . Dev Biol . 2018 ; 444 Suppl 1(Suppl 1):S337-S351. doi: 10.1016/j.ydbio.2018.09.017 OpenUrl CrossRef ↵ Lorè NI , Sipione B , He G , Strug LJ , Atamni HJ , Dorman A , et al. Collaborative Cross Mice Yield Genetic Modifiers for Pseudomonas aeruginosa Infection in Human Lung Disease . mBio . 2020 ; 11 ( 2 ): e00097 – 20 . Published 2020 Mar 3. doi: 10.1128/mBio.00097-20 OpenUrl CrossRef PubMed ↵ Lui VC , Cheng WW , Leon TY , Lau DKC , Garcia-Barcelo MM , Miao XP , et al. Perturbation of hoxb5 signaling in vagal neural crests down-regulates ret leading to intestinal hypoganglionosis in mice [published correction appears in Gastroenterology. 2008 Aug;135(2):713. Garcia-Bareclo, Maria-Mercedes [corrected to Garcia-Barcelo, Maria-Mercedes]] . Gastroenterology . 2008 ; 134 ( 4 ): 1104 – 1115 . doi: 10.1053/j.gastro.2008.01.028 OpenUrl CrossRef PubMed Web of Science ↵ Majd H , Samuel RM , Cesiulis A , Ramirez JT , Kalantari A , Barber K , et al. Engrafted nitrergic neurons derived from hPSCs improve gut dysmotility in mice . Nature. Published online June 25 , 2025 . doi: 10.1038/s41586-025-09208-3 OpenUrl CrossRef ↵ Maka M , Stolt CC , Wegner M . Identification of Sox8 as a modifier gene in a mouse model of Hirschsprung disease reveals underlying molecular defect . Dev Biol . 2005 ; 277 ( 1 ): 155 – 169 . doi: 10.1016/j.ydbio.2004.09.014 OpenUrl CrossRef PubMed Web of Science ↵ Mallory SB , Wiener E , Nordlund JJ . Waardenburg’s syndrome with Hirschsprung’s disease: a neural crest defect . Pediatr Dermatol . 1986 ; 3 ( 2 ): 119 – 124 . doi: 10.1111/j.1525-1470.1986.tb00501.x OpenUrl CrossRef PubMed ↵ May-Zhang AA , Tycksen E , Southard-Smith AN , Deal KK , Benthal JT , Buehler DP , et al. Combinatorial Transcriptional Profiling of Mouse and Human Enteric Neurons Identifies Shared and Disparate Subtypes In Situ . Gastroenterology . 2021 ; 160 ( 3 ): 755 – 770 .e26. doi: 10.1053/j.gastro.2020.09.032 OpenUrl CrossRef PubMed McCallion AS , Stames E , Conlon RA , Chakravarti A . Phenotype variation in two-locus mouse models of Hirschsprung disease: tissue-specific interaction between Ret and Ednrb . Proc Natl Acad Sci U S A . 2003 ; 100 ( 4 ): 1826 – 1831 . doi: 10.1073/pnas.0337540100 OpenUrl Abstract / FREE Full Text ↵ McLaren W , Gil L , Hunt SE , Riat HS , Ritchie GRS , Thormann A , et al. The Ensembl Variant Effect Predictor . Genome Biol . 2016 ; 17 ( 1 ): 122 . Published 2016 Jun 6. doi: 10.1186/s13059-016-0974-4 OpenUrl CrossRef PubMed ↵ Memic F , Knoflach V , Morarach K , Sadler R , Laranjeira C , Hjerling-Leffler J , et al. Transcription and Signaling Regulators in Developing Neuronal Subtypes of Mouse and Human Enteric Nervous System . Gastroenterology . 2018 ; 154 ( 3 ): 624 – 636 . doi: 10.1053/j.gastro.2017.10.005 OpenUrl CrossRef PubMed ↵ Mo J , Long R , Fantauzzo KA . Pdgfra and Pdgfrb Genetically Interact in the Murine Neural Crest Cell Lineage to Regulate Migration and Proliferation . Front Physiol . 2020 ; 11 : 588901 . Published 2020 Nov 2. doi: 10.3389/fphys.2020.588901 OpenUrl CrossRef PubMed ↵ Nagashimada M , Ohta H , Li C , Nakao K , Uesaka T , Brunet JF , et al. Autonomic neurocristopathy-associated mutations in PHOX2B dysregulate Sox10 expression . J Clin Invest . 2012 ; 122 ( 9 ): 3145 – 3158 . doi: 10.1172/JCI63401 OpenUrl CrossRef PubMed Web of Science ↵ Nagy N , Goldstein AM . Enteric nervous system development: A crest cell’s journey from neural tube to colon . Semin Cell Dev Biol . 2017 ; 66 : 94 – 106 . doi: 10.1016/j.semcdb.2017.01.006 OpenUrl CrossRef PubMed ↵ Owens SE , Broman KW , Wiltshire T , Elmore JB , Bradley KM , Smith JR , et al. Genome-wide linkage identifies novel modifier loci of aganglionosis in the Sox10Dom model of Hirschsprung disease . Hum Mol Genet . 2005 ; 14 ( 11 ): 1549 – 1558 . doi: 10.1093/hmg/ddi163 OpenUrl CrossRef PubMed Web of Science ↵ Parisi MA , Kapur RP , Neilson I , Hofstra RMW , Holloway LW , Michaelis RC , et al. Hydrocephalus and intestinal aganglionosis: is L1CAM a modifier gene in Hirschsprung disease? Am J Med Genet . 2002 ; 108 ( 1 ): 51 – 56 . doi: 10.1002/ajmg.10185 OpenUrl CrossRef PubMed Web of Science ↵ Pattyn A , Morin X , Cremer H , Goridis C , Brunet JF . The homeobox gene Phox2b is essential for the development of autonomic neural crest derivatives . Nature . 1999 ; 399 ( 6734 ):366-370. doi: 10.1038/20700 OpenUrl CrossRef PubMed Web of Science ↵ Perez G , Barber GP , Benet-Pages A , Casper J , Clawson H , Diekhans M , et al. The UCSC Genome Browser database: 2025 update . Nucleic Acids Res. Published online October 26 , 2024 . doi: 10.1093/nar/gkae974 OpenUrl CrossRef PubMed ↵ Pingault V , Bondurand N , Kuhlbrodt K , Goerich DE , Préhu MO , Puliti A , et al. SOX10 mutations in patients with Waardenburg-Hirschsprung disease . Nat Genet . 1998 ; 18 ( 2 ): 171 – 173 . doi: 10.1038/ng0298-171 OpenUrl CrossRef PubMed Web of Science ↵ Sato Y , Heuckeroth RO . Retinoic acid regulates murine enteric nervous system precursor proliferation, enhances neuronal precursor differentiation, and reduces neurite growth in vitro . Dev Biol . 2008 ; 320 ( 1 ): 185 – 198 . doi: 10.1016/j.ydbio.2008.05.524 OpenUrl CrossRef PubMed Web of Science ↵ Shih CH , Fay J . Cis -regulatory variants affect gene expression dynamics in yeast . Elife . 2021 ; 10 : e68469 . Published 2021 Aug 9. doi: 10.7554/eLife.68469 OpenUrl CrossRef PubMed ↵ Shull LC , Lencer ES , Kim HM , Goyama S , Kurokawa M , Costello JC , et al. PRDM paralogs antagonistically balance Wnt/β-catenin activity during craniofacial chondrocyte differentiation . Development . 2022 ; 149 ( 4 ):dev200082. doi: 10.1242/dev.200082 OpenUrl CrossRef ↵ Simpson MJ , Zhang DC , Mariani M , Landman KA , Newgreen DF . Cell proliferation drives neural crest cell invasion of the intestine . Dev Biol . 2007 ; 302 ( 2 ): 553 – 568 . doi: 10.1016/j.ydbio.2006.10.017 OpenUrl CrossRef PubMed Web of Science ↵ Southard-Smith EM , Angrist M , Ellison JS , Agarwala R , Baxevanis AD , Chakravarti A , et al. The Sox10(Dom) mouse: modeling the genetic variation of Waardenburg-Shah (WS4) syndrome . Genome Res . 1999 ; 9 ( 3 ): 215 – 225 . OpenUrl Abstract / FREE Full Text ↵ Stanchina L , Baral V , Robert F , Pingault V , Lemort N , Pachnis V , et al. Interactions between Sox10, Edn3 and Ednrb during enteric nervous system and melanocyte development . Dev Biol . 2006 ; 295 ( 1 ): 232 – 249 . doi: 10.1016/j.ydbio.2006.03.031 OpenUrl CrossRef PubMed Web of Science ↵ Stavely R , Hotta R , Guyer RA , Picard N , Rahman AA , Omer M , et al. A distinct transcriptome characterizes neural crest-derived cells at the migratory wavefront during enteric nervous system development . Development . 2023 ; 150 ( 5 ):dev201090. doi: 10.1242/dev.201090 OpenUrl CrossRef ↵ Stuart T , Srivastava A , Madad S , Lareau CA , Satija R. Single-cell chromatin state analysis with Signac [published correction appears in Nat Methods . 2022 Feb; 19 ( 2 ): 257 . doi: 10.1038/s41592-022-01393-7]. Nat Methods . 2021 ;18( 11 ):1333-1341. doi: 10.1038/s41592-021-01282-5 OpenUrl CrossRef ↵ Tang CS , Gui H , Kapoor A , Kim JH , Luzón-Toro B , Pelet A , et al. Trans-ethnic meta-analysis of genome-wide association studies for Hirschsprung disease . Hum Mol Genet . 2016 ; 25 ( 23 ): 5265 – 5275 . doi: 10.1093/hmg/ddw333 OpenUrl CrossRef PubMed ↵ Taylor MK , Yeager K , Morrison SJ . Physiological Notch signaling promotes gliogenesis in the developing peripheral and central nervous systems . Development . 2007 ; 134 ( 13 ): 2435 – 2447 . doi: 10.1242/dev.005520 OpenUrl Abstract / FREE Full Text ↵ Uribe RA , Hong SS , Bronner ME . Retinoic acid temporally orchestrates colonization of the gut by vagal neural crest cells . Dev Biol . 2018 ; 433 ( 1 ): 17 – 32 . doi: 10.1016/j.ydbio.2017.10.021 OpenUrl CrossRef PubMed Wakamatsu Y , Uchikawa M . The many faces of Sox2 function in neural crest development . Dev Growth Differ . 2021 ; 63 ( 1 ): 93 – 99 . doi: 10.1111/dgd.12705 OpenUrl CrossRef PubMed ↵ Wang F , Hu W , Xian J , Ohnuma S , Brenton JD . The Xenopus Tgfbi is required for embryogenesis through regulation of canonical Wnt signalling . Dev Biol . 2013 ; 379 ( 1 ): 16 – 27 . doi: 10.1016/j.ydbio.2012.11.010 OpenUrl CrossRef PubMed ↵ Wang B , Fang X , Sun X , Du C , Zhou L , Lv X , et al. m 6 A demethylase ALKBH5 suppresses proliferation and migration of enteric neural crest cells by regulating TAGLN in Hirschsprung’s disease . Life Sci . 2021 ; 278 : 119577 . doi: 10.1016/j.lfs.2021.119577 OpenUrl CrossRef ↵ Wright CM , Schneider S , Smith-Edwards KM , Mafra F , Leembruggen AJ , Gonzalez MV , et al. scRNA-Seq Reveals New Enteric Nervous System Roles for GDNF, NRTN, and TBX3 . Cell Mol Gastroenterol Hepatol . 2021 ; 11 ( 5 ): 1548 – 1592 .e1. doi: 10.1016/j.jcmgh.2020.12.014 OpenUrl CrossRef ↵ Wu K , Katiyar S , Li A , Liu M , Ju X , Popov VM , et al. Dachshund inhibits oncogene-induced breast cancer cellular migration and invasion through suppression of interleukin-8 . Proc Natl Acad Sci U S A . 2008 ; 105 ( 19 ): 6924 – 6929 . doi: 10.1073/pnas.0802085105 OpenUrl Abstract / FREE Full Text ↵ Yang J , Zaitlen NA , Goddard ME , Visscher PM , Price AL . Advantages and pitfalls in the application of mixed-model association methods . Nat Genet . 2014 ; 46 ( 2 ): 100 – 106 . doi: 10.1038/ng.2876 OpenUrl CrossRef PubMed ↵ Young HM , Jones BR , McKeown SJ . The projections of early enteric neurons are influenced by the direction of neural crest cell migration . J Neurosci . 2002 ; 22 ( 14 ): 6005 – 6018 . doi: 10.1523/JNEUROSCI.22-14-06005.2002 OpenUrl Abstract / FREE Full Text ↵ Young HM , Bergner AJ , Anderson RB , Enomoto H , Milbrandt J , Newgreen DF , et al. Dynamics of neural crest-derived cell migration in the embryonic mouse gut . Dev Biol . 2004 ; 270 ( 2 ): 455 – 473 . doi: 10.1016/j.ydbio.2004.03.015 OpenUrl CrossRef PubMed Web of Science ↵ Young HM , Bergner AJ , Simpson MJ , McKeown SJ , Hao MM , Anderson CR , et al. Colonizing while migrating: how do individual enteric neural crest cells behave? . BMC Biol . 2014 ; 12 : 23 . Published 2014 Mar 26. doi: 10.1186/1741-7007-12-23 OpenUrl CrossRef PubMed ↵ Zhang Y , Liu T , Meyer CA , Eeckhoute J , Johnson DS , Bernstein BE , et al. Model-based analysis of ChIP-Seq (MACS) . Genome Biol . 2008 ; 9 ( 9 ): R137 . doi: 10.1186/gb-2008-9-9-r137 OpenUrl CrossRef PubMed ↵ Zhang D , Rollo BN , Nagy N , Stamp L , Newgreen DF . The enteric neural crest progressively loses capacity to form enteric nervous system . Dev Biol . 2019 ; 446 ( 1 ): 34 – 42 . doi: 10.1016/j.ydbio.2018.11.017 OpenUrl CrossRef PubMed ↵ Zhao L , Song W , Chen YG . Mesenchymal-epithelial interaction regulates gastrointestinal tract development in mouse embryos . Cell Rep . 2022 ; 40 ( 2 ): 111053 . doi: 10.1016/j.celrep.2022.111053 OpenUrl CrossRef PubMed ↵ Zheng GX , Terry JM , Belgrader P , Ryvkin P , Bent ZW , Wilson R , et al. Massively parallel digital transcriptional profiling of single cells . Nat Commun . 2017 ; 8 : 14049 . Published 2017 Jan 16. doi: 10.1038/ncomms14049 OpenUrl CrossRef PubMed Zhou X , Stephens M . Genome-wide efficient mixed-model analysis for association studies . Nat Genet . 2012 ; 44 ( 7 ): 821 – 824 . Published 2012 Jun 17. doi: 10.1038/ng.2310 OpenUrl CrossRef PubMed ↵ Zhu JJ , Kam MK , Garcia-Barceló MM , Tam PK , Lui VC . HOXB5 binds to multi-species conserved sequence (MCS+9.7) of RET gene and regulates RET expression . Int J Biochem Cell Biol . 2014 ; 51 : 142 – 149 . doi: 10.1016/j.biocel.2014.04.013 OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted November 04, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Combinatorial multiomic analysis from a pedigree of Sox10Dom Hirschsprung mice implicates Dach1 as a modifier of Enteric Nervous System development Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Combinatorial multiomic analysis from a pedigree of Sox10 Dom Hirschsprung mice implicates Dach1 as a modifier of Enteric Nervous System development Joseph T. Benthal , Justin A. Avila , Jeffrey R. Smith , E. Michelle Southard-Smith bioRxiv 2025.11.02.686123; doi: https://doi.org/10.1101/2025.11.02.686123 Share This Article: Copy Citation Tools Combinatorial multiomic analysis from a pedigree of Sox10 Dom Hirschsprung mice implicates Dach1 as a modifier of Enteric Nervous System development Joseph T. Benthal , Justin A. Avila , Jeffrey R. Smith , E. Michelle Southard-Smith bioRxiv 2025.11.02.686123; doi: https://doi.org/10.1101/2025.11.02.686123 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Genomics Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41910) Biophysics (21436) Cancer Biology (18576) Cell Biology (25480) Clinical Trials (138) Developmental Biology (13368) Ecology (19887) Epidemiology (2067) Evolutionary Biology (24302) Genetics (15598) Genomics (22482) Immunology (17726) Microbiology (40360) Molecular Biology (17163) Neuroscience (88534) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4821) Physiology (7637) Plant Biology (15129) Scientific Communication and Education (2045) Synthetic Biology (4290) Systems Biology (9817) Zoology (2269)