Skip to main page content
U.S. flag

An official website of the United States government

Dot gov

The .gov means it’s official.
Federal government websites often end in .gov or .mil. Before sharing sensitive information, make sure you’re on a federal government site.

Https

The site is secure.
The https:// ensures that you are connecting to the official website and that any information you provide is encrypted and transmitted securely.

Access keys NCBI Homepage MyNCBI Homepage Main Content Main Navigation
. 2025;19(12):101590.
doi: 10.1016/j.jcmgh.2025.101590. Epub 2025 Jul 21.

Single Cell Profiling in the Sox10Dom Hirschsprung Mouse Implicates Hox Genes in Enteric Neuron Trajectory Allocation

Justin A Avila  1 , Joseph T Benthal  2 , Jenny C Schafer  3 , David K Flaherty  4 , E Michelle Southard-Smith  5
Affiliations

Single Cell Profiling in the Sox10Dom Hirschsprung Mouse Implicates Hox Genes in Enteric Neuron Trajectory Allocation

Justin A Avila et al. Cell Mol Gastroenterol Hepatol. 2025.

Abstract

Background & aims: Enteric nervous system (ENS) development requires migration, proliferation, and differentiation of progenitors for normal gastrointestinal (GI) motility. Sox10 deficit causes aganglionosis, modeling Hirschsprung disease (HSCR), and disrupts ratios of postnatal enteric neurons in proximal ganglionated bowel. How Sox10 deficiency alters enteric neuron ratios is unclear. Sox10's prominent expression in enteric neural crest-derived progenitors (ENCPs) and lack of this gene in mature enteric neurons led us to examine Sox10Dom effects in early ENS development.

Methods: Immunohistochemistry localized SOX10 in the developing ENS relative to HuC/D. ENS progenitors, developing neurons, and enteric glia were isolated from Sox10+/+ and Sox10Dom/+ littermates for single-cell RNA sequencing (scRNA-seq). scRNA-seq data was processed to identify cell type-specific markers, differentially expressed genes, cell fate trajectories, and gene regulatory network activity between genotypes. Hybridization chain reaction (HCR) coupled with immunohistochemistry validated expression changes.

Results: SOX10 protein was detected in early ENS neurons. scRNA-seq profiles detected three neuronal trajectories emerging via two transition pathways accompanied by elevated activity of Hox gene regulatory networks (GRN). Sox10Dom/+ scRNA-seq profiles exhibited a novel progenitor cluster, reduced numbers of cells in transitional states, and shifts in cell abundance between neuronal trajectories. Hoxa6 was differentially expressed in the neuronal trajectories impacted in Sox10Dom/+ mutants, and HCR identified altered Hoxa6 expression in early developing neurons of Sox10Dom/+ ENS.

Conclusions: Sox10Dom/+ mutation shifts enteric neuron types by altering neuronal trajectories early in ENS development. Multiple neurogenic transcription factors are reduced in Sox10Dom/+ scRNA-seq profiles. This work is the first to correlate changes in Hox expression, notably Hoxa6, with alterations in enteric neuron trajectories.

Keywords: Enteric Nervous System; Hirschsprung Disease; Hoxa6; Neural Crest; Sox10.

PubMed Disclaimer

Figures

