The enteric NCCs are a multipotent cell population that originates in the neural tube and migrates throughout the embryo, proliferates and eventually differentiates to give rise to the neurons and glia of the enteric nervous system . The development of ENS is complex directed by a number of cellular and molecular processes. To better understand the molecular mechanisms of enteric crest differentiation from precursors to a mature, fully functional ENS, we generated a transcriptomic profile of the zebrafish ENS at 7dpf. By 7dpf, the neural crest cells have already migrated into and populated the larval zebrafish intestine. At this stage in development, enteric NCCs have differentiated or are in the process of differentiating into the neurons and glia of the ENS [49, 51]. The composition of the ENS is likely to be a combination of fully differentiated enteric neurons and glia as well as enteric neural crest cells undergoing proliferation and differentiation. Using a transgenic line marking the phox2b expressing neural crest cells we aimed to build a lineage-specific transcriptional profile of progenitors and progenitor derivatives of the ENS.
Previous microarray experiments comparing animals with and without a complete ENS or cultured neurospheres, and RNA-seq expression analysis on culture-propagated human embryonic stem cells, or anterior intestinal tracts have revealed a number of genes important for ENS transcriptome development [1, 95–97]. However, here we expand these previous attempts to obtain a more comprehensive unbiased analysis of the enteric transcriptome. Our dataset includes low abundance and variable transcripts expressed in the ENS along the entire intestine during normal development in zebrafish. This expands our knowledge of the vertebrate ENS in itself, and opens up molecular comparisons between the ENS in higher and lower vertebrates by providing in depth gene expression profiles of enteric neural crest and their derivatives in this species. To complement this analysis we also generated a transcriptomic profile of the cells surrounding the ENS, the cellular microenvironment comprised by the endo- and mesodermally-derived gut tissue. Importantly, to minimize the need to amplify the RNA prior to sequencing and potential amplification bias prior to sequencing we dissected and FAC sorted over 4700 larval intestines. Our initial analyses suggest that this strategy allowed us to accurately and reproducibly determine transcript abundance over a wide dynamic range throughout the genome.
Our experimental strategy included using the zebrafish Tg(phox2b:EGFP) w37 line to mark GFP-positive ENS cells, the manual dissection of the intestine to separate GFP-positive ENS cells from other GFP-positive cells, the physical dissociation and FAC sorting of the intestinal cells into the GFP-positive ENS/neuronal cells, and GFP-negative non-neuronal cells, RNA-sequencing of samples with two biological replicates, and bioinformatics procedures to map reads to the zebrafish genome, and determine read counts for individual genes genome-wide. The success of our experimental strategy was validated by demonstrating that (1) ectopic GFP and endogenous phox2b transcripts almost exclusively enriched in the GFP-positive cell population, (2) housekeeping genes were similarly expressed in both cell populations, and (3) RNA-seq and qPCR results consistently showed a strong enrichment of candidate neuronal genes (phoxb, etc) and candidate non-neuronal genes (amy2, fli-1) in the expected cell populations, respectively. Furthermore, subsequent bioinformatics analyses demonstrated that expression of genes with GO terms associated with ‘neuronal’ features were enriched in the GFP-positive cell population, and genes with GO terms associated with ‘non-neuronal’ features were enriched in the GFP-negative cell population (discussed in more detail below). Thus, our datasets yielded meaningful transcriptional profiles to identify the abundance of specific ENS and non-neuronal intestinal genes genome-wide.
The GFP-positive neuronal cell population: Profiles of the ENS
Taking an RNA-seq approach, we isolated the enteric neurons from the zebrafish intestine and compared gene expression profiles between the neuronal and non-neuronal populations by carrying out a differential gene expression analysis. We identified 2561 genes (Additional file 2) that showed enrichment in the transcriptome of the enteric neuronal population as compared to the non-neuronal population of cells. These represent a diverse set of molecular and cellular genes and pathways that were further examined by various ‘pathway’ analyses.
The top 20 pathways showing the most significant enrichment activity in the neuronal population included ‘G protein-coupled receptor signaling’, ‘axon guidance’, and several well-studied neurotransmitter signaling pathways. G protein-coupled receptors have been found to be very important in neurogenesis and anomalies in members of this superfamily are known to cause aganglionosis in the intestine of mammals [98–100]. Other pathways that showed high activity in the enteric population included many of the well-studied neurotransmitter signaling pathways, the serotonin receptor, GABA receptor and NOS signaling. Serotonin receptors have been associated with abnormal motility of the gastrointestinal tract of mammals [101–103], and also associated with such diseases as inflammatory bowel disease and irritable bowel syndrome .
To find the specific functional, molecular and structural gene categories, we performed a GO enrichment analysis. The GO enrichment pathway analysis can be viewed in a hierarchy of GO terms where several terms are nested under other terms. The complete analysis gave rise to an extended hierarchy that consisted of hundreds of nodes, for example, the biological process hierarchy of the depleted list consisted of 280 nodes and 467 edges. Subsequent manual editing created a second core hierarchy reducing unnecessary complexity for easier visualization. Expectedly, our analysis revealed terms such as ‘neurogenesis’, and ‘axonogenesis’, that were enriched within the GFP-positive neuronal population. We also observed an enrichment of molecules characteristic of neurons including signal transduction, metal ion transport, neuron differentiation and cell adhesion. For example, genes such as ROBO4, SLIT2, NADL1.1, PLXN, etc. (see Additional file 2) were identified that play key roles in cell migration and axon guidance during nervous system development.
While there were common terms in both analyzed intestinal cell populations, the significance level showed clear differentiation between the two groups. For example, the term nervous system development has a p-value of 2.62E-10 in the hierarchy specific to the neuronal population while has p-value lower than <5.00E-10 in the non-neuronal GO enrichment hierarchy. Thus, the term ‘nervous system development’ is much more significant in the neuronal population as compared to the non-neuronal population consistent with a higher enrichment of genes associated with neurogenesis and neuron development in the GFP-positive population. As expected, similar terms associated with neuron differentiation, axonogenesis, neurogenesis, etc. are also found enriched in this hierarchy.
When a subset of the GO terms was grouped according to molecule type (Fig. 5) further distinct differences could be seen. The neuronal population had an overall higher number of different types of transporters and ion-channels. Ion-channel genes were the most significant increase in molecule type in the neuronal population. 8% of the genes that fell into the ‘ion-channel’ category, were genes associated with voltage-dependent calcium, potassium or sodium channels reflective of their role in generating and modulating neural activity.
Axon guidance molecules in the ENS
As enteric neuron precursors are differentiating and forming neural networks within the intestine at this stage in development it is not surprising that ‘axonal guidance signaling’, was also enriched in the neuronal dataset. Among the axon guidance molecules indicated as enriched by the pathway analysis were the Sema, Plxn, Nadl1, Fxd3 and Slit genes.
The Semaphorin (sema) family’s role in axon guidance is well documented  but recent evidence demonstrates that semaphorins play also an important role in ENS development [106–109]. In zebrafish the knockdown of Sema3C and Sema3D results in the loss of enteric neurons in a dose dependent manner [109–111]. GWAS studies have also identified polymorphisms at the Sema3C and Sema3D loci as possible at-risk alleles for HSCR .
Interestingly, our analysis suggests a novel role for the Nogo/Reticulon pathway in ENS development. The Nogo/Reticulon pathway has been implicated in the inhibition of axon growth cones and myelination . Our data demonstrates that Reticulon4 receptor (RTN4R), also known as the Nogo receptor, is upregulated in the ENS population. Concomitantly, leucine-rich repeat Ig domain-containing Nogo-interacting protein 1 (LINGO1) known to mediate the collapse of growth cones in the presence of certain myelin proteins  is also upregulated in our data. RTNR4 and LINGO1 form a complex with the p75neurotrophin receptor to mediate the inhibition of growth cones . Although a role for this axon guidance pathway has not been implicated in ENS development, reticulon4 (RTN4/Nogo) is expressed in the mammalian adult ENS . It will be exciting to see whether further studies can confirm definitive roles for LINGO1, and RTNK4 in ENS development.
In addition to the axonal guidance molecules, other gene families found within the neuronal dataset, the ADAM and neurotrophic tyrosine receptor kinase (NTRK) families are known to have a role in axon guidance. A disintegrin and metalloprotease (ADAM) is from the Metzincins superfamily of metalloproteases . This family has been shown to play a very important role in development by regulating cell migration, differentiation, cell-cell interaction, and receptor-ligand signaling [117, 118]. Interestingly, some of the axonal guidance pathway genes enhanced in our dataset have been implicated in the development of the nervous system. For example, ADAM22 is necessary in PNS development and deficiency leads to hypomyelination of peripheral nerves and ataxia . ADAM22 deficient mice display defects in the proliferation and differentiation of glia . The neurotrophic tyrosine receptor kinase (NTRK) family comprises of receptors that are required to maintain synaptic strength and plasticity in the nervous system . The three genes associated with axonal guidance from this family are NTRK1, NTRK2 and NTRK3. Similar to the ADAM family, all of these genes are highly expressed in the neuronal population in our dataset.
From our analyses we have identified both genes known to be expressed in and required for enteric neuron development, as well as candidate genes that may play a role in, enteric neuron differentiation, enteric axon guidance and connectivity. Among these genes, we found some known to be required for ENS development in other model organisms; others known to be required for neuronal development in the CNS but not directly tested for a role in ENS development; and finally genes not known to play a role in either neuronal development either in the CNS or ENS. Based on expression levels in our dataset, the latter two categories of genes may prove to be interesting candidates for further investigation. To illustrate these categories, we selected a few genes, with different molecular functions, from our enriched dataset and searched the literature for what was known about their role in ENS or neuronal development in general (Table 3).
Neuropeptides in the ENS cell population
Not unexpectedly, the neuropeptide class of molecules was also significantly enriched in the GFP-positive population. Vasoactive intestinal peptide (Vip), NeuromedinU (Nmu), Cocaine and amphetamine-regulated transcript2 (Cart2), Secretograinin II (Scg2b) and Neuromedin b (Nmba) are neuropeptide precursors, neuropeptides, or members of neuropeptide families, expressed during and required for ENS development and function [122–127].
Sectretoneurin, a functional neuropeptide expressed in the ENS, is derived from the secretograinin II gene [85, 86]. The CART genes have been associated with appetite and energy control in zebrafish , and with energy balance in chicken . Cart is expressed in a subset of Ret-positive enteric neurons in the murine intestine, and is reported to be one of the “earliest markers” of the ENS . In their study Heanue and Pachnis also showed that Cart positive neurons represented a subset of enteric neurons [1, 130, 131]. Duplication events in the zebrafish genome have resulted in 4 cart genes compared to a single gene, Cart propeptide (Cartpt), in mouse [132, 133]. We find that in our dataset only Cart2 seems to be enriched (Additional file 2: Table S1). Together the expression of cart neuropetides in our data and these other studies suggest the CART family of nueopepetides may play important roles in ENS development and/or function. A more definitive role for the CART family of neuropetides in ENS function await further studies.
Enrichment of synaptogenesis genes in the ENS
Synapse assembly and disassembly is an important part of the formation and maintenance of neural circuitry in the developing nervous system as well as plasticity in the mature nervous system [134, 135]. As expected we found many components in this category to be enriched in our dataset. Interestingly, our transcriptomic data may also show evidence that genes whose protein products form functional complexes together may have similar levels of enrichment. An example of this is proteins aligned with synaptogenesis and synaptic vesicle fusion (Table 3). Syntaxin 1 (Stx1), Syntaxin binding protein 1 (Stxbp1) and SNAP-25 are all components of the SNARE [(soluble NSF attachment protein) receptor] complex of presynaptic proteins that facilitate the fusion of synaptic vesicles with the plasma membrane during the release of neurotransmitters . Inhibitors of SNARE proteins have been shown to decrease enteric neural crest migration and ENS precursor neurite extension .
Complexin II (Cplx2), Synapsin II (Syn2), and Synaptotagmin 1 (Syt1) modulate vesicle fusion and neurotransmitter release. Syt1 is a synaptic vesicle transmembrane protein involved in Ca2+-sensitive triggered neurotransmitter release expressed in enteric neurons [89, 136, 137]. Syn2 a neuron-specific phosphoprotein, assists in clustering and recycling of synaptic vesicles . Cplx2 is a soluble pre-synaptic protein that binds the SNARE complex aids in the fusion of synaptic vesicle but negatively regulates vesicle fusion [137, 139]. In mammals the two paralogous genes stx1A and stx1B are often expressed in the same cells, and function often redundantly . Both Syntaxin 1A and 1B are expressed in the murine myenteric plexus but it is not clear whether they are functionally redundant as they localize to different regions of the neuron . Snap-25 is also expressed in the ENS of rodents and humans and localizes to the presynaptic membrane of cholinergic and nitrergic enteric neurons [136, 142–144]. During vertebrate evolution, the teleost fish lineage experienced an additional genome duplication [145, 146]. As with other genes this led to duplications and retention of additional stxbp1 and the snap-25 genes in the zebrafish genome. Thus, the zebrafish genome contains two stxbp1 paralogues, stxbp1a and 1b and two paralogues of snap-25, snap25a and 2b. Currently, it is not known which of these paralogues are required for vesicle docking and fusion in the zebrafish enteric neurons, however our data shows that all the major protein components of the presynaptic vesicle fusion are significantly enriched in our neuronal transcriptomic profile. Interestingly, these genes appear to be enriched at similar levels, between 4.36 and 6.17 with an average log2fold change expression of 5.083 (see Table 3) including stxbp1a and snap 25A but not stxbp1b and snap 25B. It is tempting to speculate that the observed expression levels may predict the specific paralogue contributions to a protein complex such as the vesicle fusion machinery in the ENS cells.
The DNA-and RNA binding proteins of the ENS
As transcription factors are generally required to drive the specification and differentiation of developmental processes we initially surveyed genes that have been described to be expressed in enteric neural crest, or enteric neurons. Our cell selection process was based on cells presumed to express phox2b (see above), so we examined the gene expression levels of phox2b and its paralogue phox2a. Both phox2a and phox2b were expressed at log2foldchange higher than 4. phox2b is known to drive the specification and differentiation of enteric neural crest and is required for the normal development of a complete and functional ENS [53–55]. Phox2a is also similar to Phox2b in its expression in the enteric neurons, cranial glia and other differentiating autonomic neurons [64, 147–149] but its expression appears to be more restricted to a subset of enteric neural crest . Both Phox2b and Phox2a have been shown to transcriptional regulate Tlx2 . Loss of Phox2b gives rise to aganglionic gastrointestinal tract  while in the absence of Phox2a some enteric neurons are present in the oesophogeal region . Whether more caudal enteric neurons are present or whether more subtle phenotypes such as the loss of specific subtypes of enteric neurons can be found in Phox2a deficient mice is not clear.
Further transcription factors that are enriched in the ENS cell population included T-cell leukemia homeobox 2 (Tlx2), Islet2a (Isl2), and Chromodomain helicase DNA binding protein 5 (Chd5). Interestingly, Tlx2 (Hox11l) has been shown to be a target for Phox2b  and Phox2a [90, 91]. Tlx2 deficiency leads to megacolon and myenteric neuronal hyperplasia in mice [90–92], which may explain its association with intestinal motility disorders like Intestinal neuronal dysplasia (IND). Isl2a, a LIM homeobox transcription factor, plays a critical role in lineage specific differentiation of motoneurons in the spinal cord [93, 153].
Chd5 encodes a chromatin remodeling protein that is specific to neurons  and promotes neurogenesis by activating the transcription of neural differentiation genes and inhibiting genes for non-neural fates [154–158]. Isl2a and Chd5 are required for neuronal differentiation in the CNS, but have not been implicated in ENS development, and may represent novel candidates required for ENS specification or differentiation.
At day 7 of larval development, some proportion of neural crest cells are differentiated and forming a functioning neuronal network within the ENS. Consistently we found, elavl3 and elavl4 two pan-neuronal RNA binding proteins expressed in postmitotic neurons were equally highly enriched with log2foldchanges of 6.9122 and 5.8957 respectively. This significant enrichment is consistent with increasing number of differentiated enteric neurons [159–163] at this time point in zebrafish development.
Signaling molecules within the ENS
An additional category of genes we selected to sample was signaling molecules. Ret tyrosine kinase receptor was significantly enriched in our neuronal population with a log2folchange of 3.7, consistent with its role in ENS development (see earlier discussion). Fibroblast growth factor 13 (Fgf13) is a members of the Fibroblast growth factor (FGF) family . Although Fgf13 shares structural similarities with the FGF family, it operates as an intracellular FGF that is neither secreted from the cell as other FGFs, nor does it activate the canonical FGF signaling pathway through FGF tyrosine kinase FGF receptors. Fgf13 is highly expressed in the embryonic CNS and has been shown to regulate neural development by stabilizing microtubules [165–167]. Mice deficient in FGF13 have defects in axonal branching, learning and memory. In addition to the CNS, Fgf13 expression is also present in the enteric ganglia of mice [136, 168]. In zebrafish there are two Fgf13 genes, fibroblast growth factors 13a and 13b (fgf13a and fgf13b, respectively) . fgf13b is expressed in the developing zebrafish ENS, and although not yet examined in the ENS, fgf13a is expressed in neurons of the CNS [170, 171]. Our data demonstrates that both fgf13a and fgf13b are expressed in transcriptome of the developing ENS, suggesting FGF13 signaling may play a role in enteric nervous system development or function.
Among the candidate genes identified from our transcriptomic profiling that could play a signaling role in ENS development and/or function is Chimerin1 (Chn1). Chn1 a GTPase-activating protein expressed in neurons in embryogenesis during a time of differentiation and synaptogenesis . To our knowledge, the expression of, or a functional role for Chn1 in ENS development has not been investigated. Intriguingly, a search of Mouse Genome Informatics (MGI) database and the International Mouse Phenotyping Consortium (IMPC) identified a transgenic reporter line, Chn1tm1.1(KOMP)Vlcg, driving LacZ expression under the Chn1 promoter. LacZ positive cells can be seen in the gastrointestinal tract in locations consistent with the position of the myenteric plexus [173–175] suggesting a potentially conserved role of chn1 in vertebrate ENS development.
Dickkopf 1 (Dkk1) belongs to the Dickkopf family of sectreted proteins that act as inhibitors of the Wnt-signaling pathway. Wnt signaling can induce synaptic protein clustering, promote the assembly development and maintenance of synapses [176–178]. Dkk1 inhibits Wnt causing the disassembly of synapses in mature neurons [176, 179]. It may be that an enrichment of Dkk1 expression in the enteric neuronal population data reflects a role for Dkk1 in synapse disassembly, an essential part of the wiring and refinement process within the developing ENS neural circuitry (Table 3) .
The GFP-negative non-neuronal cell population of the intestine
Although the goal in this study was to define the ENS (GFP positive) transcriptome, the additional transcriptome of the sorted GFP-negative non-neuronal cell population of the intestine serves not only as an important control for our experimental strategy (see above), but it provides a powerful complementary data set in itself. The non-neuronal tissue of the intestine is independently formed during development. Endodermal and mesodermal tissues form the digestive, vasculature, and muscle derivatives of the intestine, respectively. Thus, the GFP-negative cell population should be defined by entire different transcriptional signature.
Indeed, our GO term analysis found the expected enrichments for terms like proteolysis, blood vessel development and vasculature development to be very specific to the GFP-negative cell population of the intestine. By contrast, the term nervous system development had a much lower p-value (<5.00E-10) in the non-neuronal GO enrichment hierarchy than in the ENS cell population. Similarly, while the neuronal population had a higher overall number of different types of transporters and ion-channels, the non-neuronal cells had a higher number of transcriptional regulators and peptidase genes. We also observed high expression of ‘peptidases’ genes chymotrypsinogen B2, chymotrypsin-like, and several members of the carboxypeptidase, protease and matrix metallopeptidases families that reflect, molecularly define, and are consistent with the digestive nature of the large majority of the GFP-negative cell population. Interestingly, of the few ion channels expressed in the non-neuronal population, the cystic fibrosis transmembrane conductance regulator (CFTR) channel was significantly enriched. CFTR is expressed in the intestinal epithelium and is required to maintain water homeostasis in the intestine . The differential appearance of this gene in the non-neuronal population suggests our GFP-negative cell profile represents the intestinal microenvironment.
Among the pathways enriched in the GFP-negative non-neuronal population pathways, and downregulated in the neurons, are genes associated with the vasculature and tight junction formation and regulation (see Additional file 23: Table S4). The regulation of the ‘epithelial to mesenchymal transition’ and ‘axon guidance’ pathways were also downregulated in the neuronal population. The latter initially may seem counter intuitive, however, an examination of the genes within the ‘downregulated’ axon guidance category revealed genes like Shh and BMP, growth factors expressed in the endoderm and mesoderm of the developing intestinal tract and known to affect enteric neural crest migration and proliferation (see Additional file 23: Table S4) [181–189].
As indicated above the independently derived non-neuronal intestine provides a substrate for the arriving neuroectodermally derived ENS cell population. This cellular microenvironment must provide signals to guide migration, proliferation, penetration, settlement, and differentiation of the ENS cells along the entire intestine in a non-cell autonomous manner to ensure the full functional innervation along the entire developing gut. This process must be dependent on multiple interactions between these cell populations, and our data sets should have captured this molecular dialogue at a particular important and sensitive developmental time when this connection is being established. Thus, closer examination of this GFP-negative non-neuronal data set may provide a starting point to molecular dissect these little known interactions on a molecular level.