Enteric nervous system development relies on intestinal colonization by enteric neural crest-derived cells (ENCDCs). This is driven by a population of highly migratory and proliferative ENCDCs at the wavefront, but the molecular characteristics of these cells are unknown. ENCDCs from the wavefront and the trailing region were isolated and subjected to RNA-seq. Wavefront-ENCDCs were transcriptionally distinct from trailing ENCDCs, and temporal modelling confirmed their relative immaturity. This population of ENCDCs exhibited altered expression of ECM and cytoskeletal genes, consistent with a migratory phenotype. Unlike trailing ENCDCs, the wavefront lacked expression of genes related to neuronal or glial maturation. As wavefront ENCDC genes were associated with migration and developmental immaturity, the genes that remain expressed in later progenitor populations may be particularly pertinent to understanding the maintenance of ENCDC progenitor characteristics. Dusp6 expression was specifically upregulated at the wavefront. Inhibiting DUSP6 activity prevented wavefront colonization of the hindgut, and inhibited the migratory ability of post-colonized ENCDCs from midgut and postnatal neurospheres. These effects were reversed by simultaneous inhibition of ERK signaling, indicating that DUSP6-mediated ERK inhibition is required for ENCDC migration in mouse and chick.

Enteric neural crest-derived cells (ENCDCs) migrate along the length of the developing intestine to form the enteric nervous system (ENS). These cells are highly migratory and travel long distances during gut development. The impressive migratory potential of ENCDCs relies on specialized mechanisms that, when disrupted, result in failure of intestinal colonization and ENS formation. The most characterized congenital disorder of ENS malformation is Hirschsprung disease (HSCR), where ENCDCs fail to complete their craniocaudal migration of the gut. The vast number of mutations that have been identified to contribute to this multifactorial congenital disorder highlights the complexity of ENCDC migration (Kuil et al., 2021; Mederer et al., 2020).

During mouse embryogenesis, vagal level neural crest cells (NCCs) enter the foregut mesenchyme (9.5 days post-coitum in mice), when they are termed ENCDCs, and advance caudally to the distal end of the hindgut (Heanue and Pachnis, 2007). Colonization of the intestinal tract relies on several interactions with the pre-colonized gut, including GDNF- and EDN3-mediated paracrine signaling from the mesenchyme (Mwizerwa et al., 2011; Nagy and Goldstein, 2006; Uesaka and Enomoto, 2010), and physical interactions with the surrounding cellular and extracellular environment (Graham et al., 2017; Nagy et al., 2009). The speed of ENCDC migration in the developing gut is particularly pertinent, as colonization is partially dependent on the permissiveness of the tissue environment, which declines with maturation (Druckenbrod and Epstein, 2009; Hotta et al., 2010). Interestingly, key ligands, such as GDNF, can have multiple roles in migrating ENCDCs by being chemoattractive and promoting neural differentiation, with an excess of the latter arresting colonization (Mwizerwa et al., 2011; Nagy et al., 2020, 2021). The mechanisms within ENCDC populations that modulate differential responses to these exogenous signals for appropriate colonization are currently unknown.

The colonization of the gut is driven by the most caudally migrated ENCDCs (also known as the migratory wavefront) in a process of frontal expansion (Nishiyama et al., 2012; Sidebotham et al., 2002; Simpson et al., 2007). ENCDCs behind the wavefront (trailing) are relatively stationary and do not contribute to ENS expansion beyond their initial point of migration (Druckenbrod and Epstein, 2007; Nishiyama et al., 2012; Simpson et al., 2007). This suggests that wavefront ENCDCs may possess unique properties that facilitate invasion of the gut. Indeed, several studies have noted distinguishing features between the wavefront and trailing cells, including high levels of cell proliferation and a lack of neural differentiation (Nagy et al., 2012; Olden et al., 2008; Simpson et al., 2007; Young et al., 2014, 1999).

Although several studies show that the gut is predominately colonized by ENCDCs at the wavefront, it is unclear whether this population exhibits specialized properties compared with trailing ENCDCs. Many ENCDCs remain migratory, albeit to a lesser extent, for 24 h post-colonization in mice (Young et al., 2014), and transplantation of trailing ENCDCs to the wavefront can result in frontal expansion and colonization, as demonstrated in chick-quail chimera experiments (Simpson et al., 2007). Notably, NCC-derived progenitors with migratory potential can be isolated from the late fetal gut and into adulthood (Hotta et al., 2013). Therefore, it is unclear whether ENCDCs at the wavefront exhibit specialized features for migration and ENS formation compared with trailing ENCDCs and other progenitor ENCDC populations at later developmental stages.

In this study, transcriptomic analysis is performed to investigate distinguishing features between the wavefront and trailing ENCDCs that may explain their invasive capacity and ENS-forming capabilities, including developmental maturity and promigratory mechanisms. Furthermore, the transcriptomic features that persist from wavefront ENCDCs to progenitors at later development are explored.

ENCDCs at the colonizing wavefront exhibit a unique gene signature associated with promigratory mechanisms

Complex developmental mechanisms are responsible for the migration of ENCDCs through the gut. According to Young et al. (2014), the speed and directional migration of ENCDCs between 750 and 1500 μm from the wavefront are significantly slower than those within 450 μm of the most-caudal cell. Therefore, we sought to investigate whether wavefront ENCDCs exhibit specialized properties associated with migration by profiling these cell populations by RNA-seq. For these experiments, the Wnt1-tdT reporter mouse line was used at 11.5 days post coitum (dpc) (Fig. 1A). At this developmental stage, ENCDCs have densely colonized the gut proximal to the ileocecal junction and are still actively migrating as they begin their colonization of the nascent cecum and proximal hindgut (Fig. 1B). Furthermore, at 11.5 dpc the gut becomes large enough to distinguish and dissect 0.5 mm of the gut containing ‘wavefront ENCDCs’ (0-0.5 mm proximal in the wavefront) and ‘trailing ENCDCs’ (1-1.5 mm proximal to the wavefront) for downstream experiments (Fig. 1B). As demonstrated in Wnt1-tdT; NestinGFP mice, ENCDCs throughout the intestine express the stem cell marker nestin (Nes), regardless of distance from the migratory wavefront (Fig. 1C). Intestinal tissues from the wavefront and trailing regions were enzymatically digested in collagenase-dispase solution to obtain single cell suspensions. Wnt1-tdT+ cells from both groups were isolated by flow cytometry and RNA was immediately prepared from pooled flow-sorted ENCDCs for transcriptomic analysis by bulk RNA-seq (Fig. 1D). The number of trailing ENCDCs was fourfold higher than wavefront ENCDCs, which reflects the lower density of ENCDCs during active colonization (Fig. 1D). Principal component analysis (PCA) indicated that wavefront and trailing ENCDCs have distinct transcriptional profiles that are dependent on their location along the migration axis of the gut (Fig. 1E). Differential gene expression analysis revealed 49 genes upregulated and 93 genes downregulated at the wavefront (Fig. 1F, Table S1). A protein-protein interaction (PPI) network analysis was performed using the differentially expressed genes (DEGs) between wavefront and trailing ENCDCs to group the properties of gene products with predicted interactions (Fig. 1G). Five key interaction modules were identified within the genes upregulated at the wavefront and these were primarily related to the ECM and cytoskeletal structure, consistent with a promigratory role in these cells. Genes downregulated at the wavefront were associated with growth cones and neuronal cell bodies, reflecting a lack of neural differentiation relative to the trailing ENCDCs.

Fig. 1.

Wavefront ENCDCs are transcriptionally unique and exhibit a gene expression profile associated with promigratory pathways. (A) Whole mouse embryo at 11.5 dpc expressing Wnt1-tdT. Scale bar: 2 mm. (B) Partially dissected mouse embryo revealing the hairpin-loop structure of the developing gut. ENCDCs migrating in a caudal direction characteristically reach the ileocecal junction (arrow). Scale bar: 100 µm. (C-Cʹʹʹ) Dissected gut at 11.5 dpc from NestinGFP;Wnt1-tdT embryos. (C) Wnt1-tdT expression in ENCDCs is observed caudally up to the ileocecal junction (IC) with wavefront ENCDCs considered to be within 0.5 mm proximal to the wavefront and ENCDCs 1.0-1.5 mm posterior to the migrating wavefront referred to as trailing ENCDCs. (C′-C′′′) NestinGFP expression is observed in cells of the mesenteric blood vessel (C′) and colocalized with Wnt1-tdT, which is characteristic of enteric neural progenitor cells (C′′,C′′′) in the wavefront and trailing ENCDCs. Scale bars: 100 µm. (D) FACS plots of Wnt1-tdT and DAPI (dead cell) expression in single-cell suspensions generated from regions of intestine containing the wavefront and trailing ENCDCs. (E) Principal component plot of gene expression data from individual samples of wavefront and trailing ENCDCs. (F) Volcano plot of differentially expressed genes between wavefront (WF) and trailing ENCDCs from RNA-sequencing of purified ENCDC populations. (G) Node-edge network of protein-protein interactions and clustering between differentially expressed genes in wavefront and trailing ENCDCs. Interacting clusters are annotated for top associated gene ontology terms.

Fig. 1.

Wavefront ENCDCs are transcriptionally unique and exhibit a gene expression profile associated with promigratory pathways. (A) Whole mouse embryo at 11.5 dpc expressing Wnt1-tdT. Scale bar: 2 mm. (B) Partially dissected mouse embryo revealing the hairpin-loop structure of the developing gut. ENCDCs migrating in a caudal direction characteristically reach the ileocecal junction (arrow). Scale bar: 100 µm. (C-Cʹʹʹ) Dissected gut at 11.5 dpc from NestinGFP;Wnt1-tdT embryos. (C) Wnt1-tdT expression in ENCDCs is observed caudally up to the ileocecal junction (IC) with wavefront ENCDCs considered to be within 0.5 mm proximal to the wavefront and ENCDCs 1.0-1.5 mm posterior to the migrating wavefront referred to as trailing ENCDCs. (C′-C′′′) NestinGFP expression is observed in cells of the mesenteric blood vessel (C′) and colocalized with Wnt1-tdT, which is characteristic of enteric neural progenitor cells (C′′,C′′′) in the wavefront and trailing ENCDCs. Scale bars: 100 µm. (D) FACS plots of Wnt1-tdT and DAPI (dead cell) expression in single-cell suspensions generated from regions of intestine containing the wavefront and trailing ENCDCs. (E) Principal component plot of gene expression data from individual samples of wavefront and trailing ENCDCs. (F) Volcano plot of differentially expressed genes between wavefront (WF) and trailing ENCDCs from RNA-sequencing of purified ENCDC populations. (G) Node-edge network of protein-protein interactions and clustering between differentially expressed genes in wavefront and trailing ENCDCs. Interacting clusters are annotated for top associated gene ontology terms.

The gene expression profile of wavefront ENCDCs was further explored by gene set enrichment analysis (GSEA), a computational method used to determine statistical overrepresentation of defined gene sets within an input list of ranked genes between conditions (Fig. 2A). Hallmark gene sets were used to identify cellular processes that might promote the promigratory properties of ENCDCs. The sign(log2FC)×log10(P-value) ranking method was implemented for gene expression between wavefront and trailing ENCDCs. The sign function provides direction to gene regulation (+ or −) in gene lists weighted by P-values. Wavefront ENCDCs were enriched for functional processes associated with cell proliferation, including G2M checkpoint and E2F targets (Fig. 1A, Table S2), which reaffirms the highly proliferative nature of wavefront ENCDCs, as previously reported (Simpson et al., 2007). Wavefront ENCDCs were also associated with a gene signature related to epithelial-to-mesenchymal transition (EMT) (Fig. 1A). This process is essential for early neural crest cell (NCC) delamination from the dorsal tube, and these data indicate that similar properties may be maintained during ENCDC colonization of the gut. Top EMT-associated genes significantly upregulated at the wavefront included Fbn2, Itgb1, Tnc, Fbln1, Lgals1, Cald1, Flna, Lum, Wnt5a, Fn1, Col1a1, Col3a1, Cdh11 and Tgfbi. Of these genes, fibronectin (Fn1), tenascin-C (Tnc) and integrin β1 (Itgb1) have previously been shown to be crucial to ENCDC colonization (Akbareian et al., 2013; Breau et al., 2006; Watanabe et al., 2013). Expression of cadherin 11 and fibrillin 2 was further confirmed at the protein level by immunohistochemistry in Wnt1-tdT mice at 11.5 dpc (Fig. 2B,C). Cadherin 11 contributes to NCC migration in Xenopus by complexing with the fibronectin receptor syndecan 4, mediating substrate adhesion (Langhe et al., 2016). The syndecan 4 gene Sdc4 was highly expressed in the wavefront and trailing ENCDCs in our dataset, and its expression was confirmed by immunohistochemistry (Fig. 2D). Although syndecan 4 and fibrillin 2 remain expressed at 13.5 dpc in the proximal hindgut, they were not expressed as uniformly by ENCDCs as at the 11.5 dpc wavefront (Fig. S1).

Fig. 2.

Genes associated with epithelial-mesenchymal transition are upregulated in wavefront ENCDCs. (A) GSEA analysis of wavefront versus trailing ENCDC gene expression against the HALLMARK pathway database. EMT, epithelial-mesenchymal transition. Data are presented as normalized enrichment scores (NES, y axis), leading edge genes (size) and -Log10 P value (color). (B-D) Representative images of (B) cadherin 11 (CDH11), (C) fibrillin 2 (FBN2) and (D) syndecan 4 (SYND4) expression by immunohistochemistry in Wnt1-tdT-expressing mouse ENCDCs at the wavefront at 11.5 dpc. Areas outlined in B and D are shown at higher magnification on the right. Scale bars: 100 μm in B and D; 20 μm in B (inset), C and D (inset).

Fig. 2.

Genes associated with epithelial-mesenchymal transition are upregulated in wavefront ENCDCs. (A) GSEA analysis of wavefront versus trailing ENCDC gene expression against the HALLMARK pathway database. EMT, epithelial-mesenchymal transition. Data are presented as normalized enrichment scores (NES, y axis), leading edge genes (size) and -Log10 P value (color). (B-D) Representative images of (B) cadherin 11 (CDH11), (C) fibrillin 2 (FBN2) and (D) syndecan 4 (SYND4) expression by immunohistochemistry in Wnt1-tdT-expressing mouse ENCDCs at the wavefront at 11.5 dpc. Areas outlined in B and D are shown at higher magnification on the right. Scale bars: 100 μm in B and D; 20 μm in B (inset), C and D (inset).