None
Graphical abstract
Figure 1
Figure 1
Sox10-Phox2b transgene expression in emerging enteric neurons and glia facilitates capture of differentiating populations during ENS development. (A) Confocal images of 13.5 dpc fetal hindgut ∼50 microns proximal to the migrating ENP wavefront labeled by Sox10-YFP and Phox2b-CFP transgene expression detect SOX10 protein by IHC in early differentiating neurons that express HuC/D. Sox10-YFP transgene expression is already reduced in these HuC/D+ cells while SOX10 protein persists. (B) Transgenic lines crossed to enable FACS isolation of ENS populations marked by expression of Sox10-YFP in progenitors/enteric glia and Phox2b-CFP in progenitors/developing neurons. (C) Wholemount fluorescent images of YFP expression in Sox10+/+ and Sox10Dom/+ fetal intestines at 15.5 dpc. Dotted lines emphasize the hindgut borders in Sox10Dom/+ mutant where YFP cells are lacking in the colon. (C) FACS profile of ENPs sorted from 15.5 dpc fetal intestine shows gates that captured all transgene expressing cells (pink border shape) where dual transgene positivity is pseudo-colored green, Sox10-YFP+ in yellow, and Phox2b-CFP+ in blue.
Figure 2
Figure 2
Profiled ENS populations from dual transgenic strategy capture ENCPs and multiple developing neuronal trajectories. (A) UMAP plots display similar scRNA-seq clusters for both genotypes prior to data merge (Supplementary Tables 1 and 2). (B) UMAP plots for each genotype after data merge shows overlap and similar cluster distributions (Supplementary Tables 3 and 4). (C) UMAP plot showing the main ENS populations comprised of 19 total clusters (Supplementary Tables 5 and 6). Other clusters consisting of Schwann-like cells (Cluster 13 panel B) and low-quality cells with high 18s rRNA gene expression, high mitochondrial and low RNA feature counts (clusters 1, 14, 10 panel B) were removed in this view. (D) UMAP representations of marker genes for progenitors (Sox10, Sox2) and developing neurons (Elavl4, Snap25) reveal the cluster distributions of these populations. A maximum normalized expression threshold was set to 2 in these images to highlight most expression of these genes. (E) UMAP of clusters shared by Sox10+/+ and Sox10Dom/+ classified by marker gene expression and GO analysis (Supplementary Tables 7 and 8) identified two main groups: progenitors (green) and neurons (blue). (F) Dot plot of known marker genes for ENS progenitors, neuroblasts, neurons, and glia per cluster.
Figure 3
Figure 3
Differential gene expression distinguishes 3 neuronal trajectories in 15.5 dpc scRNA-seq data. (A) Dot plot displays gene expression detected in neuronal trajectories for the aggregated Sox10+/+ and Sox10Dom/+ data. (B) UMAPs highlight the specific trajectories designated here as Branch A (NC6, NC4, and NC9), Branch B (NC2 and NC0), and Branch C (NC3, NC5, and NC7).
Figure 4
Figure 4
Differential gene expression highlights transitional paths and gene signatures marking early neuronal clusters in 15.5dpc scRNA-seq data. (A) Dot plot displays gene expression detected in early neuronal branch clusters and in the transitional neuroblast populations for the aggregated Sox10+/+ and Sox10Dom/+ data. (B) UMAPs highlight clusters early in the trajectory of branch B (NC2) and branch C (NC3, NC5). (C) UMAP feature plots display expression of differentially expressed genes in early Branch B (Dlx6os1, Zfp503, Ptprk) and Branch C (Hoxa6, Adcyap1, Hoxb8) trajectories. (D) UMAPs highlight transitional neuroblast populations labeled as Path1 (NC1, NC8) and Path2 (NC11, NC10).
Figure 5
Figure 5
Deep scRNA-seq at 15.5 dpc identifies greater complexity of early ENS neurons. (A) UMAP visualization of reprocessed Morarach 15.5-dpc and 18.5-dpc data (leftmost panels). Note that sex-specific genes were removed as was done by Morarach et al, 2021. Middle panels: UMAP of reprocessed Morarach data with cluster labels from the current Sox10 +/+ 15.5-dpc scRNA-seq data overlaid after Label Transfer. Right panels: UMAP of reprocessed Morarach data with cluster labels overlaid after Label Transfer from the present studies Sox10 Dom/+ 15.5dpc scRNA-seq data. Branch labels are those designated by Morarach et al, 2021. (B) UMAP showing the dataset overlay after Harmony batch correction color-coded by dataset origin. Arrow highlights differences between datasets. (C) UMAPs displaying integration results of the 15.5-dpc dataset with either the Sox10+/+ (top row) and Sox10Dom/+ (bottom row). Labels by group show the distribution of progenitors, two neuroblast paths, and two (Morarach, left) or three (Sox10+/+ and Sox10Dom/+) neuronal branches. Sox10+/+ Branch C shares UMAP area with Morarach 15.5-dpc Branch B (top), whereas Sox10Dom/+ Branch B shares UMAP area with Morarach 15.5-dpc Branch B. (D) UMAPs of Pou3f3 and Pantr1 expression in Sox10+/+ and Sox10Dom/+ 15.5dpc scRNA-seq data. (E) UMAP of remaining Sox10+/+ cells after removal of all cells that express either Pou3f3, Pantr1, or both. (F) Plots of numbers of UMIs and Features for the Morarach data compared with the current data set shows that significantly more counts in the deeply sequenced 15.5-dpc data (nCount_RNA: unpaired Wilcoxon rank-sum P = 1.060987 × 10-133; nFeature_RNA: unpaired Wilcoxon rank-sum P = 4.927273 × 10-187)
Figure 6
Figure 6
Neuronal-specific markers persist and highlight distinct cell populations throughout ENS development. Feature plots demonstrating pairs of coexpressed markers restricted to each neuronal branch for multiple enteric single-cell data sets across later fetal stages (Morarach et al data at 15.5 dpc and 18.5 dpc) and into adulthood (May-Zhang et al, data at P42-49). Branch A (Islr2 / Cartpt), Branch B (Gch1 / Tac1), and Branch C (Calcb / Nog). ∗May-Zhang et al, 2021 snRNA-seq data only express Cartpt at appreciable levels.
Figure 7
Figure 7
Sox10Dom/+ shifts the distribution of ENS progenitors and differentiating neurons early in ENS development. (A) Separate UMAPs of Sox10+/+ and Sox10Dom/+ progenitor and neuronal clusters. (B) Bar chart displaying cell counts in every cluster between Sox10+/+ and Sox10Dom/+ profiled populations. ∗P value < .05 from Student’s t-test. (C) Relative differences in abundance of cells per cluster between genotypes. Red clusters have an FDR <0.5 and mean Log2 fold enrichment >1 (permutation test, n = 1000). Clusters to the left are under-represented in Sox10Dom/+, while clusters to the right are over-represented in Sox10Dom/+. (D) The DA-seq plot displays a DA of cell populations between genotypes. Cells of increased abundance in Sox10+/+ (blue) compared with those increased in Sox10Dom/+ (red).
Figure 8
Figure 8
Sox10Dom/+-exclusive cluster PC6 exhibits character of noncycling, migrating progenitors. (A) Gene expression of known cell cycle genes for G1, G2/M, and S phases were used in concert with Seurat’s CellCycleScoring function to assign cell cycle phases to each cell, showing 3 progenitor clusters (PC0, PC2, and PC6) as noncycling or in G1 phase. (B) Differential gene expression between PC6, the Sox10Dom/+-exclusive cluster and the 2 other G1/non-cycling progenitor clusters, as these are the most similar progenitor clusters and avoids cell cycle-specific differential gene expression in the comparative results. Genes in red are upregulated in PC6/downregulated in PC0 and PC2, whereas those in blue are downregulated in PC6/upregulated in PC0 and PC2. (C) The top 30 GO terms for significant differentially expressed genes downregulated in PC6/upregulated in PC 0 and PC 2 (Bonferroni-adjusted P value < .05, average log2(fold change) <−1). (D) The top 30 GO terms for significant differentially expressed genes upregulated in PC 6/downregulated in PC 0 and PC 2 (Bonferroni-adjusted P value < .05, average log2(fold change) >1). For (C) and (D), datapoint size corresponds to how many genes are associated with a given ontology term, the datapoints are colored by Bonferroni-adjusted P value from the ontology analysis, and the x-axis is the ratio of genes differentially expressed to the number of genes associated with a given ontology term from the database.
Figure 9
Figure 9
Sox10Dom/+ alters calculated differentiation rates between cell clusters without affecting cell cycle progression. (A) RNA Velocity plots show cluster trajectories for Sox10+/+ and Sox10Dom/+. Two distinct transitional paths, through which progenitors transition towards neuronal fates, are indicated: Path1 (NC1, NC8) and Path2 (NC11, NC10). (B) Topographical plots depict the expression of proliferative markers Pcna, Top2a, Mik67, and Ccnb1 along Ccnb1 transitional Path1 (NC1, NC8) and Path2 (NC11, NC10). (C) Tricycle analysis embedding distribution of cells in each phase of mitosis for different genotypes. Each annotation corresponds to a specific cell cycle state: 0π (G0/G1), 0.5π (S), 1π (G2), 1.5π (M), and 2π (G0/G1). (D) Density plot illustrating the distribution of cells in each phase for Sox10+/+ (blue) and Sox10Dom/+ (red). (E) Cell transition rates calculated via scVelo trajectory analysis for developing ENS clusters of Sox10+/+ and Sox10Dom/+ data are shown in positions comparable to the UMAP cluster positions. Transition rates were calculated per genotype with values shown in each box. Hues of red indicate higher transition rates, whereas lower rates are represented by blue hues relative to the calculated transition rate scale shown at bottom. For each Sox10Dom/+cluster, calculated rates are noticeably lower than those calculated for Sox10+/+. Comparing individual cluster calculated rates across genotypes for example Sox10+/+ NC9 (32.77) vs Sox10Dom/+ NC9 (20.65) illustrates these differences.
Figure 10
Figure 10
Multiple genes associated with neuronal specification in the developing ENS are significantly downregulated in Sox10Dom/+ mutants. (A) Volcano plot displaying the differentially expressed genes identified when comparing all cells from Sox10+/+ and Sox10Dom/+ data. Genes upregulated in Sox10Dom/+ are in red and genes upregulated in Sox10+/+ are shown in blue. Genes with an adjusted P value of 0 were manually changed to be 20-fold less than the lowest nonzero P value so that the labels for these genes can be shown. (B) Individual volcano plots for each neuronal branch within the Sox10Dom/+ data. The genes plotted have an average Log2 Fold Change of >0.5 and <−0.5 and a Bonferroni-adjusted P < .05. Downregulated genes in Sox10Dom/+ are in red, whereas upregulated genes in Sox10Dom/+ are in blue. Similarly, upregulation in Sox10+/+ is in red, whereas downregulation in Sox10+/+ is in blue.
Figure 11
Figure 11
Multiple gene regulatory networks that exhibit elevated activity within developing enteric neuronal trajectories are associated with Hox-related genes. Heatmap of gene regulatory network activity calculated via SCENIC for wildtype Sox10+/+ (A) and Sox10Dom/+ (B) data. Elevated activity is shown as hues of red with decreased activity indicated by hues of blue for each gene-associated regulon. Arrows emphasize Hox genes whose expression shifts markedly between genotypes.
Figure 12
Figure 12
Hox transcripts associated with active regulatory networks are differentially expressed between Sox10+/+ and Sox10Dom/+ ENS populations. (A) Feature plot UMAPs display Hoxa6 expression patterns in Sox10+/+ and Sox10Dom/+ single-cell data. A maximum normalized expression threshold was set to 0.5 in these images to highlight most Hoxa6 expression. (B) Fluorescent confocal images of HCR in situ stained for HuC/D in 15.5-dpc distal ileum from Sox10+/+ and Sox10Dom/+ fetal gut samples are shown. Insets, from regions with dim HuC/D cells indicated by arrowheads, emphasize intensity of Hoxa6 HCR signal in dim HuC/D+ cells. (C) Boxplot displaying positive, Log10-transformed HCR signal ratios (Hoxa6 / Elavl4) in wholemount fetal ileum samples based on all HuC/D+ cells (100%) to the dimmest quartile of HuC/D+ cells (25%). Comparison of Hoxa6 / Elavl4 HCR ratios based on intensity thresholds between Sox10+/+ vs Sox10Dom/+ genotype revealed a significant difference for all HuC/D+ cells (P = 1.46e-08), the 75% dimmest (P = 2.20e-19), the 50% (P = 4.82e-23) and 25% dimmest HuC/D+ cells (P = 1.06e-29). Asterisk indicates Wilcoxon nonpaired test P values between genotype (P < .05). Each datapoint represents HCR signal ratio of one ROI in a given image. (D) Violin plots show the normalized Log2-fold expression levels of Hox genes in Sox10+/+ (blue) vs Sox10Dom/+ (red) cells that depict how the Sox10Dom/+ allele alters expression levels of Hoxa6 within Path2 and Branch C, whereas other Hox genes are more broadly altered across other groups and are not restricted to Branch C.
Figure 13
Figure 13
Observations and model of Sox10-mediated enteric neurogenesis. As early as 15.5 dpc, distinct ENS progenitors, 2 neuroblast transition pathways, and 3 distinct neuronal trajectories can be identified via deep scRNA-seq. Disruption of ENS development by the Sox10Dom/+ mutant allele produces a novel progenitor population and causes a noticeable shift in abundance of cells between neuronal trajectories with decreased cell abundance in Branch C and an increase in Branch B (bold arrows). Downregulation of genes associated with neuronal specification are detected, likely contributing to the calculated decrease in transition rates between cell states.
See this image and copyright information in PMC

Update of

References

    1. Sharkey K.A., Mawe G.M. The enteric nervous system. Physiol Rev. 2023;103:1487–1564. - PMC - PubMed
    1. Drokhlyansky E., Smillie C.S., Van Wittenberghe N., et al. The human and mouse enteric nervous system at single-cell resolution. Cell. 2020;182:1606–1622.e23. - PMC - PubMed
    1. May-Zhang A.A., Tycksen E., Southard-Smith A.N., et al. Combinatorial transcriptional profiling of mouse and human enteric neurons identifies shared and disparate subtypes in situ. Gastroenterology. 2021;160:755–770.e26. - PMC - PubMed
    1. Morarach K., Mikhailova A., Knoflach V., et al. Diversification of molecularly defined myenteric neuron classes revealed by single-cell RNA sequencing. Nat Neurosci. 2021;24:34–46. - PMC - PubMed
    1. Lake J.I., Heuckeroth R.O. Enteric nervous system development: migration, differentiation, and disease. Am J Physiol Gastrointest Liver Physiol. 2013;305:G1–G24. - PMC - PubMed

MeSH terms

Substances

LinkOut - more resources