Our ENCDC data were further analyzed to determine whether the transcriptional features between wavefront and trailing ENCDCs are consistent with other migrating NCC-derived cells. GSEA was performed on ranked gene expression of our ENCDC data using signatures described previously in the wavefront and trailing axis of the cranial NCC-derived cell migratory stream during invasion of the paraxial mesoderm (HH13) and entry into the second branchial arch (HH15) in chick embryos (Morrison et al., 2017) (Table S2). Despite the limitations of cross-species comparisons, such as non-homologous gene pairings, wavefront ENCDCs had a significant negative association (bottom-ranked genes), with the trailing signature of migrating cranial NCCs in the chick (Table S2). This included an increase in expression of Plp1 and Hoxd3 in trailing cells. Conversely, wavefront ENCDCs were significantly associated with the wavefront of cranial NCCs of the paraxial mesoderm, suggesting that features of the wavefront-trailing cell axis are shared between species and migratory NCC streams.

Gene signatures of the wavefront-trailing ENCDC axis associate with developmental maturity and cell fate progression

Despite existing at the same developmental stage, we hypothesized that wavefront ENCDCs represent a transcriptionally less mature cell population than trailing ENCDCs. This was tested using existing datasets of single-cell RNA-seq gene expression obtained from mice at 9.5 to 13.5 dpc from Cao et al. (2019). NCC-derived cells were analyzed from the glia and neuron developmental trajectories in the peripheral nervous system (PNS) at these stages (Fig. 3A). Cell clustering was consistent with the original annotations of the authors and included clusters specific for enteric neuron and enteric glia trajectories (Fig. 3B and Fig. S2A). In this dataset, NCC developmental trajectories are derived from single-cell sequencing of the entire mouse. To confirm these trajectories were specific to the gut, we preformed data integration specifically with ENCDCs from mice at 15.5-18.5 dpc using single-cell RNA-seq data produced by Morarach et al. (2021) (Fig. 3C, Fig. S2B). Enteric neuron trajectory_1 and enteric glia and Schwann cell trajectory_2 from the Cao et al. (2019) dataset (Fig. 3D) were the only clusters to closely integrate with the Morarach et al. (2021) data (Fig. 3E). These cells shared transcriptional similarity with enteric glia/progenitors, neuroblasts and enteric neurons, confirming their cell identity as ENCDCs. Enteric neuron trajectory_2 and enteric glia and Schwann cell trajectory_1 from the Cao et al. (2019) dataset were omitted from analysis as we were unable to confirm their transcriptional similarity to ENCDCs provided by the Morarach et al. (2021) dataset (Fig. S2C). The developing enteric neuron or enteric glia clusters, and the adjacent NCC progenitor clusters that were enriched at 9.5 dpc, were taken for subsequent analysis (Fig. 3F). Dimensionality reduction shows gradual deviation away from the cell clusters enriched at 9.5 dpc as developmental stage progresses (Fig. 3G). To determine whether the transcriptomic features observed spatially between wavefront and trailing ENCDCs in our 11.5 dpc data represent a developmental axis of maturity, we scored the enteric neuron and enteric glia development trajectories dataset for wavefront and trailing signatures. This was performed by creating gene sets from the DEGs upregulated (wavefront signature) or downregulated (trailing signature) between wavefront and trailing ENCDCs. The AddModuleScore function of the Seurat package was used to score the expression of these signatures compared with computed sets of equivalent control genes to each cell individually using methods derived from Tirosh et al. (2016). This provides an overall estimation of wavefront and trailing signatures that account for differences in the complexity and quality of individual expression data between cells. In both the developing enteric neuron and glia data, we observed an increase in the ‘wavefront’ and ‘trailing’ signature scores over developmental stages (Fig. 3H). Differential gene expression analysis was performed for ‘wavefront’ and ‘trailing’ genes between developmental stages (Table S3). Gene expression of the DEGs in the neuronal and glial trajectory largely paralleled the direction of ‘wavefront’ and ‘trailing’ scores. In the neuron trajectory, 67% of DEGs specific to the ‘wavefront’ and detected at 9.5 dpc were downregulated at later stages, and 83% of the ‘trailing’ DEGs detected were upregulated at 13.5 dpc. In the glial trajectory, 86% of DEGs specific to the ‘wavefront’ and detected at 9.5 dpc were downregulated at later stages, and 80% of the ‘trailing’ DEGs detected were upregulated at 13.5 dpc. To further explore how the ‘wavefront’ and ‘trailing’ ENCDC-associated genes change over time, Monocle was used to model developmental pseudotime and detect corresponding differential expression in genes (Fig. 3I). Pseudotime predictions followed a similar pattern to actual developmental staging (Fig. 3G). The expression of several ‘wavefront’ and ‘trailing’ ENCDC-associated genes was significantly altered over time and typically followed a decreasing trend for ‘wavefront’ genes and an increasing trend for ‘trailing’ genes (Fig. 3J,K, Table S4). Together, this indicates that 11.5 dpc wavefront ENCDCs express numerous genes associated with developmental immaturity, while trailing ENCDCs exhibit gene signatures associated with later developmental stages. The wavefront-trailing cell axis may thus be considered as a developmental pseudotime that reflects the immaturity of wavefront ENCDCs and their gradual acquisition of maturity post-colonization at the same developmental stage. To explore this further, we repeated developmental pseudotime analysis in the ENCDCs of the neuron and glia developmental trajectories, this time modelling ENCDC progenitors in transition to enteric neuron commitment using dimensionality reduction learned from the Morarach et al. (2021) dataset (Fig. 3L). In this analysis, wavefront genes are more highly expressed by the cycling progenitor cells, which are downregulated as the ENCDCs develop into neurons (Fig. 3M, Table S4). Likewise, the ENCDCs gain the expression signature of trailing ENCDCs during development into neurons. This suggests that wavefront-associated genes are lost after commitment to neurons, which may contribute to the inability of committed neurons to migrate, as previously observed (Young et al., 2004, 2014).

Fig. 3.

Assessment of wavefront-trailing ENCDC developmental maturity and conserved transcriptomic features in migrating NCCs. (A) UMAP plot of peripheral nervous system (PNS) neural crest cells at developmental timepoints 9.5 to 13.5 dpc from single-cell RNA-seq data originally generated by Cao et al. (2019). (B) Unsupervised clustering of the Cao et al. dataset represented on UMAP. (C) UMAP representation of enteric neuron and enteric glial developmental trajectories from the Cao et al. dataset integrated with single cell RNA-seq data of ENCDCs at 15.1-18.5 dpc generated by Morarach et al. (2021). (D) Original location of enteric neuron (blue) and enteric glia (red) developmental trajectories in the Cao et al. dataset. (E) Overlay of enteric neuron (blue) and enteric glia (red) developmental trajectories from Cao et al. after integration with the Morarach et al. dataset. (F) UMAP representation of the Cao et al. PNS neural crest cell dataset highlighting clusters associated with enteric neuron (dark blue) and enteric glia (dark red) developmental trajectories and clusters of their closest associated neural crest cell progenitors (neuron, light blue; glia, light red) enriched at 9.5 dpc. (G) UMAP plots of subclusters comprising the enteric neuron trajectory and neural crest cell progenitors (top panel), and the enteric glia trajectory and neural crest cell progenitors (bottom panel). (H) Bar chart representations of average wavefront and trailing ENCDC gene set signature scores generated by the AddModuleScore function of Seurat in cells from the enteric neuron and enteric glial development subclusters. Data are presented as average±s.e.m. per actual developmental staging of cells. Kruskal–Wallis ANOVA test with a post-hoc Dunn's multiple comparison test. *P<0.05, **P<0.01 and ****P<0.0001. (I) Developmental pseudotime estimation of cell maturation in the enteric neuron and enteric glia development subclusters. (J,K) Heatmap representations of significant changes in gene expression over developmental pseudotime for (J, left panel) wavefront ENCDC-associated genes declining over pseudotime maturation in enteric neuron trajectories, (J, right panel) trailing ENCDC-associated genes increasing over pseudotime maturation in enteric neuron trajectories and (K) wavefront and trailing ENCDC-associated genes changing over pseudotime maturation in enteric glial trajectories. Annotations above heatmaps reflect modelled pseudotime, original cell clusters and actual developmental stage when cells were collected. (L) Developmental pseudotime estimation of cell maturation in enteric neuron and enteric glial trajectory clusters from the Cao et al. dataset after integration with the Morarach et al. data to improve developmental modelling. (M) Heatmap representations of significant changes in gene expression over developmental pseudotime for (left panel) wavefront ENCDC-associated genes declining over pseudotime and (right panel) trailing ENCDC-associated genes increasing over pseudotime. Annotations above heatmaps reflect modelled pseudotime, cell identity of clusters derived from Morarach et al., and actual developmental stage when cells were collected.

Fig. 3.

Assessment of wavefront-trailing ENCDC developmental maturity and conserved transcriptomic features in migrating NCCs. (A) UMAP plot of peripheral nervous system (PNS) neural crest cells at developmental timepoints 9.5 to 13.5 dpc from single-cell RNA-seq data originally generated by Cao et al. (2019). (B) Unsupervised clustering of the Cao et al. dataset represented on UMAP. (C) UMAP representation of enteric neuron and enteric glial developmental trajectories from the Cao et al. dataset integrated with single cell RNA-seq data of ENCDCs at 15.1-18.5 dpc generated by Morarach et al. (2021). (D) Original location of enteric neuron (blue) and enteric glia (red) developmental trajectories in the Cao et al. dataset. (E) Overlay of enteric neuron (blue) and enteric glia (red) developmental trajectories from Cao et al. after integration with the Morarach et al. dataset. (F) UMAP representation of the Cao et al. PNS neural crest cell dataset highlighting clusters associated with enteric neuron (dark blue) and enteric glia (dark red) developmental trajectories and clusters of their closest associated neural crest cell progenitors (neuron, light blue; glia, light red) enriched at 9.5 dpc. (G) UMAP plots of subclusters comprising the enteric neuron trajectory and neural crest cell progenitors (top panel), and the enteric glia trajectory and neural crest cell progenitors (bottom panel). (H) Bar chart representations of average wavefront and trailing ENCDC gene set signature scores generated by the AddModuleScore function of Seurat in cells from the enteric neuron and enteric glial development subclusters. Data are presented as average±s.e.m. per actual developmental staging of cells. Kruskal–Wallis ANOVA test with a post-hoc Dunn's multiple comparison test. *P<0.05, **P<0.01 and ****P<0.0001. (I) Developmental pseudotime estimation of cell maturation in the enteric neuron and enteric glia development subclusters. (J,K) Heatmap representations of significant changes in gene expression over developmental pseudotime for (J, left panel) wavefront ENCDC-associated genes declining over pseudotime maturation in enteric neuron trajectories, (J, right panel) trailing ENCDC-associated genes increasing over pseudotime maturation in enteric neuron trajectories and (K) wavefront and trailing ENCDC-associated genes changing over pseudotime maturation in enteric glial trajectories. Annotations above heatmaps reflect modelled pseudotime, original cell clusters and actual developmental stage when cells were collected. (L) Developmental pseudotime estimation of cell maturation in enteric neuron and enteric glial trajectory clusters from the Cao et al. dataset after integration with the Morarach et al. data to improve developmental modelling. (M) Heatmap representations of significant changes in gene expression over developmental pseudotime for (left panel) wavefront ENCDC-associated genes declining over pseudotime and (right panel) trailing ENCDC-associated genes increasing over pseudotime. Annotations above heatmaps reflect modelled pseudotime, cell identity of clusters derived from Morarach et al., and actual developmental stage when cells were collected.

The analysis of wavefront and trailing genes in the Cao et al. (2019) dataset appeared to suggest that wavefront ENCDCs exhibit fewer features associated with cell fate commitment compared with trailing ENCDCs. Therefore, we sought to assess the association of enteric neuronal and glial gene signatures with wavefront and trailing ENCDCs. This was performed using GSEA with gene sets of human enteric neuron and glia signatures obtained from MSigDB to examine their association with the top- and bottom-ranked genes between wavefront and trailing ENCDCs (Fig. 4A, Table S5). The bottom-ranked genes between wavefront and trailing ENCDCs were associated with both enteric glia and enteric neuron signatures, as signified by negative enrichment scores (Fig. 4A). This suggests that trailing ENCDCs become more neuronal and glial like after colonization, and this included upregulation of Tau (Mapt) and Plp1, which are enteric neuron and glial markers, respectively, according to the MSigDB annotations (Fig. 4B) and shown in previous reports (Nagy et al., 2018; Rao et al., 2015). Nevertheless, evidence of greater glial differentiation in the trailing ENCDC population may not necessarily equate to a lack of progenitor properties in ENCDCs as there is a lack of clear transcriptional segregation between populations of progenitors and glia in the embryonic gut environment, and several studies observe that neurons still arise from glial cells postnatally (Guyer et al., 2023; Laranjeira et al., 2011).

Fig. 4.

Wavefront ENCDCs are transcriptionally distinct from developing neurons and glia. (A) Enrichment of enteric neuron and enteric glia signature genes from MSigDB by GSEA in 11.5 dpc wavefront versus trailing ENCDC gene expression and postnatal ENCDC progenitors. Data are presented as normalized enrichment scores (NES, y-axis), leading-edge genes (size) and -Log10 P value (color). (B) Heatmap of leading-edge genes (GSEA) of enteric neuron and enteric glia signatures from MSigDB annotations in 11.5 dpc wavefront versus trailing ENCDC gene expression profiles presented as Z-scores. (C) Validation of TauGFP specificity to enteric neurons in the myenteric plexus of postnatal mice by labeling with the neuronal markers HuC/D and Tuj1. (D) Whole-embryo expression of TauGFP; Wnt1-tdT. Scale bar: 2 mm. (E) Representative images of Wnt1-tdT, TauGFP and merged images at 11.5 dpc. (F) Higher magnification of the area outlined in E showing trailing ENCDCs expressing TauGFP. Arrows indicate TauGFP/Wnt1;tdT+ cells. (G) Validation of Plp1GFP specificity to enteric glia in the myenteric plexus of postnatal mice by labeling with the neuronal marker HuC/D and glial marker S100B. Arrowheads indicate Plp1GFP/S100B+ cells. (H-K) Representative images of cross-sections of the embryonic gut of Wnt1-tdT; Plp1GFP mice at 11.5 dpc at the level of the migrating wavefront ENCDCs (J) and high-magnification images of Tuj1 (I) and Plp1GFP (K). Area outlined in H is shown at higher magnification in I and J. Arrows indicate Wnt1;tdT+ cells. (L) Cross-section from the 11.5 dpc gut of Wnt1-tdT; Plp1-GFP mice showing trailing ENCDCs in the colonized midgut. Area outlined in L is magnified in L′, demonstrating the glial marker Plp1-GFP in trailing ENCCs. (M) Hindgut at 13.5 dpc from a Plp1GFP mouse embryo stained for p75 to mark the ENCDCs. Area outlined in M is shown at higher magnification in M′. (N,O) Consecutive transverse sections of 13.5 dpc Plp1GFP mouse mid-hindgut stained using Tuj1 (N) and Hu (O) antibodies to detect early neuronal differentiation. (P) Quantification of the mean fluorescence intensity of Plp1GFP expression in Wnt1-tdT ENCDCs in the mouse embryonic gut at 11.5 dpc split by distance from the most caudally migrated wavefront ENCDCs in 100 µm increments. All pairwise comparisons are statistically significant at P<0.0001 unless indicated otherwise in the graph. *P<0.05; NS, not significant. Scale bars: 100 μm in C and H; 200 μm in E; 40 μm in G and I-O.

Fig. 4.

Wavefront ENCDCs are transcriptionally distinct from developing neurons and glia. (A) Enrichment of enteric neuron and enteric glia signature genes from MSigDB by GSEA in 11.5 dpc wavefront versus trailing ENCDC gene expression and postnatal ENCDC progenitors. Data are presented as normalized enrichment scores (NES, y-axis), leading-edge genes (size) and -Log10 P value (color). (B) Heatmap of leading-edge genes (GSEA) of enteric neuron and enteric glia signatures from MSigDB annotations in 11.5 dpc wavefront versus trailing ENCDC gene expression profiles presented as Z-scores. (C) Validation of TauGFP specificity to enteric neurons in the myenteric plexus of postnatal mice by labeling with the neuronal markers HuC/D and Tuj1. (D) Whole-embryo expression of TauGFP; Wnt1-tdT. Scale bar: 2 mm. (E) Representative images of Wnt1-tdT, TauGFP and merged images at 11.5 dpc. (F) Higher magnification of the area outlined in E showing trailing ENCDCs expressing TauGFP. Arrows indicate TauGFP/Wnt1;tdT+ cells. (G) Validation of Plp1GFP specificity to enteric glia in the myenteric plexus of postnatal mice by labeling with the neuronal marker HuC/D and glial marker S100B. Arrowheads indicate Plp1GFP/S100B+ cells. (H-K) Representative images of cross-sections of the embryonic gut of Wnt1-tdT; Plp1GFP mice at 11.5 dpc at the level of the migrating wavefront ENCDCs (J) and high-magnification images of Tuj1 (I) and Plp1GFP (K). Area outlined in H is shown at higher magnification in I and J. Arrows indicate Wnt1;tdT+ cells. (L) Cross-section from the 11.5 dpc gut of Wnt1-tdT; Plp1-GFP mice showing trailing ENCDCs in the colonized midgut. Area outlined in L is magnified in L′, demonstrating the glial marker Plp1-GFP in trailing ENCCs. (M) Hindgut at 13.5 dpc from a Plp1GFP mouse embryo stained for p75 to mark the ENCDCs. Area outlined in M is shown at higher magnification in M′. (N,O) Consecutive transverse sections of 13.5 dpc Plp1GFP mouse mid-hindgut stained using Tuj1 (N) and Hu (O) antibodies to detect early neuronal differentiation. (P) Quantification of the mean fluorescence intensity of Plp1GFP expression in Wnt1-tdT ENCDCs in the mouse embryonic gut at 11.5 dpc split by distance from the most caudally migrated wavefront ENCDCs in 100 µm increments. All pairwise comparisons are statistically significant at P<0.0001 unless indicated otherwise in the graph. *P<0.05; NS, not significant. Scale bars: 100 μm in C and H; 200 μm in E; 40 μm in G and I-O.

As a comparator of the GSEA method to determine enteric neuronal and enteric glial signatures, the same analysis was performed on a single-cell RNA-seq dataset of enteric neurospheres derived from the intestine of postnatal mice produced by Guyer et al. (2023). In this dataset, cells within the neurosphere expressing Plp1GFP, as an indicator of glial cell origin, were analyzed and show clear differentiation into neuronal trajectories as a model of postnatal neurogenesis (Fig. S2D). GSEA was performed as above on gene lists ranked by sign(log2FC)×log10(P-value) values between individual clusters versus the remaining dataset (Fig. 4A, Table S6). Postnatal neuroblast and neuronal clusters positively associated with the enteric neuron signature and, conversely, postnatal progenitors were positively associated with the glial signature (Fig. 4A). These data suggest that postnatal progenitors from cultured enteric neurospheres arise from glia and therefore share a similar transcriptional profile (Guyer et al., 2023; Joseph et al., 2011).

The distinctions between wavefront ENCDCs with trailing ENCDCs and postnatal neurosphere progenitors further highlight that wavefront cells appear to exist in a relatively undifferentiated cell state. To validate these data, we assessed the expression of GFP driven by Mapt (Tau) and Plp1 in transgenic Wnt1-tdT dual reporter mice. Notably, GFP fluorescence parallels gene expression data in NestinGFP;Wnt1-tdT mice and is uniformly expressed between wavefront and trailing ENCDCs (Fig. 1C). Expression of TauGFP was validated to be exclusive to enteric neurons postnatally by co-immunolabeling with the neuronal markers HuC/D and Tuj1 (TUBB3) (Fig. 4C). Embryos of TauGFP; Wnt1-tdT mice were dissected at 11.5 dpc (Fig. 4D). In the gut, TauGFP fluorescence was limited in wavefront ENCDCs, whereas several ENCDCs exhibited strong expression in the trailing region (Fig. 4E,F). Expression of Plp1GFP was confirmed to be expressed exclusively by glia in the postnatal gut by colocalization with the glial marker S100B and the absence of HuC/D (Fig. 4G). Plp1GFP expression in the postnatal gut was further validated to be expressed only by ENCDCs (Fig. 4H). Expression of GFP driven by the glial cell marker Plp1 in the wavefront ENCDCs was assessed in embryos of Plp1GFP; Wnt1-tdT mice. At the wavefront, Plp1GFP expression was low or absent in ENCDCs (Fig. 4H). Interestingly, Tuj1 could be observed in both Plp1GFP-positive and -negative cells at the wavefront (Fig. 4I,J); however, there was no overlap between the expression of Plp1GFP and neuronal marker HuC/D at 11.5 dpc (Fig. 4K) or later at 13.5 dpc (Fig. 4M-O), indicating that Plp1GFP is not expressed in maturing enteric neurons. Plp1GFP signal was increased in trailing ENCDCs compared with those at the wavefront at 11.5 dpc (Fig. 4L,L′). Quantification of the fluorescence intensity from individual ENCDCs shows a gradual increase in Plp1GFP expression starting 200 µm behind the wavefront (Fig. 4P).

Wavefront-trailing ENCDC gene signatures persist in progenitor and restricted cell populations after colonization

As the wavefront-trailing ENCDC gene expression signatures appear to reflect the progression of developmental maturity, we investigated whether these signatures could be correlated with progenitor cell populations and committed cells at later stages of development (Fig. 5A). Wavefront and trailing ENCDC signatures were examined at the single-cell level by determining module scores for individual cells in the Morarach et al. (2021) single-cell RNA-Seq dataset of ENCDC development at 15.5-18.5 dpc that includes progenitor, neuroblast and neuron populations (Fig. 5A). Scores for the wavefront signature were highest in static and cycling progenitors, while the trailing signature scores were lower in progenitor cells and higher in developing neurons (Fig. 5B,C). Genes forming the wavefront and trailing signatures were further explored within differential expression analysis between cell clusters (Table S7). Here, we observe that the majority of wavefront-associated genes that were differentially expressed are upregulated in the progenitor population rather than cells committed to a neuronal fate (Fig. 5D). Likewise, the majority of trailing ENCDC-associated genes detected in differential expression analysis were expressed by committed neurons (Fig. 5E). This included common neuronal markers Elavl4, Uchl1, Mapt and Tubb3, and neuronal specification markers Dbh, Gal, Cartpt and Slc18a3 (Fig. 5E). Several trailing ENCDC-associated genes were also confined to the progenitor cell population in this dataset at 15.5-18.5 dpc (Fig. 5E). This included expression of Plp1, Insc, Rbp1 and Gsta4 in 15.5-18.5 dpc progenitor cells. These genes are also highly expressed in postnatal glia, as identified in expression data generated from the mouse colon atlas originally produced from single-nuclei Seq of adult C57Bl/6 mice by Drokhlyansky et al. (2020) (Fig. 5F). Given that Plp1 expression is low in wavefront ENCDCs, these genes could also reflect a gradual transition of early embryonic ENCDCs to a more mature glial-like transcriptional profile.

Fig. 5.

Features of the wavefront-trailing axis of ENCDC maturation are conserved in progenitors at 15.5-18.5 dpc. (A) Unsupervised clustering of the Morarach et al. (2021) dataset containing Wnt1-tdT expressing ENCDCs at 15.5 and 18.5 dpc represented as UMAP. (B,B′) UMAP plot with individual cells colored by module scores generated by the AddModuleScore function of Seurat of the 11.5 dpc wavefront ENCDC (B) and 11.5 dpc trailing ENCDC (B′) gene signature. (C,C′) Violin plots with median values of module scores for wavefront ENCDC genes (C) and trailing ENCDC genes (C′) in ENCDCs at 15.5-18.5 dpc. (D,E) Heatmap representations of differentially expressed genes between clusters associated with (D) wavefront ENCDCs and (E) trailing ENCDCs. Annotations above heatmaps reflect general cell identity annotations. (F) Dot plot visualization of the normalized gene expression and percentage of expressing cells in the postnatal mouse colon. Data were acquired from the mouse colon atlas and were originally produced by Drokhlyansky et al. (2020), who identified genes highly expressed in postnatal enteric glia that are associated with trailing ENCDCs and the glial/progenitor population at 15.5-18.5 dpc.

Fig. 5.

Features of the wavefront-trailing axis of ENCDC maturation are conserved in progenitors at 15.5-18.5 dpc. (A) Unsupervised clustering of the Morarach et al. (2021) dataset containing Wnt1-tdT expressing ENCDCs at 15.5 and 18.5 dpc represented as UMAP. (B,B′) UMAP plot with individual cells colored by module scores generated by the AddModuleScore function of Seurat of the 11.5 dpc wavefront ENCDC (B) and 11.5 dpc trailing ENCDC (B′) gene signature. (C,C′) Violin plots with median values of module scores for wavefront ENCDC genes (C) and trailing ENCDC genes (C′) in ENCDCs at 15.5-18.5 dpc. (D,E) Heatmap representations of differentially expressed genes between clusters associated with (D) wavefront ENCDCs and (E) trailing ENCDCs. Annotations above heatmaps reflect general cell identity annotations. (F) Dot plot visualization of the normalized gene expression and percentage of expressing cells in the postnatal mouse colon. Data were acquired from the mouse colon atlas and were originally produced by Drokhlyansky et al. (2020), who identified genes highly expressed in postnatal enteric glia that are associated with trailing ENCDCs and the glial/progenitor population at 15.5-18.5 dpc.

A similar analysis was performed on the postnatal cells of enteric neurospheres derived from Plp1-expressing glia (Fig. S2D). Postnatal progenitors exhibited the highest scores for the wavefront signature and lowest scores for the trailing signature (Fig. S2E,F). Pseudotime trajectories of glial transition to enteric neurons in neurospheres in vitro paralleled the loss of wavefront-associated gene expression and increase in trailing ENCDC gene expression (Fig. S2G-I, Table S4). Together, these data indicate that, over developmental maturation, progenitor ENCDC populations show transcriptomic features more consistent with enteric glia, which are lost in those committed to becoming enteric neurons.

DUSP6 is required for migration of ENCDCs

Despite wavefront ENCDCs exhibiting highly migratory properties, several studies indicate that many ENCDCs derived from the fully colonized gut at later fetal and postnatal stages retain a migratory ability in explant culture, in vitro assays and in vivo (McKeown et al., 2017; Navoly and McCann, 2021; Stavely et al., 2021). The expression of wavefront ENCDC-associated genes that were differentially expressed in ENCDC progenitors in both the Morarach et al. (2021) fetal dataset and Guyer et al. (2023) postnatal neurosphere dataset were found to include Fbln1, Cald1, Crym, Dusp6, Cnn2, Flna, Hmcn1, Wnt5a, Itgb1, Peak1, Fn1, Peg3 and Col3a1 (Fig. 6A, Table S8). Dusp6 encodes the protein dual specificity phosphatase 6 (DUSP6), which is a MAP kinase phosphatase that acts as a negative-feedback regulator by inhibiting ERK signaling. Dusp6 was selected for further analysis from the genes associated with ENCDC progenitors due to its association with regulating GDNF-RET and ERK signaling in other cell types, and the importance of the same pathways in ENCDC colonization (Lu et al., 2009; Natarajan et al., 2002). UMAP plots indicate that Dusp6 is highly expressed within progenitor ENCDC populations both in embryonic gut post-colonization and in enteric neurospheres (Fig. 6B). Expression of DUSP6 protein was confirmed in ENCDCs at the migratory wavefront at 11.5 dpc by immunohistochemistry (Fig. 6C). ENCDCs in the trailing region at 11.5 dpc and in the proximal hindgut at 13.5 dpc have much lower DUSP6 expression (Fig. 6D-F).

Fig. 6.

Dusp6 is expressed by ENCDC progenitors. (A) Venn diagram of wavefront ENCDC-associated genes that are upregulated in ENCDC progenitors from neurospheres derived from the postnatal small bowel produced by Guyer et al. (2023) and ENCDC progenitors at 15.5-18.5 dpc from Morarach et al. (2021). (B) UMAP visualization of Sox10 and Dusp6 expression in ENCDC progenitors from mouse postnatal neurospheres and ENCDCs at 15.5-18.5 dpc. (C) Cross-sections from the gut of Wnt1-tdT mice at 11.5 dpc at the level of the (I) wavefront ENCDCs showing Wnt1-tdT expression, DUSP6 immunostaining, merged images of Wnt1-tdT and DUSP6, and merged images with DAPI. (C) Cross-sections from the gut of Wnt1-tdT mice at 11.5 dpc at the level of the (I) ceca and (II) trailing ENCDCs showing Wnt1-tdT and DAPI, DUSP6 immunostaining, and merged images of Wnt1-tdT and DUSP6 with DAPI. (D) Representative images at 13.5 dpc of cross-sections from the mouse proximal hindgut showing SOX10 and DUSP6 immunostaining. Area outlined in D is shown at higher magnification in E, demonstrating overlapping expression between SOX10+ cells and DUSP6+ cells. Arrowhead indicates a DUSP6+/Wnt1;tdT+ cell; arrows are DUSP6/Wnt1;tfT+ ENCDCs. (F) Quantification of mean intensity value in the particle expressed in gray-level units of SOX10+/DUSP6+ immunoreactive cells. ***P<0.001, ****P<0.0001 (Kruskal-Wallis test with a post-hoc Dunn's test). Scale bars: 100 μm in C and D; 20 μm in E.

Fig. 6.

Dusp6 is expressed by ENCDC progenitors. (A) Venn diagram of wavefront ENCDC-associated genes that are upregulated in ENCDC progenitors from neurospheres derived from the postnatal small bowel produced by Guyer et al. (2023) and ENCDC progenitors at 15.5-18.5 dpc from Morarach et al. (2021). (B) UMAP visualization of Sox10 and Dusp6 expression in ENCDC progenitors from mouse postnatal neurospheres and ENCDCs at 15.5-18.5 dpc. (C) Cross-sections from the gut of Wnt1-tdT mice at 11.5 dpc at the level of the (I) wavefront ENCDCs showing Wnt1-tdT expression, DUSP6 immunostaining, merged images of Wnt1-tdT and DUSP6, and merged images with DAPI. (C) Cross-sections from the gut of Wnt1-tdT mice at 11.5 dpc at the level of the (I) ceca and (II) trailing ENCDCs showing Wnt1-tdT and DAPI, DUSP6 immunostaining, and merged images of Wnt1-tdT and DUSP6 with DAPI. (D) Representative images at 13.5 dpc of cross-sections from the mouse proximal hindgut showing SOX10 and DUSP6 immunostaining. Area outlined in D is shown at higher magnification in E, demonstrating overlapping expression between SOX10+ cells and DUSP6+ cells. Arrowhead indicates a DUSP6+/Wnt1;tdT+ cell; arrows are DUSP6/Wnt1;tfT+ ENCDCs. (F) Quantification of mean intensity value in the particle expressed in gray-level units of SOX10+/DUSP6+ immunoreactive cells. ***P<0.001, ****P<0.0001 (Kruskal-Wallis test with a post-hoc Dunn's test). Scale bars: 100 μm in C and D; 20 μm in E.

To evaluate the effect of DUSP6 signaling on hindgut colonization, 11.5 dpc guts from Wnt1-tdT mice were cultured ex vivo in the presence of the specific DUSP6 inhibitor BCI (Li et al., 2012; Molina et al., 2009) or vehicle control (Fig. 7A). After 48 h, BCI inhibited ENCDC colonization of the distal hindgut (Fig. 7B). Wavefront ENCDCs were highly disorganized and failed to form a migratory chain pattern along the Tuj1+ neurites of trailing ENCDCs (Fig. 7C). Nevertheless, trailing ENCDCs or Tuj1-expressing neurons did not exhibit any obvious gross morphological changes in response to BCI (Fig. 7D). The functional role of DUSP6 signaling was further explored on avian hindgut ENS development (Fig. 7E). Intestines from the umbilicus to cloaca of E6 chick embryos were cultured in the presence or absence of BCI (Fig. 7F). Intestines cultured in a DMSO vehicle control demonstrated complete ENS colonization, as shown by SOX10 immunofluorescence. Addition of the DUSP6 inhibitor BCI to the culture media arrested ENCDC migration and led to nearly complete hindgut aganglionosis (Fig. 7F, P<0.0001, n=9 guts/group). DUSP6 expression is induced by GDNF and serves as a negative-feedback regulator of ERK signaling in other cell types (Lu et al., 2009). As DUSP6 is known to suppress ERK activity, avian hindguts were cultured in the presence of both BCI and the ERK pathway inhibitor (U0126) to determine whether ERK overactivation was responsible for the arrested ENCDC migration in BCI-treated cultures. The addition of BCI and the ERK inhibitor U0126 restored hindgut colonization by SOX10+ cells (Fig. 7F,G, P<0.01, n=9 guts/group). Previous studies demonstrate that DUSP6 can have roles in regulating proliferation and direct effects on the migratory ability of several cancer cells, including glioblastoma (Kato et al., 2020; Lee et al., 2012; Song et al., 2015; Zuchegna et al., 2020). Samples were labeled for the mitogenic marker phosphohistone H3 to assess cell proliferation (Fig. 7H). BCI significantly reduced ENCDC proliferation in the avian hindgut (P<0.001). Interestingly, however, ERK pathway inhibition restored cell migration but did not ameliorate the reduction in ENCDC proliferation (Fig. 7I, n=7 guts/group). This suggests that a mild reduction in proliferation of ENCDCs does not lead to an arrest in colonization and that DUSP6 could have a direct role in maintaining ENCDC motility via suppression of the ERK pathway.

Fig. 7.

DUSP6 is necessary for ENCDC colonization of the hindgut via suppression of ERK. (A) Cultured mouse gut explants excised from Wnt1-tdT-expressing 11.5 dpc embryos at 48 h post-culture. Explants were treated with a DMSO vehicle as a control or the DUSP6 inhibitor BCI. Scale bars: 1 mm. Arrows indicate the direction of migration. (B) Quantification of the distance of the most distally migrated ENCDCs from the ileocecal junction at 48 h post-culture in the presence of BCI or its DMSO vehicle. Mann–Whitney test, *P<0.05. n=5 gut explants/group. (C,D) Tuj1 immunohistochemistry in 11.5 dpc mouse embryonic gut explants 48 h post-culture visualizing neurons and fiber projections in the ENCDC wavefront (C) and trailing ENCDC regions (D). Arrows indicate the direction of migration. Scale bars: 50 µm. (E) Schematic diagram of chick hindgut colonization assays. (F) Longitudinal sections from organ cultures of the E6 chick gut immunolabelled for SOX10 and stained with DAPI after being cultured in the presence of DMSO (left), the DUSP inhibitor BCI (middle) and BCI+ERK inhibitor (right). Dotted lines indicate the beginning of the hindgut at the ceca. NoR, nerve of Remak; ep, epithelium. (G) Quantification of the length of SOX10+ ENCDC colonization of the hindgut from the ceca. Kruskal–Wallis test with a post-hoc Dunn's multiple comparison test performed. **P<0.01, ****P<0.0001, NS, not significant. n=9 gut explants/group. (H) Longitudinal sections of the E6 chick gut immunolabelled for SOX10 and with anti-phosphohistone H3. mp, myenteric plexus; smp, submucosal plexus. (I) Quantification of the percentage of SOX10+ ENCDCs expressing the phosphohistone H3. Statistical analysis was performed by ANOVA with a post-hoc Tukey test (R Core Team). The P-value was adjusted with Bonferroni's correction. P<0.05 was considered significant and the following further levels of significance were found: **P<0.01; ***P<0.001; ****P<0.0001. Data are mean±s.e.m. n=7 gut explants/group. Scale bars: 100 μm in F and H; 20 μm in H (inset).

Fig. 7.

DUSP6 is necessary for ENCDC colonization of the hindgut via suppression of ERK. (A) Cultured mouse gut explants excised from Wnt1-tdT-expressing 11.5 dpc embryos at 48 h post-culture. Explants were treated with a DMSO vehicle as a control or the DUSP6 inhibitor BCI. Scale bars: 1 mm. Arrows indicate the direction of migration. (B) Quantification of the distance of the most distally migrated ENCDCs from the ileocecal junction at 48 h post-culture in the presence of BCI or its DMSO vehicle. Mann–Whitney test, *P<0.05. n=5 gut explants/group. (C,D) Tuj1 immunohistochemistry in 11.5 dpc mouse embryonic gut explants 48 h post-culture visualizing neurons and fiber projections in the ENCDC wavefront (C) and trailing ENCDC regions (D). Arrows indicate the direction of migration. Scale bars: 50 µm. (E) Schematic diagram of chick hindgut colonization assays. (F) Longitudinal sections from organ cultures of the E6 chick gut immunolabelled for SOX10 and stained with DAPI after being cultured in the presence of DMSO (left), the DUSP inhibitor BCI (middle) and BCI+ERK inhibitor (right). Dotted lines indicate the beginning of the hindgut at the ceca. NoR, nerve of Remak; ep, epithelium. (G) Quantification of the length of SOX10+ ENCDC colonization of the hindgut from the ceca. Kruskal–Wallis test with a post-hoc Dunn's multiple comparison test performed. **P<0.01, ****P<0.0001, NS, not significant. n=9 gut explants/group. (H) Longitudinal sections of the E6 chick gut immunolabelled for SOX10 and with anti-phosphohistone H3. mp, myenteric plexus; smp, submucosal plexus. (I) Quantification of the percentage of SOX10+ ENCDCs expressing the phosphohistone H3. Statistical analysis was performed by ANOVA with a post-hoc Tukey test (R Core Team). The P-value was adjusted with Bonferroni's correction. P<0.05 was considered significant and the following further levels of significance were found: **P<0.01; ***P<0.001; ****P<0.0001. Data are mean±s.e.m. n=7 gut explants/group. Scale bars: 100 μm in F and H; 20 μm in H (inset).

To evaluate whether DUSP6 expression maintains the ability of ENCDCs to migrate post-colonization, the effects of BCI were examined in developmentally mature ENCDCs. Avian midguts at E6 were cultured as explants with GDNF added to induce cell migration on a fibronectin-coated surface. Midgut ENCDCs exhibit robust migration from midgut explants after 24 h of culture in response to GDNF, as previously reported (Nagy et al., 2021) (Fig. 8A, Fig. S3, Movie 1). Similar to results in the hindgut, ENCDC migration was significantly inhibited by BCI (Fig. 8B, P<0.001, n=15 guts/group), and this was partially reversed by addition of the ERK antagonist (Fig. 8C,D, P<0.01, n=15 guts/group). Live cell recordings show that, upon application of BCI, ENCDCs recede back to the explanted midgut and their migratory ability can be restored after removal of BCI from culture (Fig. S3, Movie 2). Results were further replicated in ENCDCs from postnatal neurospheres, with BCI significantly inhibiting the migration of ENCDCs on fibronectin (Fig. 8E), while non-ENCDCs were still able to migrate (Fig. 8E, arrows). In the presence of BCI, inhibition of ERK partially restored the migratory capability of ENCDCs, further validating that DUSP6 functions as a suppressor of ERK to promote cell motility in ENCDC progenitors (Fig. 8E,F).

Fig. 8.

DUSP6-mediated ERK suppression promotes cell migration in post-colonized and postnatal ENCDC progenitors. (A-C) Chick ENCDCs immunolabeled for HNK1 from midgut explants migrating on fibronectin in response to GDNF (A) and after exposure to BCI (B) or BCI and the ERK inhibitor U0126 (C). (D) Quantification of the distance of ENCDC migration from midgut explants after 24 h in culture in the presence of GDNF, GDNF+BCI, and GDNF+BCI and U0126. Statistical analysis was performed using Kruskal–Wallis test with a post-hoc Dunn's test (R Core Team). P<0.05 was considered significant and the following further levels of significance were found: **P<0.01; ***P<0.001; ****P<0.0001. Data are mean±s.e.m. n=15 gut explants/group. (E) Postnatal neurospheres from Wnt1-tdT mice cultured on fibronectin for 48 h. Neurospheres were treated with DMSO (vehicle), BCI or BCI in combination with U0126. Arrows indicate the migration of non-NCC derived cells. Scale bars: 500 µm. (F) Quantification of the area (mm2) of spread of postnatal ENCDCs from neurospheres after 48 h in culture in the presence of DMSO (vehicle), BCI or BCI and U0126. **P<0.01, ****P<0.0001. n=15 images/group.

Fig. 8.

DUSP6-mediated ERK suppression promotes cell migration in post-colonized and postnatal ENCDC progenitors. (A-C) Chick ENCDCs immunolabeled for HNK1 from midgut explants migrating on fibronectin in response to GDNF (A) and after exposure to BCI (B) or BCI and the ERK inhibitor U0126 (C). (D) Quantification of the distance of ENCDC migration from midgut explants after 24 h in culture in the presence of GDNF, GDNF+BCI, and GDNF+BCI and U0126. Statistical analysis was performed using Kruskal–Wallis test with a post-hoc Dunn's test (R Core Team). P<0.05 was considered significant and the following further levels of significance were found: **P<0.01; ***P<0.001; ****P<0.0001. Data are mean±s.e.m. n=15 gut explants/group. (E) Postnatal neurospheres from Wnt1-tdT mice cultured on fibronectin for 48 h. Neurospheres were treated with DMSO (vehicle), BCI or BCI in combination with U0126. Arrows indicate the migration of non-NCC derived cells. Scale bars: 500 µm. (F) Quantification of the area (mm2) of spread of postnatal ENCDCs from neurospheres after 48 h in culture in the presence of DMSO (vehicle), BCI or BCI and U0126. **P<0.01, ****P<0.0001. n=15 images/group.

In this study, we find that migrating wavefront ENCDCs are transcriptionally distinct from the relatively static trailing ENCDCs located just proximally, and we define the transcriptomic signature of each. Analysis of differentially expressed genes shows that ENCDCs are associated with a pro-migratory function characterized by synthesis of ECM molecules, maintenance of mesenchymal-like properties, high proliferative potential and the suppression of cell differentiation. These features distinguish wavefront ENCDCs from their trailing counterparts and collectively contribute to the underlying biology of gut colonization during enteric nervous system development.

ENCDC invasion of the intestine requires synergistic promigratory processes that include trophic support, chemotaxis, proliferation and a permissive microenvironment. In comparison to trailing ENCDCs, the gene signature of the migrating wavefront is enriched for cell proliferation processes (G2M checkpoints and Ef2 targets), confirming that wavefront ENCDCs are more proliferative than trailing ENCDCs, as has been observed in avian species (Simpson et al., 2007). Although the proliferation of wavefront ENCDCs is essential for gut colonization, several other mechanisms also contribute to the continuing migration of ENCDCs. Previously, we identified that ENCDCs can manufacture their own ECM, which impacts whether they continue or halt migration (Akbareian et al., 2013; Nagy et al., 2018). Notably, the ECM molecules identified in these previous studies, including tenascin-C, collagen-18 and agrin, were all upregulated in the wavefront at the gene expression level. Similarly, cell surface receptors and adhesion molecules interacting with the ECM are crucial for mediating the migratory properties of ENCDCs, and several differences in these factors have been identified between wavefront and trailing ENCDCs. For example, Itgb1 was upregulated in wavefront ENCDCs and has previously been shown to be crucial for gut colonization (Breau et al., 2006; Watanabe et al., 2013). We also observed increased expression of one of its ligands, Fn1, at the wavefront.

Our findings further indicate that wavefront ENCDCs possess a mesenchymal-like phenotype in comparison with trailing ENCDCs, and this includes differential expression of several ECM-interacting factors. During delamination of NCCs from the neural tube, these cells undergo EMT before entering the foregut. Our data suggest that, during gut colonization, remnants of this mesenchymal-like signature persist among the immature wavefront ENCDCs and include the expression of cadherin 11, fibrillin 2 and syndecan 4, all of which decline at later stages. Recently, it was reported that fetal ENCDCs have greater migratory potential then postnatal ENCDC progenitors, resulting in better therapeutic value for regenerative therapy (Tian et al., 2021). Therefore, defining and leveraging the properties of developing ENCDCs in postnatal progenitor cells should be explored as a method to improve cell therapy applications.

Datasets pertaining to ENCDC developmental processes were used to correlate wavefront and trailing ENCDC genes with cell maturation and differentiation. We found that wavefront-associated genes are associated with either earlier stages of ENCDC development or with progenitor cell populations. Inversely, trailing ENCDC signatures associate with later developmental stages or neuronal differentiation. In both the developing gut and in enteric neurospheres derived from postnatal intestine, the molecular distinction between neural precursors and enteric glia is ill-defined and these populations appear indistinguishable at the transcriptional level (Guyer et al., 2023; Joseph et al., 2011). The discerning features between progenitors and glia are also unclear postnatally, as lineage-tracing experiments driven by Plp1, Sox2 and Sox10 promoters show that glia can differentiate into enteric neurons, a phenomenon accentuated after injury (Belkind-Gerson et al., 2017; Gershon, 2011; Laranjeira et al., 2011). Nevertheless, our analysis indicates that wavefront ENCDCs are negatively associated with both neural and glial signatures, and therefore appear to exist in an immature state. Alternatively, progenitors derived from the postnatal intestine were negatively associated with the enteric neuron signature and positively associated with the glial signature, which is consistent with neuronal progenitors arising from postnatal glial cells (Belkind-Gerson et al., 2017; Joseph et al., 2011; Laranjeira et al., 2011).

The spatial expression of the neuronal marker Tau and glial marker Plp1, which are upregulated in trailing ENCDCs, were validated using transgenic GFP reporter mice. These data confirmed that TauGFP expression occurs predominately in trailing ENCDCs. Plp1GFP, however, is expressed at low levels in the wavefront, with levels rising along the wavefront-trailing axis of the embryonic gut, consistent with the gene expression data. Given that Plp1 expression is restricted to glia in the postnatal environment, these data appear to be indicative of a gradual cell maturation to a more glia-like phenotype in trailing ENCDCs. This may not necessarily indicate a restricted glial cell fate, however, given their ability to become enteric neurons postnatally. Previously, we demonstrated that premature differentiation of ENCDCs leads to intestinal aganglionosis (Nagy et al., 2016). Likewise, enteric neurons lose their migratory potential after neuronal differentiation (Hao et al., 2009). Therefore, the unique ability of wavefront ENCDCs to maintain an immature state throughout gut colonization is essential for successful migration and ENS formation.

We identified Dusp6 in our wavefront signature, and this gene also exhibited high expression in progenitor cell populations at later stages of development and also in postnatal-derived enteric neurospheres. DUSP6 is a MAP kinase phosphatase that acts as a negative-feedback regulator by inhibiting ERK signaling. Notably, Dusp6 expression is inducible by RET-GDNF signaling (Lu et al., 2009), a key pathway in enteric nervous system development, and ERK signaling downstream of GDNF stimulation is essential for ENCDC migration (Natarajan et al., 2002). Paradoxically, overstimulation of ENCDCs by GDNF overexpression induces premature differentiation and aganglionosis (Mwizerwa et al., 2011; Nagy et al., 2020), indicating that GDNF signaling is finely tuned for temporal and cell-specific effects on ENCDCs. In ex vivo cultures of 11.5 dpc mouse intestine, inhibition of DUSP6 by BCI significantly suppressed ENCDC migration with no obvious effects on trailing ENCDCs. Similar results were observed in avian hindgut, with BCI arresting ENCDC hindgut colonization. We further validated that inhibition of migration is dependent on ERK overactivation by reversing these effects with the ERK inhibitor, U0126. BCI was also shown to reduce cell proliferation, although addition of the ERK antagonist rescued migration defects without promoting cell proliferation, suggesting that DUSP6-ERK signaling may have a direct effect on ENCDC migration. Likewise, BCI inhibited migration of ENCDCs from a post-colonization environment, including explanted E6 avian midgut and ENCDCs in postnatal enteric neurospheres. ERK antagonism rescued migration in both models, further validating these findings. Notably, the few non-neural crest-derived cells within neurosphere cultures still demonstrated an ability to migrate in the presence of BCI, suggesting the effects are specific to ENCDCs. These data indicate that DUSP6 is an important modulator of intracellular signaling in immature ENCDCs that maintains their promigratory properties by preventing overactivation of ERK.

The present study has several limitations worth considering. Multiple analyses were performed to study the expression patterns of 11.5 dpc wavefront ENCDC-associated genes in other single-cell RNA-seq datasets. Only one of these datasets contains ENCDCs at 11.5 dpc (Cao et al., 2019) and these are present in a very limited number with even fewer likely to be wavefront ENCDCs. Therefore, these analyses may not reflect the complete spatiotemporal dynamics in gene expression that could underlie the transcriptional pattern of wavefront ENCDCs and their unique properties. A study that similarly performed transcriptomics in a spatial nature was performed on migrating chick cranial NCCs, and a cross-species comparison was performed (Morrison et al., 2017). Although similarities between the datasets are observed, a lack of conservation in genes, or their functions, between species needs to be taken into consideration when interpreting these data. Similarly, gene sets of human enteric neuron and glial markers were used to assess the overall signature of these cells types in the wavefront and trailing ENCDCs, which would not capture the genes that are distinct between these species. Nevertheless, there appears to be enough overlap of markers of enteric neurons and glia to justify this approach (Drokhlyansky et al., 2020; May-Zhang et al., 2021). Data variability between individual experimental runs of wavefront and trailing ENCDCs was also observed, which could have been caused by slight differences in embryo staging, sex, the cell isolation process or the number of cells obtained. This variation could lead to reduced power to observe changes in the expression of other genes that could be essential in ENCDC migration. Single-cell RNA-seq or similar technologies could be useful to overcome some of these limitations and may also assist in determining whether wavefront ENCDCs are a uniform or heterogeneous cell population. However, this technology is currently limited by being capable of detecting only a limited number of features compared with bulk RNA-seq, loss of low abundant transcripts, variability in cell complexity and quality, and crucially relies on an ability to provide an adequate number of cells for meaningful analysis. The latter is a major limitation for ENCDC wavefront studies given their very small numbers.

The current study was performed at 11.5 dpc, which provides enough distance between wavefront and trailing ENCDCs to reliably dissect these regions for cell isolation. Although wavefront ENCDCs at 11.5 dpc migrate considerably faster than trailing ENCDCs, this population begins to slow down from 45 μm/h at 10.5 dpc in the ileum to 32 μm/h as they approach the cecum (Druckenbrod and Epstein, 2005; Young et al., 2014). It is worth noting that transmesenteric ENCDCs (Nishiyama et al., 2012), Schwann cell precursors (Uesaka et al., 2015) and sacral NCCs (Burns et al., 2000) also contribute to the formation of the distal two-thirds of the ENS in the hindgut, whereas ENCDCs migrating through the cecum primarily form the ENS in the proximal hindgut in mice. Future studies exploring the genes driving the migration patterns in these different cell populations will be of great importance to fully comprehend the development of the ENS.

Mice and breeding

All mice were purchased from Jackson Laboratory with the exception of NestinGFP and Plp1GFP mice which were kindly donated by Jorg Dietrich MD PhD (Massachusetts General Hospital, Boston, MA, USA) (Mignone et al., 2004) and Wendy Macklin PhD (University of Colorado, Bouler, USA) (Mallon et al., 2002) (Table 1). All mice were housed and bred under specific pathogen-free conditions at the Center for Comparative Medicine animal facility at Massachusetts General Hospital (MGH). Rodents were housed in Allentown rectangular caging (160 cages per individually ventilated cage racks; which uses blower at 60 air changes/h) under a 12 h:12 h light:dark cycle from 07:00 am to 19:00. Bedding consisted of Hardwood Sanichip; with Carefresh nesting material and mice had access to Prolab Isopro RMH 3000 chow mix (ScottPharma) ad libitum. Expression of Plp1 in ENCDCs was visualized in double transgenic mice by crossing neural crest-specific Wnt1::Cre mice with glia-specific Plp1GFP mice to generate Plp1GFP;Wnt1::Cre mice. These mice were then crossed with ROSA26tdTomato reporter mice to generate Plp1GFP;Wnt1::Cre;ROSA26tdTomato (Plp1GFP;Wnt1-tdT) mice. ENCDC differentiation into neurons was visualized using the same strategy with TauGFP mice substituted for Plp1GFP mice to generate TauGFP;Wnt1::Cre;ROSA26tdTomato (TauGFP;Wnt1-tdT) mice. Likewise, ENCDC expression of nestin was performed in NestinGFP;Wnt1::Cre;ROSA26tdTomato (NestinGFP;Wnt1-tdT) mice. Experiments were approved by the MGH Institutional Animal Care and Use Committee.

Table 1.

Mice strains used for breeding in this study

Mice strains used for breeding in this study
Mice strains used for breeding in this study

Avian embryo

Fertilized White Leghorn chicken (Gallus gallus domesticus) eggs were obtained from commercial breeders (Prophyl-BIOVO, Hungary) and maintained at 37.5°C in a humidified HEKA 1+ egg incubator (Rietberg, Germany). Gut development stages were referenced to the chick embryo gut staging table (Southwell, 2006) and the ENS formation timetable (Nagy et al., 2012).

Embryonic cell isolation

Wnt1-Cre;tdTomato mouse embryos, in which all ENCDCs are fluorescently labelled, were acquired from pregnant mice at 11.5 days post fertilization. Guts were dissected from 11.5 dpc mouse embryos to isolate two cell populations: (1) wavefront cells up to 0.5 mm proximal to the wavefront; and (2) trailing ENCDCs 1.0-1.5 mm proximal to the wavefront. After dissecting these segments, tissues were enzymatically dissociated in collagenase-dispase [10 min, 37°C, dispase (250 μg/ml; StemCell Technologies) and collagenase XI (1 mg/ml; Sigma-Aldrich)]. The solution was resuspended in PBS and filtered through a 40 µm cell strainer (Corning) to obtain primary single cell suspensions for purification of ENCDCs by flow cytometry on a BD FACSAria II Cell Sorter (Becton Dickinson). Experiments were performed in duplicate, and seven embryos were used per cell-sorting experiment. All experiments were conducted at the same time of day (mornings before noon) for consistency. Cells were immediately sorted in Trizol reagent for RNA extraction processing.

RNA-seq and data processing

RNA was extracted using Trizol reagent as of manufactures instructions. Samples were submitted to the MGH Next Generation Sequencing Core facility where RNA integrity was validated using a 2100 Bioanalyzer (Agilent Technologies) microfluidics platform with the RNA 6000 Nano Kit (Agilent Technologies) and mRNA libraries were constructed using the Illumina Nextera XT kit according to the manufacturer's protocols. High-throughput sequencing was performed using a 50 bp paired-end read protocol on the Illumina MiSeq System. The dataset preprocessing was performed using the Galaxy project, open source, web-based platform for biomedical research (Afgan et al., 2018). The four samples raw reads were aligned using the STAR alignment tool version 2.5.2b (Dobin et al., 2013) and the Genome Reference Consortium Mouse Build 38 (GRCm38/mm10). Alignment quality control was performed with FastQC Version 0.11.9 (Andrews, 2010) and one sample was discarded due to low quality of reads. Expression levels were quantified using featureCounts (Liao et al., 2014). Principal component analysis (PCA) was conducted and plots generated using the default settings (removeVar=0.1) of the R package PCAtools (Blighe et al., 2019). The normalized counts and differential expression analysis was performed on Galaxy using the DeSeq2 package (Love et al., 2014). Genes with a P<0.05 and LogFC±0.5 were considered to be differentially expressed genes (DEGs). Genes upregulated in this analysis were considered the wavefront ENCDC-associated gene signature and genes downregulated were considered the trailing ENCDC-associated signature for downstream analysis.

PPI and over-representation analysis

Protein-protein interaction networks were identified from the STRING functional protein association networks database (evidence score >400) using NetworkAnalyst 3.0. Modules of interacting genes were identified by the label propagation method and annotated for over-representation of Gene Ontologies (Biological Processes) using DAVID (Dennis et al., 2003).

Gene set enrichment analysis (GSEA)

Gene set enrichment analysis (GSEA) was performed with GSEA v4.1.0 [build 27] (Broad Institute) software on gene lists pre-ranked by sign(log2FC)×log10(P-value). Genesets were obtained from MSigDB, this included HALLMARK gene sets of biological states or processes and enteric neuron and glia signatures descartes_fetal_intestine_ens_glia and descartes_fetal_intestine_ens_neuron, which are derived from Cao et al. (2020). In addition, a manually curated gene set was generated from the single-cell RNA-seq data of wavefront and trailing NCCs in the cranial neural crest cell migratory streams from embryonic chicks entering the paraxial mesoderm and branchial arches reported by Morrison et al. (2017). To generate these gene sets, the upregulated and downregulated genes between ‘front’ and ‘trail’ NCCs at HH13 and HH15 provided in Morrison et al. (2017) were taken at a cutoff of P<0.05 and Log2FC >2 (front gene set) or Log2FC <2 (trail gene set). Heatmaps for genes of interest identified in the aforementioned analysis were visualized using the web-based tool Morpheus and are presented as z-score distributions (across samples) and clustered by Euclidean distance with average linkage.

Public enteric neural crest-derived cell dataset processing

Several pre-existing datasets were used in this study to inform our analysis of the developing ENCDC transcriptome. Published single-cell RNA-Seq data during organogenesis in the mouse (Cao et al., 2019) was analyzed to determine whether the transcriptomic properties associated with spatial relation to the migrating ENCDCs of the wavefront were correlated with developmental maturity. Processed and filtered data and associated metadata (as used by Cao et al., 2019) were acquired from https://oncoscape.v3.sttrcancer.org/atlas.gs.washington.edu.mouse.rna/downloads and were analyzed using Seurat. Data were refined to focus on the neural crest trajectories of the original authors and were refined by the annotations [Neural crest (PNS neuron) trajectory 3 and Neural crest (PNS glia) trajectory 2]. Dimension reduction was performed by running the NormalizeData, ScaleData, FindVariableFeatures, RunPCA, FindNeighbors(dims=1:15) and RunUMAP(dims=1:15) commands. These values were selected using the visual aid of elbow plots and determining the first principal component where the percentage of total variation was less than 0.1% of the previous principal component, indicating a plateau in meaningful variation. Clusters were identified by FindClusters(resolution=0.8) using the Louvain algorithm. Single-cell RNA-Seq data of Wnt1-tdT cells from the gut of 15.5-18.5 dpc mice were also used (Morarach et al., 2021) (accession number SRP258962). Data were processed as described previously (Guyer et al., 2023). The datasets were filtered with removal of cells<1000 and >6000 features detected, >40,000 total UMI counts, and>5% mitochondrial genes detected. PCA and UMAP were undertaken with the top 30 principal components. Cell clusters were identified using the FindClusters command (resolution=0.5) identifying 13 cell clusters, including one cluster of contaminating non-neural crest-derived cells that was removed according to the analysis of the original authors (Morarach et al., 2021). To confirm clusters from the Cao et al. (2019) dataset exhibited an identity consistent with ENCDCs, data were integrated with the Morarach et al. (2021) dataset. Datasets were made consistent by including the intersect of genes between the two and filtering Feature_RNA>800 & nFeature_RNA<7000 & percent of mitochondrial genes<15. Data integration with anchors was performed in Seurat by running FindIntegrationAnchors(max.features=200, k.filter=50, k.anchor=3) and IntegrateData(dims=1:20), which included anchor features and genes identified in the differential expression analysis of wavefront and trailing ENCDCs at 11.5 dpc included. Dimension reduction and clustering were performed using the command above with dims=20 and resolution 0.2. Clusters to successfully integrate with the Morarach et al. (2021) dataset were considered to be bona fide ENCDCs and included clusters 13 and 9, which corresponded to the original authors Enteric neuron trajectory_1 and enteric glia and Schwann cell trajectory_2 annotations from the Cao et al. (2019) dataset. These and their adjacent NCC progenitor clusters (15 and 1) were taken for subsequent analysis to generate enteric neuron development and enteric glia development data. Dimension reduction was performed on these data as above, with dims=15. Data of Plp1GFP-expressing cells isolated from enteric neurospheres generated from postnatal small bowel were obtained from Guyer et al. (2023) and were processed as in the original publication. Differential gene expression analysis was performed between developmental stages 9.5, 10.5, 11.5 12.5 and 13.5 dpc in the Cao et al. (2019) dataset, using the FindAllMarkers() function of Seurat. Differential expression analysis was performed between clusters in the Morarach et al. (2021) and Guyer et al. (2023) datasets using the same function. These analyses were used to generate a ranked list for GSEA and to determine the differential expression of wavefront and trailing ENCDC-associated genes from the 11.5 dpc data. Data from the mouse colon atlas originally produced from single-nuclei Seq of adult C57Bl/6 mice from (Drokhlyansky et al., 2020) was downloaded from the Broad Institute single cell portal https://singlecell.broadinstitute.org/single_cell/study/SCP1038/. The processed data files for mouse large intestine (10X) data downloaded included the counts matrix (gene_sorted-mli.matrix.mtx), barcodes (mli.barcodes.tsv), genes (mli.genes.tsv) and annotations from the original authors (mli.tsne2.txt).

Wavefront and trailing ENCDC scores

Scoring of wavefront and trailing ENCDC-associated signatures was performed using the AddModuleScore function of Seurat. This function scores the expression of input gene lists by comparing expression to computed sets of equivalent control genes to each cell individually using methods derived from Tirosh et al. (2016). Gene sets were determined by genes with a P<0.05 and a LogFC>0.5 for the wavefront gene set, and LogFC<0.5 for the trailing gene set, after differential expression analysis between wavefront and trailing ENCDCs at 11.5 dpc. The AddModuleScore function was used with default parameters and control features set to 100. Cells were scored within the glia and neuron developmental trajectories generated from Cao et al. (2019) at 9.5, 10.5, 11.5, 12.5 and 13.5 dpc or between clusters for the Morarach et al. (2021) and Guyer et al. (2023) datasets.

Pseudotime analysis

Analysis of developmental pseudotime in single-cell RNA-seq datasets was performed using the R package Monocle3. Cells were clustered using the cluster_cells(reduction_method=‘UMAP’, k=40) function. Trajectories were determined by the learn_graph function and principle nodes were plotted in UMAP space to assist in the determination of root nodes for subsequent analysis. For the Cao et al. (2019) dataset, this included a root node with the majority of cells at 9.5 dpc, or within the progenitor cell cluster in the Guyer et al. (2023) dataset and integrated Cao et al. (2019) and Morarach et al. (2021) datasets. Pseuotime was culculated using the order_cells() function and subsets of the data were generated using choose_graph_segments() and selecting the first and last node in the trajectory. This included nodes consistent with the 9.5-13.5 dpc cells for the Cao et al. (2019) dataset and progenitor to neuron clusters in the Guyer et al. (2023) dataset and integrated Cao et al. (2019) and Morarach et al. (2021) datasets visualized by UMAP. A generalized linear model for gene expression changes in pseudotime was performed using the functions fit_models(model_formula_str=‘∼pseudotime’) and coefficient_table() with a q_value <= 0.05 considered significant. Genes from the 11.5 wavefront and trailing signature that were downregulated (estimate <0) or upregulated (estimate >0) over pseudotime were plotted using the ComplexHeatmap package of R.

Quantification of fluorescence intensity in ENCDCs

GFP expression driven off the Plp1 promoter was quantified as a surrogate measure of spatial Plp1 gene expression in the developing gut. The embryonic mouse gut was dissected and imaged in glass bottom petri dishes on the Keyence BZ-X800 All-in-one microscope. Imaging parameters were consistent for all experiments [for Plp1GFP: filter, BZ-X filter GFP (model OP87763); resolution, 960×720; objective lens, PlanApo_λ 4×0.20/20.00 mm; exposure time, 4 s; gain, +6 dB; mercury quantity, 20%; capturing mode, Monochrome 8bit; binning, 2×2; and for Wnt1-tdT: filter, BZ-X filter TexasRed (model OP87765); exposure time, 4.5 s, with the reminder parameters the same as Plp1GFP]. The mean fluorescence intensity of GFP from Wnt1-tdT;Plp1GFP mice was measured using previously described methods (Stavely et al., 2022) by ImageJ software. Wnt1-tdT positivity was used to define cell bodies and region of interest (ROI). A circular ROI of 5 µm in diameter was drawn in each Wnt1-tdT cell and overlaid onto images of GFP to measure Plp1 expression specifically in neural crest-derived cells. Cells were further classified by their relative distance from the most caudal wavefront ENCDCs. For quantification of DUSP6 immunoreactivity, the ROI of ENCDCs were defined by binary thresholding of Wnt1-tdT expression and mean gray values in corresponding images of DUSP6 immunofluorescence were recorded using ImageJ.

Ex vivo culture of embryonic gut

Cultures of mouse embryonic gut explants were performed similar to previously described methods (Nagy et al., 2020). Guts were dissected from 11.5 dpc Wnt1-tdT mice before being pinned at the resected margins of gut in silicone-lined cell culture dishes. The dissected gut was cultured in standard cell cultures conditions (37°C, 5%CO2) in tissue culture media containing 1% penicillin-streptomycin (Life Technologies) and 10% FBS in alpha-MEM (Life Technologies) with daily media changes. Samples were treated with either BCI hydrochloride (MedChemExpress, HY-115502A) at concentration of 20 µM in DMSO or an equivalent volume of DMSO as a vehicle control. After 48 h in culture, samples were imaged using a Keyence all-in one microscope to assess migration of the wavefront. The distance of the most caudally migrated ENCDC from the ileocecal junction was recorded. Samples were fixed in 4% PFA after 48 h in culture for processing of samples for immunohistochemistry.

Immunohistochemistry

Tissues were processed for immunohistochemistry in whole-mount and cross-section preparations a previously described (Nagy et al., 2018, 2021). Briefly, the mouse gut was fixed in 4% PFA overnight and washed in PBS for 3×10 min. For, whole-mount preparations tissues were permeabilized and blocked in PBS containing the detergent 0.1% Triton X-100 and 10% normal donkey serum for 1 h. Samples were then incubated with the antibodies anti-TUBB3 (Tuj1; a neuron-specific class III β-tubulin) conjugated with either Alexa Fluor 647 (801210, Biolegend, 1:1000) or Alexa Fluor 594 (801208, Biolegend, 1:1000), rabbit anti-S100b (Neomarkers, RB-044-A, 1:100) or human anti-Hu (Anna-1)-positive serum (1:16,000, a gift from Dr Vanda Lennon; Mayo Clinic, Rochester, MN, USA) overnight before incubating with secondary antibodies Alexa Fluor 594 donkey anti-rabbit (A21207; ThermoFisher Scientific; 1:200) or Alexa Fluor 647 goat anti-human (A21445; ThermoFisher Scientific; 1:200) before washing in PBS and mounting with Aqua-Poly/Mount (Polysciences, 18606) on glass slides. For immunocytochemistry in cross-sections, the embryonic tissue samples were fixed overnight in 4% PFA at 4°C, then infiltrated with 15% sucrose in PBS for48 h at 4°C. The sucrose was changed for 7.5% gelatin (Sigma-Aldrich, G9391) containing 15% sucrose and the tissues rapidly frozen at −50°C in 2-methylbutane (Sigma-Aldrich, M32631). Frozen sections were cut at 12 μm, collected on poly-L-lysine-coated slides (Sigma-Aldrich, P8920) and stained by immunocytochemistry, as previously described (Nagy and Goldstein, 2006). Sections were labelled with rabbit anti-DUSP6 (Abcam, ab76310; dilution), anti-fibrillin 2 polyclonal antibody (PA5-52995; ThermoFisher Scientific; 1:50), anti-syndecan-4 antibody (NB110-41551; NovusBiologicals; 1:30), anti-heparan sulfate proteoglycan 2/perlecan antibody (ab2501; Abcam; 1:100), anti-crym polyclonal antibody (PA5-53920; ThermoFisher Scientific; 1:20), anti-vimentin (sc-7557; SantaCruz Biotech; 1:100) and anti-p75 (G3231, 1:200, Promega), which recognizes the cytoplasmic domain of the mouse-specific p75 neurotrophin receptor (P75NTR) on the surface of migratory and postmigratory neural crest-derived cells; anti-HuC/D (clone 16A11, 1:50, ThermoFisher Scientific), which recognizes a neuron-specific RNA-binding protein (ELAV-like protein 4); anti-Tuj1 (TUBB3; clone B1195, 1:400, Covance), anti-Sox10 (specific for ENCDCs and neural precursors 1:30, clone A-2, Santa Cruz Biotechnology, sc-365692), anti-PHOX2B (specific for ENCDCs and glial precursors, 1:500, a gift from Dr Hideki Enomoto, Kobe University, Hyogo, Japan) and goat anti-GFP (green florescent protein, 1:200, Rockland, 600-101-215M). Secondary antibodies included Alexa Fluor 488 donkey anti-rat IgG (A21208; ThermoFisher Scientific; 1:200), Alexa Fluor 488 donkey anti-rabbit IgG (A21206; ThermoFisher Scientific; 1:200), Alexa Fluor 488 goat anti-mouse IgG (A-32723, ThermoFisher Scientific; 1:200), Alexa Fluor 594 goat anti-mouse IgG (A-32742; ThermoFisher Scientific; 1:200) and Alexa Fluor 488 donkey anti-goat IgG (A11055; ThermoFisher Scientific; 1:200).

Postnatal neurosphere migration

Neurospheres were isolated as previously described (Stavely et al., 2021). Four-week-old mice were sacrificed to obtain the small intestine. The smooth muscle and myenteric plexus layer were separated from intestinal tissue and dissociated with dispase (250 µg ml−1; StemCell Technologies) and collagenase XI (1 mg ml−1; Sigma-Aldrich) at 37°C for 1 h with gentle pipetting. The cell suspension was plated into T-25 flask at 50,000 cells ml−1 in of proliferation medium [Neurocult NSC Basal Medium with NeuroCult Proliferation Supplement for mouse and rat (StemCell Technologies)], supplemented with 20 ng ml−1 epidermal growth factor (EGF; StemCell Technologies) and 10 ng ml−1 basic fibroblast growth factor (bFGF; StemCell Technologies), 0.0002% heparin (StemCell Technologies) and 100 U ml−1 penicillin-streptomycin (Life Technologies) for 10 days to form neurospheres. Neurospheres were cultured in fibronectin-coated (1:500 for 2 h at 37°C, ThermoFisher Scientific) 24-well plates to assess cell migration. The area covered (mm2) by Wnt1-tdT expressing cells from each seeded neurosphere was measured using binary threshold analysis (Stavely et al., 2022) with ImageJ software. Neurospheres were treated with the DUSP6 inhibitor BCI (10 µM) with or without the ERK pathway inhibitor U0126 (10 µM, Cayman Chemical Company) or an equal volume of the DMSO vehicle as a control.

Chick intestinal organ culture

Organ cultures of the developing chick gut were performed using the methods previously described (Nagy et al., 2016, 2021). E6 chick intestine from umbilicus to cloaca was pinned to Sylgard and cultured for 48 h in DMEM (Sigma-Aldrich, D5796) containing 100 U/ml penicillin-G, 0.1 mg/ml streptomycin (Gibco, 15140122). Cultures were treated with the DUSP6 inhibitor BCI (20 μM) with or without ERK inhibitor U0126 (20 μM), or with DMSO as a vehicle control. Tissues were fixed and immunohistochemistry was performed using anti-SOX10 (1:30, clone: A-2, SantaCruz Biotech., sc-365692) as previously described (Nagy et al., 2021). The length of SOX10+ ENCDC colonization of the hindgut from the ceca were measured using ImageJ.

Cell migration assay

For ENCDC migration assays, post-umbilical midgut was removed from E6 (HH29) chick embryos and cultured on 20 μg/ml fibronectin-coated dishes treated with GDNF (10 ng/ml; R&D Systems, 212-GD-010; n=9) with the DUSP6 inhibitor BCI (10 μM; n=9) with or without ERK inhibitor U0126 (10 μM; n=9). Cultures were repeated with removal of GDNF after the first 12 h, followed by 8 h in the presence of BCI and incubation for another 12 h with GDNF without BCI (n=9) or addition of GDNF for 32 h (n=9). For live cell migration imaging, an Olympus Provi CM20H Incubation Monitoring System was used.

Statistical analysis

All data are presented as mean±s.e.m. unless indicated otherwise. Data were analyzed using Graphpad Prism version 9.2 software. For pairwise comparisons, a Mann–Whitney test was performed. Data with multiple groups were assessed using a Kruskal–Wallis ANOVA test with a post-hoc Dunn's multiple comparison test. For all analysis, P<0.05 was considered statistically significant unless stated otherwise.

The authors are grateful for expertise offered by the Harvard University Bauer Core and the Harvard Stem Cell Institute Center for Regenerative Medicine Flow Cytometry Core facility.

Author contributions

Conceptualization: R.S., A.M.G., N.N.; Methodology: R.S., R.H., N.P., R.A.G., A.S., E.S., N.N.; Software: R.S., R.A.G.; Validation: R.S., R.H., A.S., E.S., N.N.; Formal analysis: R.S., R.H., A.A.R., A.S.; Investigation: R.S., N.P., M.O., A.S., E.S., J.M., A.M.G., N.N.; Resources: R.S., N.N.; Data curation: R.S.; Writing - original draft: R.S., A.M.G., N.N.; Writing - review & editing: R.S., R.G., A.M.G., N.N.; Visualization: R.S., N.N.; Supervision: A.M.G., N.N.; Project administration: R.S., N.N.; Funding acquisition: R.S., A.M.G., N.N.

Funding

This work was supported by the National Institutes of Health (R01DK103785 and R21HD106036 to A.M.G., R21HD106036 to R.S., and F32DK121440 to R.A.G.); by the Research Excellence Programme of the Innovációs és Technológiai Minisztérium in Hungary, within the framework of the TKP2021-EGA-25 thematic programme of Semmelweis University [to N.N.]; and by a Hungarian Science Foundation NKFI grant [138664 to N.N.]. R.S. is supported by a Massachusetts General Hospital Future Medical Discovery award. R.A.G. is supported by the Massachusetts General Hospital Patricia K. Donahoe Resident Research Catalyst Award. Deposited in PMC for release after 12 months.

Data availability

RNA-seq data have been deposited in the Gene Expression Omnibus database under accession number GSE217757.

Afgan
,
E.
,
Baker
,
D.
,
Batut
,
B.
,
Van Den Beek
,
M.
,
Bouvier
,
D.
,
Cech
,
M.
,
Chilton
,
J.
,
Clements
,
D.
,
Coraor
,
N.
,
Grüning
,
B. A.
et al. 
(
2018
).
The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update
.
Nucleic Acids Res.
46
,
W537
-
W544
.
Akbareian
,
S. E.
,
Nagy
,
N.
,
Steiger
,
C. E.
,
Mably
,
J. D.
,
Miller
,
S. A.
,
Hotta
,
R.
,
Molnar
,
D.
and
Goldstein
,
A. M.
(
2013
).
Enteric neural crest-derived cells promote their migration by modifying their microenvironment through tenascin-C production
.
Dev. Biol.
382
,
446
-
456
.
Andrews
,
S.
(
2010
).
FastQC: A Quality Control Tool for High Throughput Sequence Data
.
Cambridge
,
United Kingdom
:
Babraham Bioinformatics, Babraham Institute
.
Belkind-Gerson
,
J.
,
Graham
,
H. K.
,
Reynolds
,
J.
,
Hotta
,
R.
,
Nagy
,
N.
,
Cheng
,
L.
,
Kamionek
,
M.
,
Shi
,
H. N.
,
Aherne
,
C. M.
and
Goldstein
,
A. M.
(
2017
).
Colitis promotes neuronal differentiation of Sox2+ and PLP1+ enteric cells
.
Sci. Rep.
7
,
2525
.
Bhave
,
S.
,
Arciero
,
E.
,
Baker
,
C.
,
Ho
,
W. L.
,
Stavely
,
R.
,
Goldstein
,
A. M.
and
Hotta
,
R.
(
2019
).
Enteric neuronal cell therapy reverses architectural changes in a novel diphtheria toxin-mediated model of colonic aganglionosis
.
Sci. Rep.
9
,
18756
.
Blighe
,
K.
,
Lewis
,
M.
,
Lun
,
A.
and
Blighe
,
M. K.
(
2019
).
Package ‘PCAtools.’
Breau
,
M. A.
,
Pietri
,
T.
,
Eder
,
O.
,
Blanche
,
M.
,
Brakebusch
,
C.
,
Fässler
,
R.
,
Thiery
,
J. P.
and
Dufour
,
S.
(
2006
).
Lack of β1 integrins in enteric neural crest cells leads to a Hirschsprung-like phenotype
.
Development
133
,
1725
-
1734
.
Burns
,
A.
,
Champeval
,
D.
and
Le Douarin
,
N.
(
2000
).
Sacral neural crest cells colonise aganglionic hindgut in vivo but fail to compensate for lack of enteric ganglia
.
Dev. Biol.
219
,
30
-
43
.
Cao
,
J.
,
Spielmann
,
M.
,
Qiu
,
X.
,
Huang
,
X.
,
Ibrahim
,
D. M.
,
Hill
,
A. J.
,
Zhang
,
F.
,
Mundlos
,
S.
,
Christiansen
,
L.
,
Steemers
,
F. J.
et al. 
(
2019
).
The single-cell transcriptional landscape of mammalian organogenesis
.
Nature
566
,
496
-
502
.
Cao
,
J.
,
O'day
,
D. R.
,
Pliner
,
H. A.
,
Kingsley
,
P. D.
,
Deng
,
M.
,
Daza
,
R. M.
,
Zager
,
M. A.
,
Aldinger
,
K. A.
,
Blecher-Gonen
,
R.
,
Zhang
,
F.
et al. 
(
2020
).
A human cell atlas of fetal gene expression
.
Science
370
,
eaba7721
.
Dennis
,
G.
,
Sherman
,
B. T.
,
Hosack
,
D. A.
,
Yang
,
J.
,
Gao
,
W.
,
Lane
,
H. C.
and
Lempicki
,
R. A.
(
2003
).
DAVID: database for annotation, visualization, and integrated discovery
.
Genome Biol.
4
,
P3
.
Dobin
,
A.
,
Davis
,
C. A.
,
Schlesinger
,
F.
,
Drenkow
,
J.
,
Zaleski
,
C.
,
Jha
,
S.
,
Batut
,
P.
,
Chaisson
,
M.
and
Gingeras
,
T. R.
(
2013
).
STAR: ultrafast universal RNA-seq aligner
.
Bioinformatics
29
,
15
-
21
.
Drokhlyansky
,
E.
,
Smillie
,
C. S.
,
Van Wittenberghe
,
N.
,
Ericsson
,
M.
,
Griffin
,
G. K.
,
Eraslan
,
G.
,
Dionne
,
D.
,
Cuoco
,
M. S.
,
Goder-Reiser
,
M. N.
,
Sharova
,
T.
et al. 
(
2020
).
The human and mouse enteric nervous system at single-cell resolution
.
Cell
182
,
1606
-
1622.e23
.
Druckenbrod
,
N. R.
and
Epstein
,
M. L.
(
2005
).
The pattern of neural crest advance in the cecum and colon
.
Dev. Biol.
287
,
125
-
133
.
Druckenbrod
,
N. R.
and
Epstein
,
M. L.
(
2007
).
Behavior of enteric neural crest-derived cells varies with respect to the migratory wavefront
.
Dev. Dyn.
236
,
84
-
92
.
Druckenbrod
,
N. R.
and
Epstein
,
M. L.
(
2009
).
Age-dependent changes in the gut environment restrict the invasion of the hindgut by enteric neural progenitors
.
Development
136
,
3195
-
3203
.
Gershon
,
M. D.
(
2011
).
Behind an enteric neuron there may lie a glial cell
.
J. Clin. Invest.
121
,
3386
-
3389
.
Graham
,
H. K.
,
Maina
,
I.
,
Goldstein
,
A. M.
and
Nagy
,
N.
(
2017
).
Intestinal smooth muscle is required for patterning the enteric nervous system
.
J. Anat.
230
,
567
-
574
.
Guyer
,
R. A.
,
Stavely
,
R.
,
Robertson
,
K.
,
Bhave
,
S.
,
Mueller
,
J.
,
Picard
,
N.
,
Hotta
,
R.
,
Kaltschmidt
,
J. A.
and
Goldstein
,
A. M.
(
2023
).
Single-cell multiome sequencing clarifies enteric glial cell diversity and identifies an intraganglionic population poised for neurogenesis
.
Cell Rep
.
(in press)
.
Hao
,
M. M.
,
Anderson
,
R. B.
,
Kobayashi
,
K.
,
Whitington
,
P. M.
and
Young
,
H. M.
(
2009
).
The migratory behavior of immature enteric neurons
.
Dev. Neurobiol.
69
,
22
-
35
.
Heanue
,
T. A.
and
Pachnis
,
V.
(
2007
).
Enteric nervous system development and Hirschsprung's disease: advances in genetic and stem cell studies
.
Nat. Rev. Neurosci.
8
,
466
-
479
.
Hotta
,
R.
,
Anderson
,
R. B.
,
Kobayashi
,
K.
,
Newgreen
,
D. F.
and
Young
,
H. M.
(
2010
).
Effects of tissue age, presence of neurones and endothelin-3 on the ability of enteric neurone precursors to colonize recipient gut: implications for cell-based therapies
.
Neurogastroentero. Motil.
22
,
331
-
386
.
Hotta
,
R.
,
Stamp
,
L. A.
,
Foong
,
J. P. P.
,
Mcconnell
,
S. N.
,
Bergner
,
A. J.
,
Anderson
,
R. B.
,
Enomoto
,
H.
,
Newgreen
,
D. F.
,
Obermayr
,
F.
,
Furness
,
J. B.
et al. 
(
2013
).
Transplanted progenitors generate functional enteric neurons in the postnatal colon
.
J. Clin. Invest.
123
,
1182
-
1191
.
Joseph
,
N. M.
,
He
,
S.
,
Quintana
,
E.
,
Kim
,
Y.-G.
,
Núñez
,
G.
and
Morrison
,
S. J.
(
2011
).
Enteric glia are multipotent in culture but primarily form glia in the adult rodent gut
.
J. Clin. Invest.
121
,
3398
-
3411
.
Kato
,
M.
,
Onoyama
,
I.
,
Yoshida
,
S.
,
Cui
,
L.
,
Kawamura
,
K.
,
Kodama
,
K.
,
Hori
,
E.
,
Matsumura
,
Y.
,
Yagi
,
H.
,
Asanoma
,
K.
et al. 
(
2020
).
Dual-specificity phosphatase 6 plays a critical role in the maintenance of a cancer stem-like cell phenotype in human endometrial cancer
.
Int. J. Cancer
147
,
1987
-
1999
.
Kuil
,
L. E.
,
Mackenzie
,
K. C.
,
Tang
,
C. S.
,
Windster
,
J. D.
,
Le
,
T. L.
,
Karim
,
A.
,
De Graaf
,
B. M.
,
Van Der Helm
,
R.
,
Van Bever
,
Y.
,
Sloots
,
C. E. J.
et al. 
(
2021
).
Size matters: large copy number losses in Hirschsprung disease patients reveal genes involved in enteric nervous system development
.
PLoS Genet.
17
,
e1009698
.
Langhe
,
R. P.
,
Gudzenko
,
T.
,
Bachmann
,
M.
,
Becker
,
S. F.
,
Gonnermann
,
C.
,
Winter
,
C.
,
Abbruzzese
,
G.
,
Alfandari
,
D.
,
Kratzer
,
M.-C.
,
Franz
,
C. M.
et al. 
(
2016
).
Cadherin-11 localizes to focal adhesions and promotes cell–substrate adhesion
.
Nat. Commun.
7
,
10909
.
Laranjeira
,
C.
,
Sandgren
,
K.
,
Kessaris
,
N.
,
Richardson
,
W.
,
Potocnik
,
A.
,
Berghe
,
P. V.
and
Pachnis
,
V.
(
2011
).
Glial cells in the mouse enteric nervous system can undergo neurogenesis in response to injury
.
J. Clin. Invest.
121
,
3412
-
3424
.
Lee
,
J. U.
,
Huang
,
S.
,
Lee
,
M. H.
,
Lee
,
S. E.
,
Ryu
,
M. J.
,
Kim
,
S. J.
,
Kim
,
Y. K.
,
Kim
,
S. Y.
,
Joung
,
K. H.
,
Kim
,
J. M.
et al. 
(
2012
).
Dual specificity phosphatase 6 as a predictor of invasiveness in papillary thyroid cancer
.
Eur. J. Endocrinol.
167
,
93
-
101
.
Li
,
G.
,
Yu
,
M.
,
Lee
,
W.-W.
,
Tsang
,
M.
,
Krishnan
,
E.
,
Weyand
,
C. M.
and
Goronzy
,
J. J.
(
2012
).
Decline in miR-181a expression with age impairs T cell receptor sensitivity by increasing DUSP6 activity
.
Nat. Med.
18
,
1518
-
1524
.
Liao
,
Y.
,
Smyth
,
G. K.
and
Shi
,
W.
(
2014
).
featureCounts: an efficient general purpose program for assigning sequence reads to genomic features
.
Bioinformatics
30
,
923
-
930
.
Love
,
M. I.
,
Huber
,
W.
and
Anders
,
S.
(
2014
).
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol.
15
,
550
.
Lu
,
B. C.
,
Cebrian
,
C.
,
Chi
,
X.
,
Kuure
,
S.
,
Kuo
,
R.
,
Bates
,
C. M.
,
Arber
,
S.
,
Hassell
,
J.
,
Macneil
,
L.
,
Hoshi
,
M.
et al. 
(
2009
).
Etv4 and Etv5 are required downstream of GDNF and Ret for kidney branching morphogenesis
.
Nat. Genet.
41
,
1295
-
1302
.
Mallon
,
B. S.
,
Shick
,
H. E.
,
Kidd
,
G. J.
and
Macklin
,
W. B.
(
2002
).
Proteolipid promoter activity distinguishes two populations of NG2-positive cells throughout neonatal cortical development
.
J. Neurosci.
22
,
876
-
885
.
May-Zhang
,
A. A.
,
Tycksen
,
E.
,
Southard-Smith
,
A. N.
,
Deal
,
K. K.
,
Benthal
,
J. T.
,
Buehler
,
D. P.
,
Adam
,
M.
,
Simmons
,
A. J.
,
Monaghan
,
J. R.
and
Matlock
,
B. K.
et al. 
(
2021
).
Combinatorial transcriptional profiling of mouse and human enteric neurons identifies shared and disparate subtypes in situ
.
Gastroenterology
160
,
755
-
770.e26
.
Mckeown
,
S. J.
,
Mohsenipour
,
M.
,
Bergner
,
A. J.
,
Young
,
H. M.
and
Stamp
,
L. A.
(
2017
).
Exposure to GDNF enhances the ability of enteric neural progenitors to generate an enteric nervous system
.
Stem Cell Rep.
8
,
476
-
488
.
Mederer
,
T.
,
Schmitteckert
,
S.
,
Volz
,
J.
,
Martínez
,
C.
,
Röth
,
R.
,
Thumberger
,
T.
,
Eckstein
,
V.
,
Scheuerer
,
J.
,
Thöni
,
C.
,
Lasitschka
,
F.
et al. 
(
2020
).
A complementary study approach unravels novel players in the pathoetiology of Hirschsprung disease
.
PLoS Genet.
16
,
e1009106
.
Mignone
,
J. L.
,
Kukekov
,
V.
,
Chiang
,
A.-S.
,
Steindler
,
D.
and
Enikolopov
,
G.
(
2004
).
Neural stem and progenitor cells in nestin-GFP transgenic mice
.
J. Comp. Neurol.
469
,
311
-
324
.
Molina
,
G.
,
Vogt
,
A.
,
Bakan
,
A.
,
Dai
,
W.
,
De Oliveira
,
P. Q.
,
Znosko
,
W.
,
Smithgall
,
T. E.
,
Bahar
,
I.
,
Lazo
,
J. S.
and
Day
,
B. W.
(
2009
).
Zebrafish chemical screening reveals an inhibitor of Dusp6 that expands cardiac cell lineages
.
Nat. Chem. Biol.
5
,
680
-
687
.
Morarach
,
K.
,
Mikhailova
,
A.
,
Knoflach
,
V.
,
Memic
,
F.
,
Kumar
,
R.
,
Li
,
W.
,
Ernfors
,
P.
and
Marklund
,
U.
(
2021
).
Diversification of molecularly defined myenteric neuron classes revealed by single-cell RNA sequencing
.
Nat. Neurosci.
24
,
34
-
46
.
Morrison
,
J. A.
,
Mclennan
,
R.
,
Wolfe
,
L. A.
,
Gogol
,
M. M.
,
Meier
,
S.
,
Mckinney
,
M. C.
,
Teddy
,
J. M.
,
Holmes
,
L.
,
Semerad
,
C. L.
,
Box
,
A. C.
et al. 
(
2017
).
Single-cell transcriptome analysis of avian neural crest migration reveals signatures of invasion and molecular transitions
.
Elife
6
,
e28415
.
Mwizerwa
,
O.
,
Das
,
P.
,
Nagy
,
N.
,
Akbareian
,
S. E.
,
Mably
,
J. D.
and
Goldstein
,
A. M.
(
2011
).
Gdnf is mitogenic, neurotrophic, and chemoattractive to enteric neural crest cells in the embryonic colon
.
Dev. Dyn.
240
,
1402
-
1411
.
Nagy
,
N.
and
Goldstein
,
A. M.
(
2006
).
Endothelin-3 regulates neural crest cell proliferation and differentiation in the hindgut enteric nervous system
.
Dev. Biol.
293
,
203
-
217
.
Nagy
,
N.
,
Mwizerwa
,
O.
,
Yaniv
,
K.
,
Carmel
,
L.
,
Pieretti-Vanmarcke
,
R.
,
Weinstein
,
B. M.
and
Goldstein
,
A. M.
(
2009
).
Endothelial cells promote migration and proliferation of enteric neural crest cells via β1 integrin signaling
.
Dev. Biol.
330
,
263
-
272
.
Nagy
,
N.
,
Burns
,
A. J.
and
Goldstein
,
A. M.
(
2012
).
Immunophenotypic characterization of enteric neural crest cells in the developing avian colorectum
.
Dev. Dyn.
241
,
842
-
851
.
Nagy
,
N.
,
Barad
,
C.
,
Graham
,
H. K.
,
Hotta
,
R.
,
Cheng
,
L. S.
,
Fejszak
,
N.
and
Goldstein
,
A. M.
(
2016
).
Sonic hedgehog controls enteric nervous system development by patterning the extracellular matrix
.
Development
143
,
264
-
275
.
Nagy
,
N.
,
Barad
,
C.
,
Hotta
,
R.
,
Bhave
,
S.
,
Arciero
,
E.
,
Dora
,
D.
and
Goldstein
,
A. M.
(
2018
).
Collagen 18 and agrin are secreted by neural crest cells to remodel their microenvironment and regulate their migration during enteric nervous system development
.
Development
145
,
dev160317
.
Nagy
,
N.
,
Guyer
,
R. A.
,
Hotta
,
R.
,
Zhang
,
D.
,
Newgreen
,
D. F.
,
Halasy
,
V.
,
Kovacs
,
T.
and
Goldstein
,
A. M.
(
2020
).
RET overactivation leads to concurrent Hirschsprung disease and intestinal ganglioneuromas
.
Development
147
,
dev190900
.
Nagy
,
N.
,
Kovacs
,
T.
,
Stavely
,
R.
,
Halasy
,
V.
,
Soos
,
A.
,
Szocs
,
E.
,
Hotta
,
R.
,
Graham
,
H.
and
Goldstein
,
A. M.
(
2021
).
Avian ceca are indispensable for hindgut enteric nervous system development
.
Development
148
,
dev199825
.
Natarajan
,
D.
,
Marcos-Gutierrez
,
C.
,
Pachnis
,
V.
and
De Graaff
,
E.
(
2002
).
Requirement of signalling by receptor tyrosine kinase RET for the directed migration of enteric nervous system progenitor cells during mammalian embryogenesis
.
Development
129
,
5151
-
5160
.
Navoly
,
G.
and
Mccann
,
C. J.
(
2021
).
Dynamic integration of enteric neural stem cells in ex vivo organotypic colon cultures
.
Sci. Rep.
11
,
15889
.
Nishiyama
,
C.
,
Uesaka
,
T.
,
Manabe
,
T.
,
Yonekura
,
Y.
,
Nagasawa
,
T.
,
Newgreen
,
D. F.
,
Young
,
H. M.
and
Enomoto
,
H.
(
2012
).
Trans-mesenteric neural crest cells are the principal source of the colonic enteric nervous system
.
Nat. Neurosci.
15
,
1211
-
1218
.
Olden
,
T.
,
Akhtar
,
T.
,
Beckman
,
S. A.
and
Wallace
,
K. N.
(
2008
).
Differentiation of the zebrafish enteric nervous system and intestinal smooth muscle
.
Genesis
46
,
484
-
498
.
Rao
,
M.
,
Nelms
,
B. D.
,
Dong
,
L.
,
Salinas-Rios
,
V.
,
Rutlin
,
M.
,
Gershon
,
M. D.
and
Corfas
,
G.
(
2015
).
Enteric glia express proteolipid protein 1 and are a transcriptionally unique population of glia in the mammalian nervous system
.
Glia
63
,
2040
-
2057
.
Sidebotham
,
E. L.
,
Woodward
,
M. N.
,
Kenny
,
S. E.
,
Lloyd
,
D. A.
,
Vaillant
,
C. R.
and
Edgar
,
D. H.
(
2002
).
Localization and endothelin-3 dependence of stem cells of the enteric nervous system in the embryonic colon
.
J. Pediatr. Surg.
37
,
145
-
150
.
Simpson
,
M. J.
,
Zhang
,
D. C.
,
Mariani
,
M.
,
Landman
,
K. A.
and
Newgreen
,
D. F.
(
2007
).
Cell proliferation drives neural crest cell invasion of the intestine
.
Dev. Biol.
302
,
553
-
568
.
Song
,
H.
,
Wu
,
C.
,
Wei
,
C.
,
Li
,
D.
,
Hua
,
K.
,
Song
,
J.
,
Xu
,
H.
,
Chen
,
L.
and
Fang
,
L.
(
2015
).
Silencing of DUSP6 gene by RNAi-mediation inhibits proliferation and growth in MDA-MB-231 breast cancer cells: an in vitro study
.
Int. J. Clin. Exp. Med.
8
,
10481
-
10490
.
Southwell
,
B. R.
(
2006
).
Staging of intestinal development in the chick embryo
.
Anat. Rec. A Discov. Mol. Cell. Evol. Biol.
288A
,
909
-
920
.
Stavely
,
R.
,
Bhave
,
S.
,
Ho
,
W. L. N.
,
Ahmed
,
M.
,
Pan
,
W.
,
Rahman
,
A. A.
,
Ulloa
,
J.
,
Bousquet
,
N.
,
Omer
,
M.
,
Guyer
,
R.
et al. 
(
2021
).
Enteric mesenchymal cells support the growth of postnatal enteric neural stem cells
.
Stem Cells
39
,
1236
-
1252
.
Stavely
,
R.
,
Hotta
,
R.
,
Picard
,
N.
,
Rahman
,
A. A.
,
Pan
,
W.
,
Bhave
,
S.
,
Omer
,
M.
,
Ho
,
W. L. N.
,
Guyer
,
R. A.
and
Goldstein
,
A. M.
(
2022
).
Schwann cells in the subcutaneous adipose tissue have neurogenic potential and can be used for regenerative therapies
.
Sci. Transl. Med.
14
,
eabl8753
.
Tian
,
D.-H.
,
Qin
,
C.-H.
,
Xu
,
W.-Y.
,
Pan
,
W.-K.
,
Zhao
,
Y.-Y.
,
Zheng
,
B.-J.
,
Chen
,
X.-L.
,
Liu
,
Y.
,
Gao
,
Y.
and
Yu
,
H.
(
2021
).
Phenotypic and functional comparison of rat enteric neural crest-derived cells during fetal and early-postnatal stages
.
Neural Regen. Res.
16
,
2310
.
Tirosh
,
I.
,
Izar
,
B.
,
Prakadan
,
S. M.
,
Wadsworth
,
M. H.
,
Treacy
,
D.
,
Trombetta
,
J. J.
,
Rotem
,
A.
,
Rodman
,
C.
,
Lian
,
C.
and
Murphy
,
G.
(
2016
).
Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq
.
Science
352
,
189
-
196
.
Uesaka
,
T.
and
Enomoto
,
H.
(
2010
).
Neural precursor death is central to the pathogenesis of intestinal aganglionosis in ret hypomorphic mice
.
J. Neurosci.
30
,
5211
.
Uesaka
,
T.
,
Nagashimada
,
M.
and
Enomoto
,
H.
(
2015
).
Neuronal differentiation in schwann cell lineage underlies postnatal neurogenesis in the enteric nervous system
.
J. Neurosci.
35
,
9879
-
9888
.
Watanabe
,
Y.
,
Broders-Bondon
,
F.
,
Baral
,
V.
,
Paul-Gilloteaux
,
P.
,
Pingault
,
V.
,
Dufour
,
S.
and
Bondurand
,
N.
(
2013
).
Sox10 and Itgb1 interaction in enteric neural crest cell migration
.
Dev. Biol.
379
,
92
-
106
.
Young
,
H. M.
,
Ciampoli
,
D.
,
Hsuan
,
J.
and
Canty
,
A. J.
(
1999
).
Expression of Ret-, p75NTR-, Phox2a-, Phox2b-, and tyrosine hydroxylase-immunoreactivity by undifferentiated neural crest-derived cells and different classes of enteric neurons in the embryonic mouse gut
.
Dev. Dyn.
216
,
137
-
152
.
Young
,
H. M.
,
Bergner
,
A. J.
,
Anderson
,
R. B.
,
Enomoto
,
H.
,
Milbrandt
,
J.
,
Newgreen
,
D. F.
and
Whitington
,
P. M.
(
2004
).
Dynamics of neural crest-derived cell migration in the embryonic mouse gut
.
Dev. Biol.
270
,
455
-
473
.
Young
,
H. M.
,
Bergner
,
A. J.
,
Simpson
,
M. J.
,
Mckeown
,
S. J.
,
Hao
,
M. M.
,
Anderson
,
C. R.
and
Enomoto
,
H.
(
2014
).
Colonizing while migrating: how do individual enteric neural crest cells behave?
BMC Biol.
12
,
23
.
Zuchegna
,
C.
,
Di Zazzo
,
E.
,
Moncharmont
,
B.
and
Messina
,
S.
(
2020
).
Dual-specificity phosphatase (DUSP6) in human glioblastoma: epithelial-to-mesenchymal transition (EMT) involvement
.
BMC Res. Notes
13
,
374
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information