BACKGROUND
the enteric nervous system , the largest division of the peripheral nervous system , is composed of a network of neurons and glia that innervate the gastrointestinal tract  <cit> .
the ens functions within the gastrointestinal tract to control motility for mixing and moving food, absorption of nutrients, epithelial secretions, and blood circulation  <cit> .
it arises from a population of multipotent neural crest cells , which migrate from the vagal, and in some species sacral, axial area to the intestine.
there they continue to proliferate and differentiate into enteric neurons and glia within the gastrointestinal wall  <cit> .
the migration, specification and differentiation of the enteric neural crest are determined by intrinsic and extracellular cues.
thus, the enteric neural crest cells interact with other cells and respond to the local environments within the gut.
defects in the migration, proliferation or differentiation of the neural crest progenitors have been shown to lead to enteric neuropathies and are implicated in a variety of gastrointestinal dysmotility disorders .
previous studies have identified many important genes and pathways required for normal ens development .
best understood are genes whose role in ens development have been elucidated by loss-of-function studies or the study of congenital birth defects .
these genes encode a variety of transcription factors, cell-cell signaling molecules, neuropeptides, biosynthesis enzymes, cell adhesion molecules, soluble growth factors, extracellular matrix and cytoskeletal proteins.
categories of genes functionally linked to ens development have been expertly collated and described in a number of recent reviews .
among the most studied genes in ens development, are the genes required for the ret signaling pathway.
mutations in the ret tyrosine kinase receptor gene, it’s co-receptor gdnf family receptor alpha  <dig> , and it’s ligand gdnf emerged as the most common causes of hirschprung disease .
in hscr patients, insufficient proliferation and migration of enteric neural crest lead to partial loss of the ens generally characterized by a complete lack of enteric neurons within the caudal region of the gastrointestinal tract .
the ret signaling pathway is most commonly affected in hscr,  <cit>  accounting for 20–50% of familial and 15–30% of sporadic cases  <cit> .
the etiology of hscr has recently been extensively covered in several reviews  <cit> .
in addition to ret, more than a dozen other genes including dna -methyltransferase  <dig> beta , endothelin converting enzyme  <dig> , endothelin  <dig> , endothelin receptor type b , glial cell derived neurotrophic factor , gfrα <dig>  kinesin family member 1b , l <dig> cell adhesion molecule , neurotrophin  <dig> , neuregulin  <dig> , neuregulin  <dig> , neurturin , neurotrophic receptor tyrosine kinase  <dig> , phox2b, prokineticin  <dig> , prokineticin receptor  <dig> , prokineticin receptor  <dig> , persephin , semaphorin 3a , semaphorin 3d , sry-box  <dig> , transcription factor  <dig> , and zinc finger e-box binding homeobox  are known to be associated with, or cause hscr .
yet, the exact genetic cause of 50% of hscr cases remain unknown  <cit> .
as the enteric neural progenitors migrate into the gut, they proliferate and differentiate into at least  <dig> different neuronal subtypes that differ in their function, physiology, morphology, and neurotransmitter expression  <cit> .
in mammals, these subtypes form ganglia that are organized into a network of two interconnecting plexi.
the myenteric plexus is embedded between the circular and longitudinal smooth muscle layers, and the submucosal plexus is situated between the muscle and the mucosal layer.
each neuronal subtype expresses different combinations of neuropeptides giving them a unique neurochemical code  <cit> .
along the length of the intestinal tract, the combinations of neuronal subtypes vary between the different regions of the gut  <cit> .
although the diversity of enteric neuronal subtypes and their neurochemical code has been well described, our understanding of the genes required for enteric neuronal specification, differentiation and for establishing neural connectivity within the ens remains rather limited  <cit> .
to fill this gap we pursued a broad strategy to identify candidate genes required to drive neural crest precursors to differentiate into enteric neurons including potential signaling molecules for the formation of neuronal subtypes and their integration into a functional neural network in zebrafish.
zebrafish embryos and larvae have emerged in recent years as a powerful complementary vertebrate model for ens development.
although some important anatomical differences to the ens in mammals exist e.g.
the neurons and glia of the zebrafish myenteric plexus do not form ganglia  <cit>  and the submucosal plexus is absent, the signaling pathways and genes associated with ens development remain conserved .
markers of differentiated enteric neurons can first be detected in the developing intestine approximately 48 hours post fertilization .
over the next 24 hours, neural crest cells continue to migrate and proliferate and by  <dig> hpf the entire length of the gut is colonized.
the first contractions within the intestine are seen at 4 days post fertilization  and by  <dig> dpf the larval gut is functional and the larvae are capable of feeding  <cit> .
in this study, 7 days post fertilization  larvae from a transgenic zebrafish line were used to isolate enteric neuronal cells from the intestine.
7dpf represented the time in development at which a subset of neurons has fully differentiated; yet the ens is not fully developed and likely to retain some precursors.
to track and distinguish precursor cells and differentiating enteric neurons we made use of a zebrafish transgenic line in which egfp is expressed under the regulatory/promoter region of the transcription factor phox2b, tg
w <dig>  <cit> .
this strategy was chosen because phox2b is required for the specification and development of the neural crest-derived autonomic nervous system including the enteric neurons and glia in various model systems .
loss of phox2b during development leads to loss of enteric ganglia in mice  <cit>  while in heterozygous animals a normal gi tract is observed  <cit> .
in humans, mutations in phox2b are linked to a loss of enteric neurons in the gastrointestinal tract leading to hscr, and many cases have been associated with congenital central hypoventilation syndrome  <cit> .
determination of the genetic processes underlying the syndrome revealed that mutations of the phox2b gene are principally responsible for the broad range of symptoms encountered by disrupting early development of autonomic neurons .
importantly, phox2b has been shown to regulate ret expression  <cit>  and the genetic interaction between ret and phox2b has been demonstrated to be critical for normal ens development  <cit> .
in mice phox2b expression is turned on as the neural crest cells  leave the neural tube, remains expressed as neural crest migrate into the intestine and differentiate into neurons and glia and continues to be expressed in ens derivatives into adulthood  <cit> .
our goal was to generate a comprehensive transcriptional profile of enteric neurons along the entire intestinal tract during normal development.
to improve the chances of obtaining discrete transcript counts and unbiased detection of novel or low-abundance transcripts across a broader dynamic range  <cit> , endogenous enteric neurons and enteric neural crest progenitors were transcriptionally profiled using rna-seq.
complimentary to profiling the enteric neurons, we also utilized rna-seq to transcriptionally profile their surrounding non- neural microenvironment  <cit> , the intestine including the mucosa, musculature and associated vasculature.
thus, this strategy enabled us to perform a comparative gene expression study between the neuronal cell and the non- neuronal cell population of the intestine.
our analyses identified scores of candidate genes that may function as cell fate determinants during ens development, may constitute useful markers to distinguish different ens subtypes, and define the surrounding cellular microenvironment of the ens on a molecular level.
a subset of the genes has been validated by expression studies, and selected for ongoing functional studies.
the result of these studies will not only expand our knowledge of genes required for normal ens development and a better understanding of the cellular and molecular basis of the complex process of development of the ens but may help inform stem cell therapies aimed at directing lineage specific enteric neuron development.
methods
animal husbandry and line
the embryos and larvae were obtained from a transgenic zebrafish line as described previously and provided by dr. a nechiporuk , in which green fluorescent protein  is expressed under the regulatory elements of the transcription factor phox2b tg
w <dig>  <cit> .
the embryos and larvae were raised at  <dig>  °c in a light/dark cycle according to standard zebrafish husbandry protocols and staged to seven days post fertilization   <cit> .
larvae were prescreened for gfp expression and anesthetized prior to dissection.
all animal experiments were performed in accordance with the u.s. laws, guidelines and policies of laboratory animal care and use and approved by the institutional animal care and use committee and institutional biosafety committee at iowa state university.
tissue collection and rna isolation
to acquire the enteric neuronal population, larvae from the transgenic line tg
w <dig> were dissected at 7 days post fertilization  to isolate the entire gut.
guts were dissociated into cells using papain digestion.
for dissociation, 20 μl of  <dig> mgp/ml papain was added to the dissected guts along with 360 μl of a solution of hbss,  <dig> % of hepes and 20 μl of an activation agent according to the manufacturer’s instructions .
the samples were incubated at 37 °c for a total of 30 minute with trituration of the sample every 10 minutes.
a final incubation was done for 13 minute with 20 μl of dnase i. the samples were resuspended in 500 μl of  <dig> % bsa in pbs.
fluorescent activated cell sorting  was done on each sample using a bd facsaria iii available at the iowa state flow cytometry facility.
gfp-positive cells were detected through a 525/50 nm bandpass filter, and simultaneously checked for viability, and sorted directly into lysis buffer .
both, the gfp positive neuronal cells and the gfp negative non-neuronal cells were collected directly into commercially available lysis buffer and stored at -80 °c until separated cell populations from multiple facs sorting experiments could be pooled and processed.
rna isolation was done using pooled samples of  <dig> – <dig>  cells and the ambion™ ‘s rnaqueous®- micro total rna isolation kit.
the protocol was followed as per kit instructions with the following modifications.
a total volume of 20 μl of elution buffer was used and the column was incubated at room temperature for 10 minutes before the final spin to collect the eluted rna.
a total of  <dig> larval intestines were dissected to collect  <dig>  gfp positive cells to obtain sufficient rna for two biological replicates of the gfp-positive cell population.
<dig>  gfp negative cells were collected from the same  <dig> larval intestine samples for the two biological replicates of the gfp-negative cell population.
the rna quality check was checked using an agilent  <dig> bioanalyzer.
rna integrity number s of  <dig> - <dig>  were measured for all rna samples.
libraries were prepared using illumina’s truseqv <dig> kit .
sequencing and raw read processing
the sequencing was done using an illumina hiseq  <dig> sequencing system at the dna facility of the iowa state university office of biotechnology.
deep sequencing of the mrna was performed with 100 bp paired-end reads.
each biological replicate was considered as a separate sample.
preprocessing of the raw reads was done using trimmomatic,  <cit>  to remove adapters and the low quality reads.
the preprocessing was done using the parameters leading: <dig> trailing: <dig> slidingwindow:4: <dig> minlen: <dig>  the removal of the bases was based on the phred score where the threshold was  <dig>  a window size of  <dig> was used which removed the bases with an average phred score of less than  <dig>  the quality of the trimmed reads was analyzed and visualized using fastx-toolkit  <cit> .
alignment, transcript abundance and differential gene expression analysis
the trimmed reads were mapped to the zv <dig> danio rerio reference genome using tophat <dig>  <cit>  which uses bowtie <dig>  <cit>  as the mapping tools.
the version zv <dig>  of gtf  was used for annotation of the mapped reads.
the expression levels of the mapped reads were obtained using featurecounts  <cit>  to specifically count the transcripts mapping to exons.
for the analyses shown in figs.
<dig> and 3a and b, we used cufflink to get the fragments per kilobase of transcript per million mapped reads   <cit> .
an r package deseq  <cit>  was used to perform the differential gene expression analysis between the gfp-positive enteric neuron gut population and the gfp-negative non-neuronal gut population.
the software considered the biological replicates from both the cell populations during the analysis to provide an expression level for the genes.
the differentially expressed genes were filtered using the false discovery rate  adjusted p-value ≤  <dig> .
the significantly upregulated and downregulated genes were chosen based on a log2foldchange value of  <dig> .
the heatmap was generated using the r function “pheatmap”.
for the correlation analysis, the spearman correlation coefficient for pairwise samples was calculated using r function “cor”.
the scatter plots were visualized using r function “plot”.
the expression value was log-transformed using log <dig> .
gene ontology enrichment and pathway analysis
a gene ontology  enrichment analysis was performed using the cytoscape plugin bingo  <cit>  with the adjusted p-value <  <dig>  where the fdr correction was done by benjamin and hochberg’s method.
pathway analysis was performed to identify pathways most enriched in the enriched and diminished genes in the neuronal population using qiagen’s ingenuity pathway analysis  .
to perform ipa pathway analysis human and mouse orthologues of the identified zebrafish ensembl genes were converted and ipa pathway analysis was performed on 44% of the genes from the complete data set.
the datasets analyzed were the list of significantly enhanced genes; and the list of significantly diminished genes.
the percentage of molecules was calculated for the analysis using a total of  <dig> genes from the enhanced list and a total of  <dig> genes from the depleted list.
qpcr
eight candidate genes were chosen for qpcr.
qpcr was performed on the same cdna libraries used for illlumina sequencing to validate the rnaseq results.
primers were designed using ncbi primer designing tool  and synthesized at the iowa state dna facility .
primers to chn <dig>  elavl <dig>  chata, phox2bb, syn2a, amy2a, fli1a, and hprt <dig> were designed to span exon-exon junction and give a single product of 150–200 bp.
individual qpcr reactions were carried out with the standard steponeplus  cycling protocol using perfecta sybr® green fastmix, rox  according to manufactures instructions.
for each qpcr run, the cycle threshold  were exported into microsoft excel for statistical analysis.
we used 2-ΔΔct method of relative quantification to estimate copy number of selected genes  <cit> .
the expression levels were normalized to zebrafish  hypoxanthine phosphoribosyltransferase  <dig>  as a reference gene  <cit> .
RESULTS
transcriptional profiling of ens cells and the zebrafish intestine
to identify genes required for the specification and differentiation of enteric neurons from neural crest progenitors we isolated and carried out rna-seq transcriptional profiling on a cell lineage expressing the transcription factor phox2b.
phox2b has been shown to be expressed in early enteric neural crest progenitors as they migrate into the gut  <cit> .
as development proceeds phox2b expression remains high in a subset of the differentiated neurons into adulthood and decreases to lower levels in the developing glia  <cit> .
phox2b expression is retained in a subset of adult ens derivatives.
therefore, a transgenic reporter line driving expression of gfp under the phox2b promoter, tg
w <dig>  <cit>  seemed a good candidate line to identify enteric neuron precursors and enteric neurons as gfp-positive cells within the intestine, and to facilitate separation of the neuronal from non-neuronal cell populations of the intestine for subsequent transcriptional profiling.
indeed, personal observations confirmed that gfp expression appears in the enteric neural crest as they migrate into the intestine and continues throughout development and into adulthood in the zebrafish tg
w <dig> line .
thus, gfp-positive cells of the intestine mark both enteric neurons and enteric neuron precursors in 7 days post fertilization  larvae and may facilitate separation of the ens cells from other intestine tissues by fluorescence activated cell sorting .
however, gfp expression is not restricted to ens in  <dig> dpf larvae.
endogenous phox2b, and therefore gfp, in the zebrafish tg
w <dig> line is also expressed in cells of the hindbrain and cranial ganglia as seen in fig.
1a.
to ensure that only gfp-positive cells  that comprise the ens were profiled, the intestinal tracts of 7dpf zebrafish larvae were manually dissected away from the remaining larval body prior to dissociation and facs sorting .fig.
<dig> image and flow chart illustrating source of gfp labeled cells, sample preparation and bioinformatics analysis.
a image of a 7 days post fertilization  zebrafish larvae of the transgenic line, tg.
the boxed area shows the intestine with gfp-positive enteric neurons .
b schematic diagram illustrating the experimental approach.
intestines were dissected out to avoid contamination with gfp-positive cells expressing phox2b in the hindbrain or the spinal cord.
the intestines were dissociated using papain digestion followed by fluorescent activated cell sorting .
rna isolated from the sorted cell populations was used to construct cdna libraries.
raw reads from illumina sequencing were processed and analyzed using various bioinformatics analyses programs
at this stage in development, the enteric neuron precursors have migrated to the end of the intestine, and a sufficient number of progenitors have differentiated into different ens subtypes to form a functioning ens.
enteric neurons present at this time reflect at least some of the diversity of enteric neuron subtypes  <cit> .
after enzymatically treating the dissected intestines with papain, the dissociated cells were fac sorted and viable gfp positive cells were collected .
cells from whole intestines of wild type zebrafish were used as a negative control to establish a threshold for gfp expression.
viable gfp-negative cells representing the remaining cell types in the intestine were collected simultaneously and processed in the same manner as the gfp-positive cells.
to minimize potential changes in gene expression patterns, to obtain sufficient amounts of rna to sequence, and to limit amplification steps, dissection times were limited and multiple independent cell sorts were pooled.
a total of  <dig> larval intestines were dissected to collect  <dig>  gfp positive cells to obtain sufficient rna for two biological replicates of the gfp-positive cell populations.
<dig>  gfp negative cells were collected from the same larval intestine samples for the two biological replicates of the gfp-negative cell populations.
libraries were generated from the pooled samples using unamplified rna.
biological duplicates of both the gfp-positive  and gfp-negative  cells were sequenced using illumina hiseq  <dig>  following sequencing, a developed pipeline was implemented to process the raw rna-seq reads .
figure 1b outlines the basic workflow from tissue collection to data processing.
differential gene expression analysis between neuronal and non-neuronal cells of the intestine
to identify genes that may be required for the differentiation and assembly of enteric neuron circuitry we looked for differences in expression profiles between the sorted neuronal  and non-neuronal  populations of intestinal cells.
for differential gene expression analysis the raw reads were trimmed using trimmomatic to remove the adapters and the low quality reads, and mapped to the annotated d.rerio reference genome  .
to validate that our experimental approach separating the ens/neuronal cells from non-neuronal cells by dissociation and fac sorting, and that our differential gene expression analysis yielded meaningful results, we tested for several predicted outcomes.
first, gfp mrna should be entirely restricted to the two neuronal gfp-positive cell samples and absent in the two non-neuronal gfp-negative samples.
indeed, high levels of gfp  were found in the two neuronal, and miniscule amounts  in the two non-neuronal samples .
therefore, dissociation and sorting successfully separated the two intestinal cell populations leading to a > <dig> fold enrichment of gfp in the neuronal versus non-neuronal samples.
second, we would expect that endogenous phox2b mrna should also be enriched in the gfp-positive cells.
indeed, endogenous phox2b was highly enriched in the gfp-positive cells , a ~ <dig> fold enrichment compared to the gfp-negative samples while high and low expressing house-keeping genes like actb, gapdh, and hypoxanthine phosphoribosyltransferase  <dig>  showed similar high and low expression in both sorted cell populations.
this demonstrates that gfp expression in the intestine in the zebrafish tg
w <dig> line is indeed labeling the phox2b-expressing ens cell population.fig.
<dig> experimental validation of reliable separation of the neuronal ens from the non-neuronal intestinal cell population.
a the bar graph  and corresponding data in table  show the expression levels in fpkm of five selected genes in the two neuronal replicate  and the two non-neuronal replicate  cell populations.
gfp transcript levels show the level of gfp transcript driven from the phox2b transgene.
the expression values of phox2bb measure the endogenous transcript levels of the phox2b gene.
expression levels for the three house-keeping genes actb <dig>  gapdh, and hprt <dig> are also shown.
qpcr validation of expression levels of  <dig> genes comparing the gfp-positive neuronal populations  versus the gfp-negative non-neuronal .
chn <dig>  elavl <dig>  chata, phox2bb, and syn2a were more highly expressed in the neuronal population, whereas amy2a and fli1a were more highly expressed in the non-neuronal population.
changes in gene expression were calculated using 2− ΔΔct method with hprt <dig> as the house- keeping gene.
the mean relative expression of the neuronal population was significantly different from that of the non-neuronal population
to verify that our estimated abundance levels for transcripts by rna-seq are consistent with the actual expression levels in our samples we determined the relative abundance levels of a few candidate genes using semi-quantitative pcr .
based on our experimental strategy we would predict that candidate genes that are known to control ens development and differentiation are entirely expressed or strongly enriched in the gfp-positive cells.
conversely, we would expect that that known genes involved in patterning and differentiation of the non-neuronal intestine should be strongly enriched in the gfp-negative samples.
indeed, relative expression levels of the candidate genes using semi-quantitative pcr  correlated with the rnaseq results, and showed the predicted differences in both cell populations .
elav like neuron-specific rna binding protein  <dig> , a rna binding protein found in post mitotic neurons and specific only to neurons, choline o-acetyltransferase a , whose gene product is required for biosynthesis of the neurotransmitter acetylcholine, phox2bb, a dna-binding protein which is associated with the development of noradrenergic neurons, chimerin  <dig> , a gene that is primarily expressed in the neurons and is involved in neuronal signal-transduction and synapsin 2a , a synaptic gene that has been involved in neurotransmitter release and synaptogenesis, were all found to be highly expressed in the neuronal population consistent with their roles in neurogenesis and synaptogenesis.
whereas amylase alpha 2a , one of the genes that is highly involved in proper digestion by helping in breakdown of polysaccharides and oligosaccharides and fli- <dig> proto-oncogene, ets transcription factor a , showed higher levels in the non-neuronal cell population.
hrpt <dig> was chosen as a more reliable house keeping gene  <cit> .
our results showed that there was a high correlation of relative expression levels between the rna-seq results and the qpcr for all eight genes.
furthermore, the expression for both groups of marker genes are almost entirely separated demonstrating that the ens, and the non-neuronal cell populations have been efficiently separated by dissociation and sorting with a minimum of cross contamination.
all of the above indicates that a genome scale differential expression analysis between gfp-positive and gfp-negative cell populations will yield meaningful transcriptional profiles of neuronal/ens cells and non-neuronal cells of the zebrafish intestine, respectively.
to do so we compared the genes expressed in the neuronal and non-neuronal cell populations of the intestine.
<dig> differentially expressed genes  were identified  by benjamin and hochberg’s test ≤  <dig>  and log2foldchange ≥  <dig> ) .
of these deg’s,  <dig> genes were found to be significantly enriched in the neuronal cell population .
for ease in this manuscript, we refer to genes expressed at significantly higher levels in the neuronal population in comparison to the non-neuronal population as ‘enriched’ and to genes most significantly lower in the neuronal population as ‘depleted’.
both the neuronal and non-neuronal cell populations showed a high level of similarity in the fpkms between duplicates prompting us to expand our analysis to compare all the genes expressed above fpkm of  <dig> in our data sets.
to analyze how reproducible the estimated gene expression levels between the biological duplicates were, we calculated the spearman correlation coefficient for all expressed genes .
a comparison of the biological replicates of the neuronal and non-neuronal cells yielded average coefficients of  <dig>  and  <dig>  respectively demonstrating a high correlation of gene expression levels between the duplicates .
the high level of concordance between biological replicates confirms the quality and consistency of our samples and dataset.fig.
<dig> differential gene expression analysis of 7dpf danio rerio enteric neurons.
a and b. plot showing the correlation between the two biological replicates of the gfp-positive neuronal population  and of the gfp-negative non-neuronal population .
the x-axis represents the 1st biological replicate while the y-axis represents the 2nd biological replicate.
the correlation coefficient was calculated using fpkm.
the correlation coefficient for the neuronal and non-neuronal cell population is  <dig>  and  <dig> , respectively.
c dispersion plot showing the log2foldchange of individual reads as compared to the mean normalized counts using deseq with the adjusted p-value ≤  <dig> .
the red dots represent the upregulated/enriched and the blue dots represent downregulated/depleted genes in the neuronal cells of the zebrafish intestine.
d. heat map of the normalized counts for the  <dig> degs between neuronal and non-neuronal cell populations.
the thresholds selected are fdr ≤  <dig>  and log2foldchange ≥  <dig> .
the first two columns represent biological replicates  of the non-neuronal cell population and the second two columns represent biological replicates  of the neuronal population.
each row denotes a single gene and their expression pattern across the two different samples.
the blue regions in the neuronal population represent the  <dig> % of the upregulated genes and the red region represents the downregulated genes.
the dendogram on the x-axis shows the grouping of same samples  while the y-axis shows clustering of  similar expression patterns
to illustrate the differentially expressed genes and their expression pattern visually, we generated a dispersion plot  heat map  showing the clustering of differentially expressed gene expression profiles for the neuronal and non-neuronal samples .
the heat map shows the relative expression of the genes based on their mapped read counts.
enriched genes include genes such as phox2b, elavl <dig>  ret and gng <dig>  previously associated with ens development  <cit> .
interestingly, hscr related genes such as gfra <dig> , ret, phox2b, and zeb <dig> that cause aganglionosis when knocked down in zebrafish or mutated in mice are also found enriched  <cit> .
other elevated genes in the neuronal cell population are known to impact the ens functionally without aganglionosis including ascl <dig>  cadherin egf lag seven-pass g-type receptor  <dig> , frizzled class receptor  <dig> , gfra <dig>  hand <dig>  hoxb <dig> , kinesin family members ), neural adhesion molecules , ntrk <dig> , solute carrier family  <dig> member  <dig> , t-cell leukemia homeobox  <dig>   <cit> .
among the genes significantly depleted in the neuronal cell population, we find genes that lead to intestinal aganglionosis endothelin converting enzyme  <dig> , indian hedgehog homolog a , integrin, beta 1b.
<dig> , itgb.
<dig>  and genes that cause alterations in ens function like distal-less homeobox 2a , erb-b <dig> receptor tyrosine kinase 3a , noggin , sonic hedgehog b , and zic family member  <dig>   <cit> .
therefore, through the de analysis we were able to identify previously characterized as well as many novel genes that are specifically expressed in the neuronal versus non-neuronal cell population verifying the quality of our samples and datasets.
these findings suggest that the identified enriched and depleted genes constitute a comprehensive and valuable resource for discovery and subsequent targeted functional analysis of enteric neuronal and non-neuronal cell types within the intestine, respectively.
go and molecular enrichment analysis within the two intestinal cell populations
to identify specific functional, molecular as well as structural gene categories within the enriched or depleted gene sets, a gene ontology  enrichment analysis was performed.
among genes associated with biological processes,  <dig> go terms were enriched within the gfp-positive, neuronal cell population, while  <dig> go terms were enriched among genes within the gfp-negative, non-neuronal cell population.
we generated comprehensive go term hierarchies where each node represents the number of genes associated with a particular go term, and the intensity of the orange to yellow to white color represents the significance from high to low, respectively.
biological process hierarchies were generated for both neuronal  and non-neuronal populations .
similar analyses were performed for molecular function and cellular component go categories in both the neuronal and non-neuronal populations , additional files 12–14: figures s11-s <dig> , additional files 15–17: figures s14-s <dig> , and additional files 18–20: figure s17-s <dig> .
some of the most significant go term nodes were selected to generate a core hierarchy for easier visualization .
<dig> terms found in both sets of analyses were included to illustrate the differences in gene categories between the neuronal and non-neuronal population.
interestingly the results show several of the go terms and their associated genes have reciprocal expression patterns between the neuronal and non-neuronal populations.
for example, the go term ‘nervous system development’ is highly significant  among genes enriched in the neuronal cell population with a p-value of  <dig> e- <dig>  while showing little or no significance in the depleted list  .
conversely, terms like ‘organ development’  and ‘regulation of primary metabolic processes’  had low or no significance in the neuronal population but were highly significantly enriched in the non-neuronal set .
the list for the selected go terms present in the highlighted pink and green boxes as well as the terms specific to either population along with their p-values and percentage of associated genes corresponding to fig.
<dig> can be found in table  <dig>  additional file 21: table s <dig> gives the p-value and percentage of genes for the additional top  <dig> go term categories in the analyses.fig.
<dig> gene ontology  enrichment analysis of the differentially expressed genes in the neuronal cells.
go enrichment analysis was performed using the enriched and the depleted genes from the neuronal cells.
figures illustrate a subset of go terms common to both the neuronal and non-neuronal population biological properties hierarchies  and a few specific to one or the other hierarchy  and .
the hierarchies of the go terms are based on the number of genes and the significance for the category biological process .
the circles/nodes represent the different go terms, where the size of the circle represent the number of genes associated with the term, i.e., the larger the circles, the more the number of genes associated with the go term.
the color scheme represents the significance: white represents the least, and orange represents the most significant terms.
the white, yellow and orange circles represent the terms common to both the neuronal and non-neuronal networks, while the red and blue circumferences of the circles represent terms specific to either the neuronal or non-neuronal data sets respectively.
a core hierarchy of go terms defining biological processes in enriched genes comprising the ‘neuronal’ gene set.
<dig> genes are included in this analysis.
b. core hierarchy of go terms defining biological processes in the depleted genes comprising the ‘non-neuronal’ gene set.
<dig> genes are included in this analysis.
the green and pink boxes highlight examples of sub hierarchies within the larger hierarchies of biological process go terms common to both populations but differentially enriched between the two.
sub-hierarchies of biological process go terms are significantly enriched in the neuronal population, ‘ion transport’ and ‘nervous system development’ but not in the non-neuronal population .
examples of go term sub-hierarchies enriched in the non-neuronal population include regulation of ‘metabolic and biosynthetic processes’ and ‘organ development’ .
the corresponding p-values and percentage of genes associated with the go terms in the green and pink boxes as well as the specific terms in either population are shown in table 1
the ‘--‘ indicates n/a in the table
the terms specific to a single category have ' --' in both the, p-value and percentage of genes columns
in addition to shared go terms found in both the enriched and depleted gene sets, several significant go terms were identified belonging specifically to one set or the other.
in the neuronal population the go terms ‘neurogenesis’, ‘neuron development’, and ‘axonogenesis’ were highly significant as expected, but these go terms were not present in the non-neuronal population .
similarly, go terms such as ‘proteolysis’, ‘vasculature development’ and ‘blood vessel morphogenesis’ were significantly enriched in the non-neuronal but absent in the neuronal cells .
in addition, the most significantly enriched go terms of the neuronal population emphasized the differentiation states and the establishment of neural circuitry between enteric neurons.
terms significantly enriched include terms relating to ‘synaptic function’ and ‘neurotransmission’, and are consistent with the onset of a functioning ens at this stage in development.
to gain further insights into the diverse functional molecular categories between the two cell populations, a differentially expressed gene  analysis was performed using qiagen’s ingenuity pathway analysis  software.
the identified functional categories and their percentage based on the total number of genes in each population are shown in the pie-charts in fig.
<dig>  this analysis has important limitations.
the first pie-chart  shows the categories of enriched genes according to their molecular function and their corresponding percentages within the neuronal population.
here, the total number of identified genes was  <dig>  of which  <dig> were mapped to the ipa knowledge base, and  <dig> were subsequently used for the analysis shown.
the second pie-chart  shows a similar analysis in the depleted genes, with an initial input of  <dig> genes, of which  <dig> were mapped and  <dig> were found to be analysis ready.
therefore, of the  <dig> gene orthologues included in the analysis  <dig> genes were assigned to functional designated categories.
the remaining mapped genes, though equally interesting, could not be formally assigned to a functional category due to limitations in ipa category designations and were classified as ‘others’ .
when the molecular categories within the enhanced and the depleted genes were examined separately, the neuronal population showed 8% ‘ion-channels’ as compared to 1% in the non-neuronal population.
in contrast, there were 5% ‘peptidases’ in the non-neuronal population as compared to 2% in the neuronal population.
additionally, other categories such as ‘growth factors’, ‘transcriptional regulators’ and ‘cytokines’ were enriched in the non-neuronal population  as compared to the neuronal population  .
thus, this analysis categorized genes in both cell populations based on their association with biological or molecular functions indicating significant shifts in several categories and identifying associated genes involved in specification and differentiation of cell types among both cell populations.fig.
<dig> gene function category enrichment of the significantly differentially expressed genes in the neuronal and non-neuronal gene populations: pie-chart showing the most functionally enriched terms in the significantly differentially expressed genes ipa analysis.
a gene function categories enriched in the neuronal population.
b gene category functions enriched in the non-neuronal population.
different colors represent different molecular functions and the size of the pie slice, the percentage of genes representing this category.
each category represents differentially expressed genes grouped by functional category.
for example, ion-channels represent 8% of the deg in the neuronal population  but only  <dig> % of the genes in the non-neuronal population .
the surprisingly high percentage of the genes are listed as ‘other’ because they are not represented within the functional categories defined by the ipa software parameters.
number of genes associated with each category is available in additional file 22: table s3
to obtain more information on signaling pathways that can be associated with our datasets, and to gain insight into the genes comprising these pathways, additional analysis was carried out.
ipa pathway analysis was used to determine the molecular pathways that were highly represented in either the neuronal or non-neuronal dataset.
the final analysis was done using  <dig> enriched genes and  <dig> of depleted genes of the deg genes between the neuronal and non-neuronal populations.
figure  <dig> shows the top  <dig> enriched canonical pathways in the enhanced as well as the depleted list of genes.
the top  <dig> canonical pathways along with their –log and associated genes are shown in table  <dig>  the complete information about the pathways can be found in additional file 23: table s <dig>  as expected, ‘g-protein coupled receptor-signaling’ pathway, ‘camp-mediated signaling’, ‘nnos signaling’ pathway and ‘axon guidance’ were among the highly enriched pathways in the neuronal population .
these pathway predictions implicated many additional genes as up-regulated, and, consistently, the transcript levels within our neuronal transcriptome confirmed these expectations .
conversely, some of the significant pathways in the depleted list comprised of ‘il- <dig> signaling pathway’, ‘tec kinase signaling’, and ‘inhibition of matrix metalloprotease’.
the first two pathways are known to regulate the immune system .
interestingly, the third identified pathway involving matrix metalloproteases is known to induce cell migration  <cit> .
thus, the enrichment of inhibition of this pathway may suggest a mechanism for cessation of cell migration.
this potential function would be consistent with this developmental stage in which a proportion of enteric neural crest cells has completed the migratory phase, has differentiated and is functionally capable of driving intestinal motility.fig.
<dig> pathway analysis of the differentially expressed genes in the neuronal and non-neuronal cells of the zebrafish intestine.
pathway analysis of the upregulated, downregulated, and the complete list of significantly differentially expressed genes was done using ipa.
the purple bar graph represent the top  <dig> canonical pathways predicted to be most likely to be activated in the neuronal genes using the genes most significantly enriched in the neuronal population while the blue bar graph represents the top  <dig> canonical pathways predicted to be most likely to be activated in the neuronal genes using the genes most significantly downregulated in the non-neuronal population
differentially expressed genes may drive development and differentiation of the ens
in examining our transcriptomic data, we noted many genes to be highly enriched within the neuronal population compared to the non-neuronal population.
these genes could be further divided in distinct categories according to their possible function as candidates for ens development and differentiation .
examples of molecular categories represented in the data were dna and rna binding proteins, signaling molecules, neuropeptides, and proteins involved in synaptogenesis.
interestingly but not surprising there is a significant enrichment in several neuropeptide genes such as vasoactive intestinal peptide   <cit>  and secretogranin ii , both of which are implicated in ens function.
the cytokine gene, secretogranin ii , from which the neuropeptide secretoneurin is derived, exhibits the highest enriched fold change in expression in the neuronal population.
though, secretoneurin is expressed in rat and human myenteric and submucosal ganglia  <cit>  and has been implicated in the regulation of motility of the gastrointestinal tract  <cit> , no study involving loss of function in the ens has been reported  <cit> .table  <dig> potential candiate genes with their log2foldchange
phox2aa
chd5
isl2a
elavl4a
cart2
fgf13b
chn1
reta
dkk1b
apreviously described genes in the enteric nervous system
potential novel candidate genes are bolded
genes such as synapsin ii,  and snap <dig>  required for synaptic vesicle formation and fusion and neurotransmitter release  <cit> , are similarly enriched, consistent with the establishment of synapses and a functioning neural network.
the transcription factor tlx <dig>  required for ens development and function , is enriched, as is islet 2a , a transcription factor known to play a role in motoneuron specification in the cns  <cit> .
the newly identified ens related genes can be subdivided into three main categories:  genes known to be expressed by neurons outside the ens but not yet associated with ens development;  genes reported to be expressed in the ens but for which functional studies have not yet been reported; and  novel genes that have not yet been associated with expression or function during ens development.
a more detailed discussion of table  <dig> follows below.
thus, our approach successfully identified genes and pathways known to be involved in neurogenesis and synaptogenesis, functions that are expected to be present in the ens cell population.
furthermore, this suggests also that our dataset uncovers reliably novel genes or pathways that have not been implicated in ens development yet, and can be used as a comprehensive resource to start further investigations into their roles in ens development and differentiation.
discussion
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  <cit> .
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  <cit> .
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 .
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  <dig> 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
w <dig> 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  ectopic gfp and endogenous phox2b transcripts almost exclusively enriched in the gfp-positive cell population,  housekeeping genes were similarly expressed in both cell populations, and  rna-seq and qpcr results consistently showed a strong enrichment of candidate neuronal genes  and candidate non-neuronal genes  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 .
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  <dig> genes  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  <dig> 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 .
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 , and also associated with such diseases as inflammatory bowel disease and irritable bowel syndrome  <cit> .
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  <dig> nodes and  <dig> 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 robo <dig>  slit <dig>  nadl <dig> , plxn, etc.
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  <dig> e- <dig> in the hierarchy specific to the neuronal population while has p-value lower than < <dig> e- <dig> 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  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, nadl <dig>  fxd <dig> and slit genes.
the semaphorin  family’s role in axon guidance is well documented  <cit>  but recent evidence demonstrates that semaphorins play also an important role in ens development .
in zebrafish the knockdown of sema3c and sema3d results in the loss of enteric neurons in a dose dependent manner .
gwas studies have also identified polymorphisms at the sema3c and sema3d loci as possible at-risk alleles for hscr  <cit> .
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  <cit> .
our data demonstrates that reticulon <dig> receptor , also known as the nogo receptor, is upregulated in the ens population.
concomitantly, leucine-rich repeat ig domain-containing nogo-interacting protein  <dig>  known to mediate the collapse of growth cones in the presence of certain myelin proteins  <cit>  is also upregulated in our data.
rtnr <dig> and lingo <dig> form a complex with the p75neurotrophin receptor to mediate the inhibition of growth cones  <cit> .
although a role for this axon guidance pathway has not been implicated in ens development, reticulon <dig>  is expressed in the mammalian adult ens  <cit> .
it will be exciting to see whether further studies can confirm definitive roles for lingo <dig>  and rtnk <dig> 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  families are known to have a role in axon guidance.
a disintegrin and metalloprotease  is from the metzincins superfamily of metalloproteases  <cit> .
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  <cit> .
interestingly, some of the axonal guidance pathway genes enhanced in our dataset have been implicated in the development of the nervous system.
for example, adam <dig> is necessary in pns development and deficiency leads to hypomyelination of peripheral nerves and ataxia  <cit> .
adam <dig> deficient mice display defects in the proliferation and differentiation of glia  <cit> .
the neurotrophic tyrosine receptor kinase  family comprises of receptors that are required to maintain synaptic strength and plasticity in the nervous system  <cit> .
the three genes associated with axonal guidance from this family are ntrk <dig>  ntrk <dig> and ntrk <dig>  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 .
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 , neuromedinu , cocaine and amphetamine-regulated transcript <dig> , secretograinin ii  and neuromedin b  are neuropeptide precursors, neuropeptides, or members of neuropeptide families, expressed during and required for ens development and function .
sectretoneurin, a functional neuropeptide expressed in the ens, is derived from the secretograinin ii gene  <cit> .
the cart genes have been associated with appetite and energy control in zebrafish  <cit> , and with energy balance in chicken  <cit> .
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  <cit> .
in their study heanue and pachnis also showed that cart positive neurons represented a subset of enteric neurons  <cit> .
duplication events in the zebrafish genome have resulted in  <dig> cart genes compared to a single gene, cart propeptide , in mouse  <cit> .
we find that in our dataset only cart <dig> seems to be enriched .
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  <cit> .
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 .
syntaxin  <dig> , syntaxin binding protein  <dig>  and snap- <dig> are all components of the snare  complex of presynaptic proteins that facilitate the fusion of synaptic vesicles with the plasma membrane during the release of neurotransmitters  <cit> .
inhibitors of snare proteins have been shown to decrease enteric neural crest migration and ens precursor neurite extension  <cit> .
complexin ii , synapsin ii , and synaptotagmin  <dig>  modulate vesicle fusion and neurotransmitter release.
syt <dig> is a synaptic vesicle transmembrane protein involved in ca2+-sensitive triggered neurotransmitter release expressed in enteric neurons  <cit> .
syn <dig> a neuron-specific phosphoprotein, assists in clustering and recycling of synaptic vesicles  <cit> .
cplx <dig> is a soluble pre-synaptic protein that binds the snare complex aids in the fusion of synaptic vesicle but negatively regulates vesicle fusion  <cit> .
in mammals the two paralogous genes stx1a and stx1b are often expressed in the same cells, and function often redundantly  <cit> .
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  <cit> .
snap- <dig> is also expressed in the ens of rodents and humans and localizes to the presynaptic membrane of cholinergic and nitrergic enteric neurons .
during vertebrate evolution, the teleost fish lineage experienced an additional genome duplication  <cit> .
as with other genes this led to duplications and retention of additional stxbp <dig> and the snap- <dig> genes in the zebrafish genome.
thus, the zebrafish genome contains two stxbp <dig> paralogues, stxbp1a and 1b and two paralogues of snap- <dig>  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  <dig>  and  <dig>  with an average log2fold change expression of  <dig>   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 , so we examined the gene expression levels of phox2b and its paralogue phox2a.
both phox2a and phox2b were expressed at log2foldchange higher than  <dig>  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 .
phox2a is also similar to phox2b in its expression in the enteric neurons, cranial glia and other differentiating autonomic neurons  but its expression appears to be more restricted to a subset of enteric neural crest  <cit> .
both phox2b and phox2a have been shown to transcriptional regulate tlx <dig>  <cit> .
loss of phox2b gives rise to aganglionic gastrointestinal tract  <cit>  while in the absence of phox2a some enteric neurons are present in the oesophogeal region  <cit> .
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  <dig> , islet2a , and chromodomain helicase dna binding protein  <dig> .
interestingly, tlx <dig>  has been shown to be a target for phox2b  <cit>  and phox2a  <cit> .
tlx <dig> deficiency leads to megacolon and myenteric neuronal hyperplasia in mice , which may explain its association with intestinal motility disorders like intestinal neuronal dysplasia .
isl2a, a lim homeobox transcription factor, plays a critical role in lineage specific differentiation of motoneurons in the spinal cord  <cit> .
chd <dig> encodes a chromatin remodeling protein that is specific to neurons  <cit>  and promotes neurogenesis by activating the transcription of neural differentiation genes and inhibiting genes for non-neural fates .
isl2a and chd <dig> 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  <dig> of larval development, some proportion of neural crest cells are differentiated and forming a functioning neuronal network within the ens.
consistently we found, elavl <dig> and elavl <dig> two pan-neuronal rna binding proteins expressed in postmitotic neurons were equally highly enriched with log2foldchanges of  <dig>  and  <dig>  respectively.
this significant enrichment is consistent with increasing number of differentiated enteric neurons  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  <dig> , consistent with its role in ens development .
fibroblast growth factor  <dig>  is a members of the fibroblast growth factor  family  <cit> .
although fgf <dig> 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.
fgf <dig> is highly expressed in the embryonic cns and has been shown to regulate neural development by stabilizing microtubules .
mice deficient in fgf <dig> have defects in axonal branching, learning and memory.
in addition to the cns, fgf <dig> expression is also present in the enteric ganglia of mice  <cit> .
in zebrafish there are two fgf <dig> genes, fibroblast growth factors 13a and 13b   <cit> .
fgf13b is expressed in the developing zebrafish ens, and although not yet examined in the ens, fgf13a is expressed in neurons of the cns  <cit> .
our data demonstrates that both fgf13a and fgf13b are expressed in transcriptome of the developing ens, suggesting fgf <dig> 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 chimerin <dig> .
chn <dig> a gtpase-activating protein expressed in neurons in embryogenesis during a time of differentiation and synaptogenesis  <cit> .
to our knowledge, the expression of, or a functional role for chn <dig> in ens development has not been investigated.
intriguingly, a search of mouse genome informatics  database and the international mouse phenotyping consortium  identified a transgenic reporter line, chn1tm <dig> vlcg, driving lacz expression under the chn <dig> promoter.
lacz positive cells can be seen in the gastrointestinal tract in locations consistent with the position of the myenteric plexus  suggesting a potentially conserved role of chn <dig> in vertebrate ens development.
dickkopf  <dig>  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 .
dkk <dig> inhibits wnt causing the disassembly of synapses in mature neurons  <cit> .
it may be that an enrichment of dkk <dig> expression in the enteric neuronal population data reflects a role for dkk <dig> in synapse disassembly, an essential part of the wiring and refinement process within the developing ens neural circuitry   <cit> .
the gfp-negative non-neuronal cell population of the intestine
although the goal in this study was to define the ens  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 , 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  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 b <dig>  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  channel was significantly enriched.
cftr is expressed in the intestinal epithelium and is required to maintain water homeostasis in the intestine  <cit> .
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 .
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  .
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.
CONCLUSIONS
we report a comprehensive gene expression profile of a population of enteric neural crest progenitors and their enteric nervous system derivatives during normal development.
we also report the gene expression profile of the cells that constitute the microenvironment that the enteric crest and derivative reside in, the intestinal tract.
our results confirm previously reported genes and pathways known to be required for ens development and function, as well as novel genes and pathways, suggesting that this dataset provides valuable new insights into the genetic, cellular and molecular mechanisms driving the development and maintenance of a functioning ens.
additional files
additional file 1: figure s <dig>  fluorescent activated cell sorting  of the dissociated zebrafish.
fluorescent activated cell  analysis of sorted gfp-positive and gfp-negative cells.
viable cells separated into gfp-positive cells  and gfp-negative cells .
the table below gives an approximate survival and sorting rate of the sample.
this particular analysis had  <dig> % of live cells, among which  <dig> % were positive and  <dig> % were negative.
additional file 2: table s <dig>  list of differentially expressed genes between the two cell populations of the intestine.
the list in the first sheet represents all the genes filtered based on their adjusted p-value  while the second sheet shows the genes further filtered based on log2foldchange .
sheets  <dig> and  <dig> show the separated gene lists from the second sheet into upregulated genes in the neuronal cells and downregulated genes in the neuronal cell population, respectively.
the scatter plot and the heat map in the main paper  were generated using the data provided in this table.
additional file 3: figure s <dig>  the complete hierarchy showing the biological processes category of go enrichment in neurons.
go enrichment analysis of biological processes was performed with  <dig> genes which resulted in  <dig> enriched go terms  and  <dig> edges.
the size of the nodes refers to the number of genes associated with the go term while the color represents the significance of association with the go term.
the low significant nodes to the higher significance ranges from white to dark orange nodes.
additional file 4: figure s <dig>  top  <dig> significant go terms from the biological process category in the neurons.
a bar graph illustrating the first significant  <dig> biological processes go term nodes based on the p-value.
the y-axis shows the top  <dig> go terms for the biological processes, while the x-axis shows the corresponding –log.
the specific values can be found in the additional file  <dig>
additional file 5: figure s <dig>  top  <dig> go categories based on the highest percentages of genes associated with the go terms.
a bar graph illustrating the first  <dig> nodes based on the percentage of genes associated with a specific go term.
the y-axis shows the top  <dig> go terms for the biological processes, while the x-axis shows the corresponding percentage of genes.
the specific values can be found in the additional file  <dig>
additional file 6: figure s <dig>  the complete hierarchy showing the biological processes category of go enrichment in non-neurons.
a go enrichment analysis of biological processes in the non-neuronal genes.
the analysis was performed with  <dig> genes which resulted in  <dig> enriched go terms  and  <dig> edges.
the size of the nodes refers to the number of genes associated with the go term while the color represents the significance of association with the go term.
the low significant nodes to the higher significance ranges from white to dark orange nodes.
additional file 7: figure s <dig>  top  <dig> significant go terms from the biological process category in the non-neurons.
a bar graph plotted with the top  <dig> significant go terms nodes based on the p-value.
the y-axis shows the top  <dig> go terms for the biological processes, while the x-axis shows the corresponding –log.
the specific values can be found in the additional file  <dig>
additional file 8: figure s <dig>  top  <dig> go categories with the higher percentage of genes associated with a specific go term in the non-neurons.
a bar graph plotted with the first  <dig> nodes representing the highest percentage of genes associated with a specific go term.
the y-axis shows the top  <dig> go terms for the biological processes, while the x-axis shows the corresponding percentage of genes.
the specific values can be found in the additional file  <dig>
additional file 9: figure s <dig>  the complete hierarchy showing the molecular function category of go enrichment in neurons.
a go enrichment analysis of specific molecular function the neuronal genes.
the analysis was performed with  <dig> genes which resulted in  <dig> enriched go terms  and  <dig> edges.
the size of the nodes refers to the number of genes associated with the go term while the color represents the significance of association with the go term.
the low significant nodes to the higher significance ranges from white to dark orange nodes.
additional file 10: figure s <dig>  top  <dig> significant go terms from the molecular function category in neurons.
a bar graph was plotted with the first  <dig> significant go term nodes based on the p-value.
the y-axis shows the top  <dig> go terms for the molecular function, while the x-axis shows the corresponding –log.
the specific values can be found in the additional file  <dig>
additional file 11: figure s <dig>   top  <dig> go categories with the highest percentage of neuron associated genes.
a bar graph plotted with the first  <dig> nodes based on the percentage of genes associated with a particular go term.
the y-axis shows the top  <dig> go terms for the molecular function, while the x-axis shows the corresponding percentage of genes.
the specific values can be found in the additional file  <dig>
additional file 12: figure s <dig>  the complete hierarchy showing the cellular component category of go enrichment in neurons.
a graph of the cellular component go term analysis of the neuronal genes.
the analysis was performed with  <dig> genes which resulted in  <dig> enriched go terms  and  <dig> edges.
the size of the nodes refers to the number of genes associated with the go term while the color represents the significance of association with the go term.
the low significant nodes to the higher significance ranges from white to dark orange nodes.
additional file 13: figure s <dig>  top  <dig> significant go terms from the cellular component category in neurons.
a bar graph with the top  <dig> most significant nodes based on the p-value.
the y-axis shows the top  <dig> go terms for the cellular component, while the x-axis shows the corresponding –log.
the specific values can be found in the additional file  <dig>
additional file 14: figure s <dig>  top go categories with the largest number of gene association in neurons.
a bar graph was plotted with the top nodes based on the percentage of genes associated with a particular go term.
the y-axis shows the top  <dig> go terms for the molecular function, while the x-axis shows the corresponding percentage of genes.
the small percentage of gene associated with go terms beyond the 41st category were too low to be plotted.
the specific values for the  <dig> nodes can be found in the additional file  <dig>
additional file 15: figure s <dig>  the complete hierarchy showing the molecular function category of go enrichment in non-neurons.
a go enrichment analysis of molecular function in the non-neuronal genes.
the analysis was performed with  <dig> genes which resulted in  <dig> enriched go terms  and  <dig> edges.
the size of the nodes refers to the number of genes associated with the go term while the color represents the significance of association with the go term.
the low significant nodes to the higher significance ranges from white to dark orange nodes.
additional file 16: figure s <dig>  top  <dig> significant go terms from the molecular function category in the non-neurons.
a bar graph was plotted with the top  <dig> significant go term nodes based on the p-value.
the y-axis shows the top  <dig> go terms for the molecular function, while the x-axis shows the corresponding –log.
the specific values can be found in the additional file  <dig>
additional file 17: figure s <dig>  top  <dig> go categories with the largest number of gene associated with go terms in the non-neurons.
a bar graph plotted with the top  <dig> nodes representing the highest percentage of genes associated with a specific go term based on the percentage of genes associated.
the y-axis shows the top  <dig> go terms for the molecular function, while the x-axis shows the corresponding percentage of genes.
the specific values can be found in the additional file  <dig>
additional file 18: figure s <dig>  the complete hierarchy showing the cellular component category of go enrichment in non-neurons.
a go enrichment analysis of cellular function in the non-neuronal genes.
the analysis was performed with  <dig> genes which resulted in  <dig> enriched go terms  and  <dig> edges.
the size of the nodes refers to the number of genes associated with the go term while the color represents the significance of association with the go term.
the low significant nodes to the higher significance ranges from white to dark orange nodes.
additional file 19: figure s <dig>  top significant go terms from the cellular component category in non-neurons.
to identify and visualize some of the top significant go terms, a bar graph was plotted with the top nodes based on the p-value.
the y-axis shows the top  <dig> go terms for the cellular component, while the x-axis shows the corresponding –log.
the nodes with significant p-values beyond the 25th category were too low to be plotted.
the specific values for the  <dig> nodes can be found in the additional file  <dig>
additional file 20: figure s <dig>  top go categories with the highest number of gene association in non-neurons.
a bar graph was plotted with the top nodes based on the highest number of genes associated with a specific go term.
the y-axis shows the top  <dig> go terms for the molecular function, while the x-axis shows the corresponding percentage of genes.
the nodes with significant p-values beyond the 25th category were too low to be plotted.
the specific values for the  <dig> nodes can be found in the additional file  <dig>
additional file 21: table s <dig>  list of go terms and their corresponding values.
table showing the top go terms in the three categories; biological processes, molecular function and cellular components in both the cell populations  with their corresponding p-value as well as the percentage of genes associated from the data set.
the values have been plotted and the bar graphs can be seen in the additional files  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>
additional file 22: table s <dig>  list of molecular functions and their associated values for the present dataset.
table of all the values of the molecular functional enrichment analysis performed on the dataset.
the analysis was accomplished with  <dig> genes.
the table gives the number of “molecules”  associated with the particular molecular function as well as the percentage of genes that was used to create the pie charts  in both the intestinal cell populations.
additional file 23: table s <dig>  list of the canonical pathways in the two intestinal cell populations.
pathway analysis on the genes enriched in both the populations to identify the most significant pathways.
the analysis was accomplished using  <dig> genes.
the table provides the top canonical pathways with their respective –log and the genes from the dataset that are known to be associated with that pathway.
abbreviations
bppase pair
d.p.fdays post, fertilization
degdifferentially expressed genes
ensenteric nervous system
facsfluorescence, activated cell sorting
fdrfalse discovery rate
fpkmfragments per kilobase of transcript per million mapped reads
gfpgreen fluorescent protein
gfp + vegreen fluorescent protein positive
gfp-vegreen fluorescent protein negative
gigastrointestinal
gogene ontology
gtpgene transfer format
gwasgenome wide association studies
h.p.fhours post fertilization
hscrhirschprungs
ipaingenuity pathway analysis
nccneural crest cells
qpcrquantitative polymerase chain reaction

 the enteric nervous system, the largest division of the peripheral nervous system, is composed of a network of neurons and glia that innervate the gastrointestinal tract. the ens functions within the gastrointestinal tract to control motility for mixing and moving food, absorption of nutrients, epithelial secretions, and blood circulation. it arises from a population of multipotent neural crest cells, which migrate from the vagal, and in some species sacral, axial area to the intestine. there they continue to proliferate and differentiate into enteric neurons and glia within the gastrointestinal wall. the migration, specification and differentiation of the enteric neural crest are determined by intrinsic and extracellular cues. thus, the enteric neural crest cells interact with other cells and respond to the local environments within the gut. defects in the migration, proliferation or differentiation of the neural crest progenitors have been shown to lead to enteric neuropathies and are implicated in a variety of gastrointestinal dysmotility disorders. previous studies have identified many important genes and pathways required for normal ens development. best understood are genes whose role in ens development have been elucidated by loss-of-function studies or the study of congenital birth defects. these genes encode a variety of transcription factors, cell-cell signaling molecules, neuropeptides, biosynthesis enzymes, cell adhesion molecules, soluble growth factors, extracellular matrix and cytoskeletal proteins. categories of genes functionally linked to ens development have been expertly collated and described in a number of recent reviews. among the most studied genes in ens development, are the genes required for the ret signaling pathway. mutations in the ret tyrosine kinase receptor gene, it’s co-receptor gdnf family receptor alpha, and it’s ligand gdnf emerged as the most common causes of hirschprung disease. in hscr patients, insufficient proliferation and migration of enteric neural crest lead to partial loss of the ens generally characterized by a complete lack of enteric neurons within the caudal region of the gastrointestinal tract. the ret signaling pathway is most commonly affected in hscr, accounting for 20–50% of familial and 15–30% of sporadic cases. the etiology of hscr has recently been extensively covered in several reviews. in addition to ret, more than a dozen other genes including dna -methyltransferase beta, endothelin converting enzyme, endothelin, endothelin receptor type b, glial cell derived neurotrophic factor, gfrα kinesin family member 1b, l cell adhesion molecule, neurotrophin, neuregulin, neuregulin, neurturin, neurotrophic receptor tyrosine kinase, phox2b, prokineticin, prokineticin receptor, prokineticin receptor, persephin, semaphorin 3a, semaphorin 3d, sry-box, transcription factor, and zinc finger e-box binding homeobox are known to be associated with, or cause hscr. yet, the exact genetic cause of 50% of hscr cases remain unknown. as the enteric neural progenitors migrate into the gut, they proliferate and differentiate into at least different neuronal subtypes that differ in their function, physiology, morphology, and neurotransmitter expression. in mammals, these subtypes form ganglia that are organized into a network of two interconnecting plexi. the myenteric plexus is embedded between the circular and longitudinal smooth muscle layers, and the submucosal plexus is situated between the muscle and the mucosal layer. each neuronal subtype expresses different combinations of neuropeptides giving them a unique neurochemical code. along the length of the intestinal tract, the combinations of neuronal subtypes vary between the different regions of the gut. although the diversity of enteric neuronal subtypes and their neurochemical code has been well described, our understanding of the genes required for enteric neuronal specification, differentiation and for establishing neural connectivity within the ens remains rather limited. to fill this gap we pursued a broad strategy to identify candidate genes required to drive neural crest precursors to differentiate into enteric neurons including potential signaling molecules for the formation of neuronal subtypes and their integration into a functional neural network in zebrafish. zebrafish embryos and larvae have emerged in recent years as a powerful complementary vertebrate model for ens development. although some important anatomical differences to the ens in mammals exist e.g. the neurons and glia of the zebrafish myenteric plexus do not form ganglia and the submucosal plexus is absent, the signaling pathways and genes associated with ens development remain conserved. markers of differentiated enteric neurons can first be detected in the developing intestine approximately 48 hours post fertilization. over the next 24 hours, neural crest cells continue to migrate and proliferate and by hpf the entire length of the gut is colonized. the first contractions within the intestine are seen at 4 days post fertilization and by dpf the larval gut is functional and the larvae are capable of feeding. in this study, 7 days post fertilization larvae from a transgenic zebrafish line were used to isolate enteric neuronal cells from the intestine. 7dpf represented the time in development at which a subset of neurons has fully differentiated; yet the ens is not fully developed and likely to retain some precursors. to track and distinguish precursor cells and differentiating enteric neurons we made use of a zebrafish transgenic line in which egfp is expressed under the regulatory/promoter region of the transcription factor phox2b, tg w. this strategy was chosen because phox2b is required for the specification and development of the neural crest-derived autonomic nervous system including the enteric neurons and glia in various model systems. loss of phox2b during development leads to loss of enteric ganglia in mice while in heterozygous animals a normal gi tract is observed. in humans, mutations in phox2b are linked to a loss of enteric neurons in the gastrointestinal tract leading to hscr, and many cases have been associated with congenital central hypoventilation syndrome. determination of the genetic processes underlying the syndrome revealed that mutations of the phox2b gene are principally responsible for the broad range of symptoms encountered by disrupting early development of autonomic neurons. importantly, phox2b has been shown to regulate ret expression and the genetic interaction between ret and phox2b has been demonstrated to be critical for normal ens development. in mice phox2b expression is turned on as the neural crest cells leave the neural tube, remains expressed as neural crest migrate into the intestine and differentiate into neurons and glia and continues to be expressed in ens derivatives into adulthood. our goal was to generate a comprehensive transcriptional profile of enteric neurons along the entire intestinal tract during normal development. to improve the chances of obtaining discrete transcript counts and unbiased detection of novel or low-abundance transcripts across a broader dynamic range, endogenous enteric neurons and enteric neural crest progenitors were transcriptionally profiled using rna-seq. complimentary to profiling the enteric neurons, we also utilized rna-seq to transcriptionally profile their surrounding non- neural microenvironment, the intestine including the mucosa, musculature and associated vasculature. thus, this strategy enabled us to perform a comparative gene expression study between the neuronal cell and the non- neuronal cell population of the intestine. our analyses identified scores of candidate genes that may function as cell fate determinants during ens development, may constitute useful markers to distinguish different ens subtypes, and define the surrounding cellular microenvironment of the ens on a molecular level. a subset of the genes has been validated by expression studies, and selected for ongoing functional studies. the result of these studies will not only expand our knowledge of genes required for normal ens development and a better understanding of the cellular and molecular basis of the complex process of development of the ens but may help inform stem cell therapies aimed at directing lineage specific enteric neuron development. methods animal husbandry and line the embryos and larvae were obtained from a transgenic zebrafish line as described previously and provided by dr. a nechiporuk, in which green fluorescent protein is expressed under the regulatory elements of the transcription factor phox2b tg w. the embryos and larvae were raised at  °c in a light/dark cycle according to standard zebrafish husbandry protocols and staged to seven days post fertilization. larvae were prescreened for gfp expression and anesthetized prior to dissection. all animal experiments were performed in accordance with the u.s. laws, guidelines and policies of laboratory animal care and use and approved by the institutional animal care and use committee and institutional biosafety committee at iowa state university. tissue collection and rna isolation to acquire the enteric neuronal population, larvae from the transgenic line tg w were dissected at 7 days post fertilization to isolate the entire gut. guts were dissociated into cells using papain digestion. for dissociation, 20 μl of mgp/ml papain was added to the dissected guts along with 360 μl of a solution of hbss, % of hepes and 20 μl of an activation agent according to the manufacturer’s instructions. the samples were incubated at 37 °c for a total of 30 minute with trituration of the sample every 10 minutes. a final incubation was done for 13 minute with 20 μl of dnase i. the samples were resuspended in 500 μl of % bsa in pbs. fluorescent activated cell sorting was done on each sample using a bd facsaria iii available at the iowa state flow cytometry facility. gfp-positive cells were detected through a 525/50 nm bandpass filter, and simultaneously checked for viability, and sorted directly into lysis buffer. both, the gfp positive neuronal cells and the gfp negative non-neuronal cells were collected directly into commercially available lysis buffer and stored at -80 °c until separated cell populations from multiple facs sorting experiments could be pooled and processed. rna isolation was done using pooled samples of – cells and the ambion™ ‘s rnaqueous®- micro total rna isolation kit. the protocol was followed as per kit instructions with the following modifications. a total volume of 20 μl of elution buffer was used and the column was incubated at room temperature for 10 minutes before the final spin to collect the eluted rna. a total of larval intestines were dissected to collect gfp positive cells to obtain sufficient rna for two biological replicates of the gfp-positive cell population. gfp negative cells were collected from the same larval intestine samples for the two biological replicates of the gfp-negative cell population. the rna quality check was checked using an agilent bioanalyzer. rna integrity number s of - were measured for all rna samples. libraries were prepared using illumina’s truseqv kit. sequencing and raw read processing the sequencing was done using an illumina hiseq sequencing system at the dna facility of the iowa state university office of biotechnology. deep sequencing of the mrna was performed with 100 bp paired-end reads. each biological replicate was considered as a separate sample. preprocessing of the raw reads was done using trimmomatic, to remove adapters and the low quality reads. the preprocessing was done using the parameters leading: trailing: slidingwindow:4: minlen: the removal of the bases was based on the phred score where the threshold was a window size of was used which removed the bases with an average phred score of less than the quality of the trimmed reads was analyzed and visualized using fastx-toolkit. alignment, transcript abundance and differential gene expression analysis the trimmed reads were mapped to the zv danio rerio reference genome using tophat which uses bowtie as the mapping tools. the version zv of gtf was used for annotation of the mapped reads. the expression levels of the mapped reads were obtained using featurecounts to specifically count the transcripts mapping to exons. for the analyses shown in figs.  and 3a and b, we used cufflink to get the fragments per kilobase of transcript per million mapped reads. an r package deseq was used to perform the differential gene expression analysis between the gfp-positive enteric neuron gut population and the gfp-negative non-neuronal gut population. the software considered the biological replicates from both the cell populations during the analysis to provide an expression level for the genes. the differentially expressed genes were filtered using the false discovery rate adjusted p-value ≤ . the significantly upregulated and downregulated genes were chosen based on a log2foldchange value of. the heatmap was generated using the r function “pheatmap”. for the correlation analysis, the spearman correlation coefficient for pairwise samples was calculated using r function “cor”. the scatter plots were visualized using r function “plot”. the expression value was log-transformed using log. gene ontology enrichment and pathway analysis a gene ontology enrichment analysis was performed using the cytoscape plugin bingo with the adjusted p-value <  where the fdr correction was done by benjamin and hochberg’s method. pathway analysis was performed to identify pathways most enriched in the enriched and diminished genes in the neuronal population using qiagen’s ingenuity pathway analysis. to perform ipa pathway analysis human and mouse orthologues of the identified zebrafish ensembl genes were converted and ipa pathway analysis was performed on 44% of the genes from the complete data set. the datasets analyzed were the list of significantly enhanced genes; and the list of significantly diminished genes. the percentage of molecules was calculated for the analysis using a total of genes from the enhanced list and a total of genes from the depleted list. qpcr eight candidate genes were chosen for qpcr. qpcr was performed on the same cdna libraries used for illlumina sequencing to validate the rnaseq results. primers were designed using ncbi primer designing tool and synthesized at the iowa state dna facility. primers to chn elavl chata, phox2bb, syn2a, amy2a, fli1a, and hprt were designed to span exon-exon junction and give a single product of 150–200 bp. individual qpcr reactions were carried out with the standard steponeplus cycling protocol using perfecta sybr® green fastmix, rox according to manufactures instructions. for each qpcr run, the cycle threshold were exported into microsoft excel for statistical analysis. we used 2-ΔΔct method of relative quantification to estimate copy number of selected genes. the expression levels were normalized to zebrafish hypoxanthine phosphoribosyltransferase as a reference gene.  transcriptional profiling of ens cells and the zebrafish intestine to identify genes required for the specification and differentiation of enteric neurons from neural crest progenitors we isolated and carried out rna-seq transcriptional profiling on a cell lineage expressing the transcription factor phox2b. phox2b has been shown to be expressed in early enteric neural crest progenitors as they migrate into the gut. as development proceeds phox2b expression remains high in a subset of the differentiated neurons into adulthood and decreases to lower levels in the developing glia. phox2b expression is retained in a subset of adult ens derivatives. therefore, a transgenic reporter line driving expression of gfp under the phox2b promoter, tg w seemed a good candidate line to identify enteric neuron precursors and enteric neurons as gfp-positive cells within the intestine, and to facilitate separation of the neuronal from non-neuronal cell populations of the intestine for subsequent transcriptional profiling. indeed, personal observations confirmed that gfp expression appears in the enteric neural crest as they migrate into the intestine and continues throughout development and into adulthood in the zebrafish tg w line. thus, gfp-positive cells of the intestine mark both enteric neurons and enteric neuron precursors in 7 days post fertilization larvae and may facilitate separation of the ens cells from other intestine tissues by fluorescence activated cell sorting. however, gfp expression is not restricted to ens in dpf larvae. endogenous phox2b, and therefore gfp, in the zebrafish tg w line is also expressed in cells of the hindbrain and cranial ganglia as seen in fig. 1a. to ensure that only gfp-positive cells that comprise the ens were profiled, the intestinal tracts of 7dpf zebrafish larvae were manually dissected away from the remaining larval body prior to dissociation and facs sorting.fig. image and flow chart illustrating source of gfp labeled cells, sample preparation and bioinformatics analysis. a image of a 7 days post fertilization zebrafish larvae of the transgenic line, tg. the boxed area shows the intestine with gfp-positive enteric neurons. b schematic diagram illustrating the experimental approach. intestines were dissected out to avoid contamination with gfp-positive cells expressing phox2b in the hindbrain or the spinal cord. the intestines were dissociated using papain digestion followed by fluorescent activated cell sorting. rna isolated from the sorted cell populations was used to construct cdna libraries. raw reads from illumina sequencing were processed and analyzed using various bioinformatics analyses programs at this stage in development, the enteric neuron precursors have migrated to the end of the intestine, and a sufficient number of progenitors have differentiated into different ens subtypes to form a functioning ens. enteric neurons present at this time reflect at least some of the diversity of enteric neuron subtypes. after enzymatically treating the dissected intestines with papain, the dissociated cells were fac sorted and viable gfp positive cells were collected. cells from whole intestines of wild type zebrafish were used as a negative control to establish a threshold for gfp expression. viable gfp-negative cells representing the remaining cell types in the intestine were collected simultaneously and processed in the same manner as the gfp-positive cells. to minimize potential changes in gene expression patterns, to obtain sufficient amounts of rna to sequence, and to limit amplification steps, dissection times were limited and multiple independent cell sorts were pooled. a total of larval intestines were dissected to collect gfp positive cells to obtain sufficient rna for two biological replicates of the gfp-positive cell populations. gfp negative cells were collected from the same larval intestine samples for the two biological replicates of the gfp-negative cell populations. libraries were generated from the pooled samples using unamplified rna. biological duplicates of both the gfp-positive and gfp-negative cells were sequenced using illumina hiseq following sequencing, a developed pipeline was implemented to process the raw rna-seq reads. figure 1b outlines the basic workflow from tissue collection to data processing. differential gene expression analysis between neuronal and non-neuronal cells of the intestine to identify genes that may be required for the differentiation and assembly of enteric neuron circuitry we looked for differences in expression profiles between the sorted neuronal and non-neuronal populations of intestinal cells. for differential gene expression analysis the raw reads were trimmed using trimmomatic to remove the adapters and the low quality reads, and mapped to the annotated d.rerio reference genome. to validate that our experimental approach separating the ens/neuronal cells from non-neuronal cells by dissociation and fac sorting, and that our differential gene expression analysis yielded meaningful results, we tested for several predicted outcomes. first, gfp mrna should be entirely restricted to the two neuronal gfp-positive cell samples and absent in the two non-neuronal gfp-negative samples. indeed, high levels of gfp were found in the two neuronal, and miniscule amounts in the two non-neuronal samples. therefore, dissociation and sorting successfully separated the two intestinal cell populations leading to a > fold enrichment of gfp in the neuronal versus non-neuronal samples. second, we would expect that endogenous phox2b mrna should also be enriched in the gfp-positive cells. indeed, endogenous phox2b was highly enriched in the gfp-positive cells, a ~ fold enrichment compared to the gfp-negative samples while high and low expressing house-keeping genes like actb, gapdh, and hypoxanthine phosphoribosyltransferase showed similar high and low expression in both sorted cell populations. this demonstrates that gfp expression in the intestine in the zebrafish tg w line is indeed labeling the phox2b-expressing ens cell population.fig. experimental validation of reliable separation of the neuronal ens from the non-neuronal intestinal cell population. a the bar graph and corresponding data in table show the expression levels in fpkm of five selected genes in the two neuronal replicate and the two non-neuronal replicate cell populations. gfp transcript levels show the level of gfp transcript driven from the phox2b transgene. the expression values of phox2bb measure the endogenous transcript levels of the phox2b gene. expression levels for the three house-keeping genes actb gapdh, and hprt are also shown. qpcr validation of expression levels of genes comparing the gfp-positive neuronal populations versus the gfp-negative non-neuronal. chn elavl chata, phox2bb, and syn2a were more highly expressed in the neuronal population, whereas amy2a and fli1a were more highly expressed in the non-neuronal population. changes in gene expression were calculated using 2− ΔΔct method with hprt as the house- keeping gene. the mean relative expression of the neuronal population was significantly different from that of the non-neuronal population to verify that our estimated abundance levels for transcripts by rna-seq are consistent with the actual expression levels in our samples we determined the relative abundance levels of a few candidate genes using semi-quantitative pcr. based on our experimental strategy we would predict that candidate genes that are known to control ens development and differentiation are entirely expressed or strongly enriched in the gfp-positive cells. conversely, we would expect that that known genes involved in patterning and differentiation of the non-neuronal intestine should be strongly enriched in the gfp-negative samples. indeed, relative expression levels of the candidate genes using semi-quantitative pcr correlated with the rnaseq results, and showed the predicted differences in both cell populations. elav like neuron-specific rna binding protein, a rna binding protein found in post mitotic neurons and specific only to neurons, choline o-acetyltransferase a, whose gene product is required for biosynthesis of the neurotransmitter acetylcholine, phox2bb, a dna-binding protein which is associated with the development of noradrenergic neurons, chimerin, a gene that is primarily expressed in the neurons and is involved in neuronal signal-transduction and synapsin 2a, a synaptic gene that has been involved in neurotransmitter release and synaptogenesis, were all found to be highly expressed in the neuronal population consistent with their roles in neurogenesis and synaptogenesis. whereas amylase alpha 2a, one of the genes that is highly involved in proper digestion by helping in breakdown of polysaccharides and oligosaccharides and fli- proto-oncogene, ets transcription factor a, showed higher levels in the non-neuronal cell population. hrpt was chosen as a more reliable house keeping gene. our results showed that there was a high correlation of relative expression levels between the rna-seq results and the qpcr for all eight genes. furthermore, the expression for both groups of marker genes are almost entirely separated demonstrating that the ens, and the non-neuronal cell populations have been efficiently separated by dissociation and sorting with a minimum of cross contamination. all of the above indicates that a genome scale differential expression analysis between gfp-positive and gfp-negative cell populations will yield meaningful transcriptional profiles of neuronal/ens cells and non-neuronal cells of the zebrafish intestine, respectively. to do so we compared the genes expressed in the neuronal and non-neuronal cell populations of the intestine. differentially expressed genes were identified by benjamin and hochberg’s test ≤  and log2foldchange ≥  ). of these deg’s, genes were found to be significantly enriched in the neuronal cell population. for ease in this manuscript, we refer to genes expressed at significantly higher levels in the neuronal population in comparison to the non-neuronal population as ‘enriched’ and to genes most significantly lower in the neuronal population as ‘depleted’. both the neuronal and non-neuronal cell populations showed a high level of similarity in the fpkms between duplicates prompting us to expand our analysis to compare all the genes expressed above fpkm of in our data sets. to analyze how reproducible the estimated gene expression levels between the biological duplicates were, we calculated the spearman correlation coefficient for all expressed genes. a comparison of the biological replicates of the neuronal and non-neuronal cells yielded average coefficients of and respectively demonstrating a high correlation of gene expression levels between the duplicates. the high level of concordance between biological replicates confirms the quality and consistency of our samples and dataset.fig. differential gene expression analysis of 7dpf danio rerio enteric neurons. a and b. plot showing the correlation between the two biological replicates of the gfp-positive neuronal population and of the gfp-negative non-neuronal population. the x-axis represents the 1st biological replicate while the y-axis represents the 2nd biological replicate. the correlation coefficient was calculated using fpkm. the correlation coefficient for the neuronal and non-neuronal cell population is and, respectively. c dispersion plot showing the log2foldchange of individual reads as compared to the mean normalized counts using deseq with the adjusted p-value ≤ . the red dots represent the upregulated/enriched and the blue dots represent downregulated/depleted genes in the neuronal cells of the zebrafish intestine. d. heat map of the normalized counts for the degs between neuronal and non-neuronal cell populations. the thresholds selected are fdr ≤  and log2foldchange ≥ . the first two columns represent biological replicates of the non-neuronal cell population and the second two columns represent biological replicates of the neuronal population. each row denotes a single gene and their expression pattern across the two different samples. the blue regions in the neuronal population represent the % of the upregulated genes and the red region represents the downregulated genes. the dendogram on the x-axis shows the grouping of same samples while the y-axis shows clustering of similar expression patterns to illustrate the differentially expressed genes and their expression pattern visually, we generated a dispersion plot heat map  showing the clustering of differentially expressed gene expression profiles for the neuronal and non-neuronal samples. the heat map shows the relative expression of the genes based on their mapped read counts. enriched genes include genes such as phox2b, elavl ret and gng previously associated with ens development. interestingly, hscr related genes such as gfra, ret, phox2b, and zeb that cause aganglionosis when knocked down in zebrafish or mutated in mice are also found enriched. other elevated genes in the neuronal cell population are known to impact the ens functionally without aganglionosis including ascl cadherin egf lag seven-pass g-type receptor, frizzled class receptor, gfra hand hoxb, kinesin family members ), neural adhesion molecules, ntrk, solute carrier family member, t-cell leukemia homeobox. among the genes significantly depleted in the neuronal cell population, we find genes that lead to intestinal aganglionosis endothelin converting enzyme, indian hedgehog homolog a, integrin, beta 1b., itgb. and genes that cause alterations in ens function like distal-less homeobox 2a, erb-b receptor tyrosine kinase 3a, noggin, sonic hedgehog b, and zic family member. therefore, through the de analysis we were able to identify previously characterized as well as many novel genes that are specifically expressed in the neuronal versus non-neuronal cell population verifying the quality of our samples and datasets. these findings suggest that the identified enriched and depleted genes constitute a comprehensive and valuable resource for discovery and subsequent targeted functional analysis of enteric neuronal and non-neuronal cell types within the intestine, respectively. go and molecular enrichment analysis within the two intestinal cell populations to identify specific functional, molecular as well as structural gene categories within the enriched or depleted gene sets, a gene ontology enrichment analysis was performed. among genes associated with biological processes, go terms were enriched within the gfp-positive, neuronal cell population, while go terms were enriched among genes within the gfp-negative, non-neuronal cell population. we generated comprehensive go term hierarchies where each node represents the number of genes associated with a particular go term, and the intensity of the orange to yellow to white color represents the significance from high to low, respectively. biological process hierarchies were generated for both neuronal and non-neuronal populations. similar analyses were performed for molecular function and cellular component go categories in both the neuronal and non-neuronal populations, additional files 12–14: figures s11-s, additional files 15–17: figures s14-s, and additional files 18–20: figure s17-s. some of the most significant go term nodes were selected to generate a core hierarchy for easier visualization. terms found in both sets of analyses were included to illustrate the differences in gene categories between the neuronal and non-neuronal population. interestingly the results show several of the go terms and their associated genes have reciprocal expression patterns between the neuronal and non-neuronal populations. for example, the go term ‘nervous system development’ is highly significant among genes enriched in the neuronal cell population with a p-value of e- while showing little or no significance in the depleted list. conversely, terms like ‘organ development’ and ‘regulation of primary metabolic processes’ had low or no significance in the neuronal population but were highly significantly enriched in the non-neuronal set. the list for the selected go terms present in the highlighted pink and green boxes as well as the terms specific to either population along with their p-values and percentage of associated genes corresponding to fig.  can be found in table  additional file 21: table s gives the p-value and percentage of genes for the additional top go term categories in the analyses.fig. gene ontology enrichment analysis of the differentially expressed genes in the neuronal cells. go enrichment analysis was performed using the enriched and the depleted genes from the neuronal cells. figures illustrate a subset of go terms common to both the neuronal and non-neuronal population biological properties hierarchies and a few specific to one or the other hierarchy and. the hierarchies of the go terms are based on the number of genes and the significance for the category biological process. the circles/nodes represent the different go terms, where the size of the circle represent the number of genes associated with the term, i.e., the larger the circles, the more the number of genes associated with the go term. the color scheme represents the significance: white represents the least, and orange represents the most significant terms. the white, yellow and orange circles represent the terms common to both the neuronal and non-neuronal networks, while the red and blue circumferences of the circles represent terms specific to either the neuronal or non-neuronal data sets respectively. a core hierarchy of go terms defining biological processes in enriched genes comprising the ‘neuronal’ gene set. genes are included in this analysis. b. core hierarchy of go terms defining biological processes in the depleted genes comprising the ‘non-neuronal’ gene set. genes are included in this analysis. the green and pink boxes highlight examples of sub hierarchies within the larger hierarchies of biological process go terms common to both populations but differentially enriched between the two. sub-hierarchies of biological process go terms are significantly enriched in the neuronal population, ‘ion transport’ and ‘nervous system development’ but not in the non-neuronal population. examples of go term sub-hierarchies enriched in the non-neuronal population include regulation of ‘metabolic and biosynthetic processes’ and ‘organ development’. the corresponding p-values and percentage of genes associated with the go terms in the green and pink boxes as well as the specific terms in either population are shown in table 1 the ‘--‘ indicates n/a in the table the terms specific to a single category have' --' in both the, p-value and percentage of genes columns in addition to shared go terms found in both the enriched and depleted gene sets, several significant go terms were identified belonging specifically to one set or the other. in the neuronal population the go terms ‘neurogenesis’, ‘neuron development’, and ‘axonogenesis’ were highly significant as expected, but these go terms were not present in the non-neuronal population. similarly, go terms such as ‘proteolysis’, ‘vasculature development’ and ‘blood vessel morphogenesis’ were significantly enriched in the non-neuronal but absent in the neuronal cells. in addition, the most significantly enriched go terms of the neuronal population emphasized the differentiation states and the establishment of neural circuitry between enteric neurons. terms significantly enriched include terms relating to ‘synaptic function’ and ‘neurotransmission’, and are consistent with the onset of a functioning ens at this stage in development. to gain further insights into the diverse functional molecular categories between the two cell populations, a differentially expressed gene analysis was performed using qiagen’s ingenuity pathway analysis software. the identified functional categories and their percentage based on the total number of genes in each population are shown in the pie-charts in fig.  this analysis has important limitations. the first pie-chart shows the categories of enriched genes according to their molecular function and their corresponding percentages within the neuronal population. here, the total number of identified genes was of which were mapped to the ipa knowledge base, and were subsequently used for the analysis shown. the second pie-chart shows a similar analysis in the depleted genes, with an initial input of genes, of which were mapped and were found to be analysis ready. therefore, of the gene orthologues included in the analysis genes were assigned to functional designated categories. the remaining mapped genes, though equally interesting, could not be formally assigned to a functional category due to limitations in ipa category designations and were classified as ‘others’. when the molecular categories within the enhanced and the depleted genes were examined separately, the neuronal population showed 8% ‘ion-channels’ as compared to 1% in the non-neuronal population. in contrast, there were 5% ‘peptidases’ in the non-neuronal population as compared to 2% in the neuronal population. additionally, other categories such as ‘growth factors’, ‘transcriptional regulators’ and ‘cytokines’ were enriched in the non-neuronal population as compared to the neuronal population. thus, this analysis categorized genes in both cell populations based on their association with biological or molecular functions indicating significant shifts in several categories and identifying associated genes involved in specification and differentiation of cell types among both cell populations.fig. gene function category enrichment of the significantly differentially expressed genes in the neuronal and non-neuronal gene populations: pie-chart showing the most functionally enriched terms in the significantly differentially expressed genes ipa analysis. a gene function categories enriched in the neuronal population. b gene category functions enriched in the non-neuronal population. different colors represent different molecular functions and the size of the pie slice, the percentage of genes representing this category. each category represents differentially expressed genes grouped by functional category. for example, ion-channels represent 8% of the deg in the neuronal population but only % of the genes in the non-neuronal population. the surprisingly high percentage of the genes are listed as ‘other’ because they are not represented within the functional categories defined by the ipa software parameters. number of genes associated with each category is available in additional file 22: table s3 to obtain more information on signaling pathways that can be associated with our datasets, and to gain insight into the genes comprising these pathways, additional analysis was carried out. ipa pathway analysis was used to determine the molecular pathways that were highly represented in either the neuronal or non-neuronal dataset. the final analysis was done using enriched genes and of depleted genes of the deg genes between the neuronal and non-neuronal populations. figure  shows the top enriched canonical pathways in the enhanced as well as the depleted list of genes. the top canonical pathways along with their –log and associated genes are shown in table  the complete information about the pathways can be found in additional file 23: table s as expected, ‘g-protein coupled receptor-signaling’ pathway, ‘camp-mediated signaling’, ‘nnos signaling’ pathway and ‘axon guidance’ were among the highly enriched pathways in the neuronal population. these pathway predictions implicated many additional genes as up-regulated, and, consistently, the transcript levels within our neuronal transcriptome confirmed these expectations. conversely, some of the significant pathways in the depleted list comprised of ‘il- signaling pathway’, ‘tec kinase signaling’, and ‘inhibition of matrix metalloprotease’. the first two pathways are known to regulate the immune system. interestingly, the third identified pathway involving matrix metalloproteases is known to induce cell migration. thus, the enrichment of inhibition of this pathway may suggest a mechanism for cessation of cell migration. this potential function would be consistent with this developmental stage in which a proportion of enteric neural crest cells has completed the migratory phase, has differentiated and is functionally capable of driving intestinal motility.fig. pathway analysis of the differentially expressed genes in the neuronal and non-neuronal cells of the zebrafish intestine. pathway analysis of the upregulated, downregulated, and the complete list of significantly differentially expressed genes was done using ipa. the purple bar graph represent the top canonical pathways predicted to be most likely to be activated in the neuronal genes using the genes most significantly enriched in the neuronal population while the blue bar graph represents the top canonical pathways predicted to be most likely to be activated in the neuronal genes using the genes most significantly downregulated in the non-neuronal population differentially expressed genes may drive development and differentiation of the ens in examining our transcriptomic data, we noted many genes to be highly enriched within the neuronal population compared to the non-neuronal population. these genes could be further divided in distinct categories according to their possible function as candidates for ens development and differentiation. examples of molecular categories represented in the data were dna and rna binding proteins, signaling molecules, neuropeptides, and proteins involved in synaptogenesis. interestingly but not surprising there is a significant enrichment in several neuropeptide genes such as vasoactive intestinal peptide and secretogranin ii, both of which are implicated in ens function. the cytokine gene, secretogranin ii, from which the neuropeptide secretoneurin is derived, exhibits the highest enriched fold change in expression in the neuronal population. though, secretoneurin is expressed in rat and human myenteric and submucosal ganglia and has been implicated in the regulation of motility of the gastrointestinal tract, no study involving loss of function in the ens has been reported.table potential candiate genes with their log2foldchange phox2aa chd5 isl2a elavl4a cart2 fgf13b chn1 reta dkk1b apreviously described genes in the enteric nervous system potential novel candidate genes are bolded genes such as synapsin ii, and snap required for synaptic vesicle formation and fusion and neurotransmitter release, are similarly enriched, consistent with the establishment of synapses and a functioning neural network. the transcription factor tlx required for ens development and function, is enriched, as is islet 2a, a transcription factor known to play a role in motoneuron specification in the cns. the newly identified ens related genes can be subdivided into three main categories: genes known to be expressed by neurons outside the ens but not yet associated with ens development; genes reported to be expressed in the ens but for which functional studies have not yet been reported; and novel genes that have not yet been associated with expression or function during ens development. a more detailed discussion of table  follows below. thus, our approach successfully identified genes and pathways known to be involved in neurogenesis and synaptogenesis, functions that are expected to be present in the ens cell population. furthermore, this suggests also that our dataset uncovers reliably novel genes or pathways that have not been implicated in ens development yet, and can be used as a comprehensive resource to start further investigations into their roles in ens development and differentiation. discussion 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. 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. 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 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 w 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 ectopic gfp and endogenous phox2b transcripts almost exclusively enriched in the gfp-positive cell population, housekeeping genes were similarly expressed in both cell populations, and rna-seq and qpcr results consistently showed a strong enrichment of candidate neuronal genes and candidate non-neuronal genes 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. 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 genes 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 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. 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, 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 nodes and 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 robo slit nadl, plxn, etc. 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 e- in the hierarchy specific to the neuronal population while has p-value lower than < e- 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 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, nadl fxd and slit genes. the semaphorin family’s role in axon guidance is well documented but recent evidence demonstrates that semaphorins play also an important role in ens development. in zebrafish the knockdown of sema3c and sema3d results in the loss of enteric neurons in a dose dependent manner. 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 reticulon receptor, also known as the nogo receptor, is upregulated in the ens population. concomitantly, leucine-rich repeat ig domain-containing nogo-interacting protein known to mediate the collapse of growth cones in the presence of certain myelin proteins is also upregulated in our data. rtnr and lingo 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, reticulon is expressed in the mammalian adult ens. it will be exciting to see whether further studies can confirm definitive roles for lingo and rtnk 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 families are known to have a role in axon guidance. a disintegrin and metalloprotease 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. interestingly, some of the axonal guidance pathway genes enhanced in our dataset have been implicated in the development of the nervous system. for example, adam is necessary in pns development and deficiency leads to hypomyelination of peripheral nerves and ataxia. adam deficient mice display defects in the proliferation and differentiation of glia. the neurotrophic tyrosine receptor kinase 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 ntrk ntrk and ntrk 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. 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, neuromedinu, cocaine and amphetamine-regulated transcript, secretograinin ii and neuromedin b are neuropeptide precursors, neuropeptides, or members of neuropeptide families, expressed during and required for ens development and function. sectretoneurin, a functional neuropeptide expressed in the ens, is derived from the secretograinin ii gene. 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. duplication events in the zebrafish genome have resulted in cart genes compared to a single gene, cart propeptide, in mouse. we find that in our dataset only cart seems to be enriched. 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. 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. syntaxin, syntaxin binding protein and snap- are all components of the snare 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, synapsin ii, and synaptotagmin  modulate vesicle fusion and neurotransmitter release. syt is a synaptic vesicle transmembrane protein involved in ca2+-sensitive triggered neurotransmitter release expressed in enteric neurons. syn a neuron-specific phosphoprotein, assists in clustering and recycling of synaptic vesicles. cplx is a soluble pre-synaptic protein that binds the snare complex aids in the fusion of synaptic vesicle but negatively regulates vesicle fusion. 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- is also expressed in the ens of rodents and humans and localizes to the presynaptic membrane of cholinergic and nitrergic enteric neurons. during vertebrate evolution, the teleost fish lineage experienced an additional genome duplication. as with other genes this led to duplications and retention of additional stxbp and the snap- genes in the zebrafish genome. thus, the zebrafish genome contains two stxbp paralogues, stxbp1a and 1b and two paralogues of snap- 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 and with an average log2fold change expression of 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, so we examined the gene expression levels of phox2b and its paralogue phox2a. both phox2a and phox2b were expressed at log2foldchange higher than 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. phox2a is also similar to phox2b in its expression in the enteric neurons, cranial glia and other differentiating autonomic neurons 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 tlx. 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, islet2a, and chromodomain helicase dna binding protein. interestingly, tlx has been shown to be a target for phox2b and phox2a. tlx deficiency leads to megacolon and myenteric neuronal hyperplasia in mice, which may explain its association with intestinal motility disorders like intestinal neuronal dysplasia. isl2a, a lim homeobox transcription factor, plays a critical role in lineage specific differentiation of motoneurons in the spinal cord. chd 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. isl2a and chd 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 of larval development, some proportion of neural crest cells are differentiated and forming a functioning neuronal network within the ens. consistently we found, elavl and elavl two pan-neuronal rna binding proteins expressed in postmitotic neurons were equally highly enriched with log2foldchanges of and respectively. this significant enrichment is consistent with increasing number of differentiated enteric neurons 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, consistent with its role in ens development. fibroblast growth factor is a members of the fibroblast growth factor family. although fgf 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. fgf is highly expressed in the embryonic cns and has been shown to regulate neural development by stabilizing microtubules. mice deficient in fgf have defects in axonal branching, learning and memory. in addition to the cns, fgf expression is also present in the enteric ganglia of mice. in zebrafish there are two fgf genes, fibroblast growth factors 13a and 13b. fgf13b is expressed in the developing zebrafish ens, and although not yet examined in the ens, fgf13a is expressed in neurons of the cns. our data demonstrates that both fgf13a and fgf13b are expressed in transcriptome of the developing ens, suggesting fgf 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 chimerin. chn 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 chn in ens development has not been investigated. intriguingly, a search of mouse genome informatics database and the international mouse phenotyping consortium identified a transgenic reporter line, chn1tm vlcg, driving lacz expression under the chn promoter. lacz positive cells can be seen in the gastrointestinal tract in locations consistent with the position of the myenteric plexus suggesting a potentially conserved role of chn in vertebrate ens development. dickkopf  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. dkk inhibits wnt causing the disassembly of synapses in mature neurons. it may be that an enrichment of dkk expression in the enteric neuronal population data reflects a role for dkk in synapse disassembly, an essential part of the wiring and refinement process within the developing ens neural circuitry. the gfp-negative non-neuronal cell population of the intestine although the goal in this study was to define the ens 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, 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 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 b 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 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. 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. 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.  we report a comprehensive gene expression profile of a population of enteric neural crest progenitors and their enteric nervous system derivatives during normal development. we also report the gene expression profile of the cells that constitute the microenvironment that the enteric crest and derivative reside in, the intestinal tract. our results confirm previously reported genes and pathways known to be required for ens development and function, as well as novel genes and pathways, suggesting that this dataset provides valuable new insights into the genetic, cellular and molecular mechanisms driving the development and maintenance of a functioning ens. additional files additional file 1: figure s fluorescent activated cell sorting of the dissociated zebrafish. fluorescent activated cell analysis of sorted gfp-positive and gfp-negative cells. viable cells separated into gfp-positive cells and gfp-negative cells. the table below gives an approximate survival and sorting rate of the sample. this particular analysis had % of live cells, among which % were positive and % were negative. additional file 2: table s list of differentially expressed genes between the two cell populations of the intestine. the list in the first sheet represents all the genes filtered based on their adjusted p-value while the second sheet shows the genes further filtered based on log2foldchange. sheets and show the separated gene lists from the second sheet into upregulated genes in the neuronal cells and downregulated genes in the neuronal cell population, respectively. the scatter plot and the heat map in the main paper were generated using the data provided in this table. additional file 3: figure s the complete hierarchy showing the biological processes category of go enrichment in neurons. go enrichment analysis of biological processes was performed with genes which resulted in enriched go terms and edges. the size of the nodes refers to the number of genes associated with the go term while the color represents the significance of association with the go term. the low significant nodes to the higher significance ranges from white to dark orange nodes. additional file 4: figure s top significant go terms from the biological process category in the neurons. a bar graph illustrating the first significant biological processes go term nodes based on the p-value. the y-axis shows the top go terms for the biological processes, while the x-axis shows the corresponding –log. the specific values can be found in the additional file additional file 5: figure s top go categories based on the highest percentages of genes associated with the go terms. a bar graph illustrating the first nodes based on the percentage of genes associated with a specific go term. the y-axis shows the top go terms for the biological processes, while the x-axis shows the corresponding percentage of genes. the specific values can be found in the additional file additional file 6: figure s the complete hierarchy showing the biological processes category of go enrichment in non-neurons. a go enrichment analysis of biological processes in the non-neuronal genes. the analysis was performed with genes which resulted in enriched go terms and edges. the size of the nodes refers to the number of genes associated with the go term while the color represents the significance of association with the go term. the low significant nodes to the higher significance ranges from white to dark orange nodes. additional file 7: figure s top significant go terms from the biological process category in the non-neurons. a bar graph plotted with the top significant go terms nodes based on the p-value. the y-axis shows the top go terms for the biological processes, while the x-axis shows the corresponding –log. the specific values can be found in the additional file additional file 8: figure s top go categories with the higher percentage of genes associated with a specific go term in the non-neurons. a bar graph plotted with the first nodes representing the highest percentage of genes associated with a specific go term. the y-axis shows the top go terms for the biological processes, while the x-axis shows the corresponding percentage of genes. the specific values can be found in the additional file additional file 9: figure s the complete hierarchy showing the molecular function category of go enrichment in neurons. a go enrichment analysis of specific molecular function the neuronal genes. the analysis was performed with genes which resulted in enriched go terms and edges. the size of the nodes refers to the number of genes associated with the go term while the color represents the significance of association with the go term. the low significant nodes to the higher significance ranges from white to dark orange nodes. additional file 10: figure s top significant go terms from the molecular function category in neurons. a bar graph was plotted with the first significant go term nodes based on the p-value. the y-axis shows the top go terms for the molecular function, while the x-axis shows the corresponding –log. the specific values can be found in the additional file additional file 11: figure s top go categories with the highest percentage of neuron associated genes. a bar graph plotted with the first nodes based on the percentage of genes associated with a particular go term. the y-axis shows the top go terms for the molecular function, while the x-axis shows the corresponding percentage of genes. the specific values can be found in the additional file additional file 12: figure s the complete hierarchy showing the cellular component category of go enrichment in neurons. a graph of the cellular component go term analysis of the neuronal genes. the analysis was performed with genes which resulted in enriched go terms and edges. the size of the nodes refers to the number of genes associated with the go term while the color represents the significance of association with the go term. the low significant nodes to the higher significance ranges from white to dark orange nodes. additional file 13: figure s top significant go terms from the cellular component category in neurons. a bar graph with the top most significant nodes based on the p-value. the y-axis shows the top go terms for the cellular component, while the x-axis shows the corresponding –log. the specific values can be found in the additional file additional file 14: figure s top go categories with the largest number of gene association in neurons. a bar graph was plotted with the top nodes based on the percentage of genes associated with a particular go term. the y-axis shows the top go terms for the molecular function, while the x-axis shows the corresponding percentage of genes. the small percentage of gene associated with go terms beyond the 41st category were too low to be plotted. the specific values for the nodes can be found in the additional file additional file 15: figure s the complete hierarchy showing the molecular function category of go enrichment in non-neurons. a go enrichment analysis of molecular function in the non-neuronal genes. the analysis was performed with genes which resulted in enriched go terms and edges. the size of the nodes refers to the number of genes associated with the go term while the color represents the significance of association with the go term. the low significant nodes to the higher significance ranges from white to dark orange nodes. additional file 16: figure s top significant go terms from the molecular function category in the non-neurons. a bar graph was plotted with the top significant go term nodes based on the p-value. the y-axis shows the top go terms for the molecular function, while the x-axis shows the corresponding –log. the specific values can be found in the additional file additional file 17: figure s top go categories with the largest number of gene associated with go terms in the non-neurons. a bar graph plotted with the top nodes representing the highest percentage of genes associated with a specific go term based on the percentage of genes associated. the y-axis shows the top go terms for the molecular function, while the x-axis shows the corresponding percentage of genes. the specific values can be found in the additional file additional file 18: figure s the complete hierarchy showing the cellular component category of go enrichment in non-neurons. a go enrichment analysis of cellular function in the non-neuronal genes. the analysis was performed with genes which resulted in enriched go terms and edges. the size of the nodes refers to the number of genes associated with the go term while the color represents the significance of association with the go term. the low significant nodes to the higher significance ranges from white to dark orange nodes. additional file 19: figure s top significant go terms from the cellular component category in non-neurons. to identify and visualize some of the top significant go terms, a bar graph was plotted with the top nodes based on the p-value. the y-axis shows the top go terms for the cellular component, while the x-axis shows the corresponding –log. the nodes with significant p-values beyond the 25th category were too low to be plotted. the specific values for the nodes can be found in the additional file additional file 20: figure s top go categories with the highest number of gene association in non-neurons. a bar graph was plotted with the top nodes based on the highest number of genes associated with a specific go term. the y-axis shows the top go terms for the molecular function, while the x-axis shows the corresponding percentage of genes. the nodes with significant p-values beyond the 25th category were too low to be plotted. the specific values for the nodes can be found in the additional file additional file 21: table s list of go terms and their corresponding values. table showing the top go terms in the three categories; biological processes, molecular function and cellular components in both the cell populations with their corresponding p-value as well as the percentage of genes associated from the data set. the values have been plotted and the bar graphs can be seen in the additional files additional file 22: table s list of molecular functions and their associated values for the present dataset. table of all the values of the molecular functional enrichment analysis performed on the dataset. the analysis was accomplished with genes. the table gives the number of “molecules” associated with the particular molecular function as well as the percentage of genes that was used to create the pie charts in both the intestinal cell populations. additional file 23: table s list of the canonical pathways in the two intestinal cell populations. pathway analysis on the genes enriched in both the populations to identify the most significant pathways. the analysis was accomplished using genes. the table provides the top canonical pathways with their respective –log and the genes from the dataset that are known to be associated with that pathway. abbreviations bppase pair d.p.fdays post, fertilization degdifferentially expressed genes ensenteric nervous system facsfluorescence, activated cell sorting fdrfalse discovery rate fpkmfragments per kilobase of transcript per million mapped reads gfpgreen fluorescent protein gfp + vegreen fluorescent protein positive gfp-vegreen fluorescent protein negative gigastrointestinal gogene ontology gtpgene transfer format gwasgenome wide association studies h.p.fhours post fertilization hscrhirschprungs ipaingenuity pathway analysis nccneural crest cells qpcrquantitative polymerase chain reaction 
pathway analysis pointed to a high level of g-protein coupled pathway activation, and identified novel roles for candidate pathways such as the nogo/reticulon axon guidance pathway in ens development. despite its critical function, many of the molecular mechanisms that direct the development of the ens and the elaboration of neural network connections remain unknown. thousands of zebrafish intestines were manually dissected from a transgenic line expressing green fluorescent protein under the phox2b regulatory elements [tg our results confirm previously implicated genes and pathways required for ens development, and also identify scores of novel candidate genes and pathways. analysis revealed genes and pathways that may function in ens cell determination, genes that may be identifiers of different ens subtypes, and genes that define the non-neural cellular microenvironment of the ens. motility and the coordination of moving food through the gastrointestinal tract rely on a complex network of neurons known as the enteric nervous system. differential expression analysis between the two cell populations revealed the expected neuronal nature of the phox2b expressing lineage including the enrichment for genes required for neurogenesis and synaptogenesis, and identified many novel genes not previously associated with ens development. because ens development is tightly linked to its environment, the transcriptional landscape of the cellular environment of the intestine was also analyzed. thus, our dataset suggests various potential mechanisms that drive ens development facilitating characterization and discovery of novel therapeutic strategies to improve gastrointestinal disorders. the goal of this study was to transcriptionally identify molecular pathways and candidate genes that drive specification, differentiation and the neural circuitry of specific neural progenitors, the phox2b expressing ens cell lineage, during normal enteric nervous system development. we report the comprehensive gene expression profiles of a lineage-specific population of enteric progenitors, their derivatives, and their microenvironment during normal enteric nervous system development. enteric nervous systemneural cresttranscriptomerna-sequencingzebrafishhirschsprungsphox2bhttp://dx.doi.org/ /100001024roy j. carver charitable trustid0etabg19824kuhlman julie a. http://dx.doi.org/ /100009227iowa state universityid0e2hbg19826kuhlman julie a. http://dx.doi.org/ /100005280american association of university womenid0epgbg19825roy-carson sweta issue-copyright-statement© the author 2017
 motility and the coordination of moving food through the gastrointestinal tract rely on a complex network of neurons known as the enteric nervous system. despite its critical function, many of the molecular mechanisms that direct the development of the ens and the elaboration of neural network connections remain unknown. the goal of this study was to transcriptionally identify molecular pathways and candidate genes that drive specification, differentiation and the neural circuitry of specific neural progenitors, the phox2b expressing ens cell lineage, during normal enteric nervous system development. because ens development is tightly linked to its environment, the transcriptional landscape of the cellular environment of the intestine was also analyzed.  thousands of zebrafish intestines were manually dissected from a transgenic line expressing green fluorescent protein under the phox2b regulatory elements [tg w37]. fluorescence-activated cell sorting was used to separate gfp-positive phox2b expressing ens progenitor and derivatives from gfp-negative intestinal cells. rna-seq was performed to obtain accurate, reproducible transcriptional profiles and the unbiased detection of low level transcripts. analysis revealed genes and pathways that may function in ens cell determination, genes that may be identifiers of different ens subtypes, and genes that define the non-neural cellular microenvironment of the ens. differential expression analysis between the two cell populations revealed the expected neuronal nature of the phox2b expressing lineage including the enrichment for genes required for neurogenesis and synaptogenesis, and identified many novel genes not previously associated with ens development. pathway analysis pointed to a high level of g-protein coupled pathway activation, and identified novel roles for candidate pathways such as the nogo/reticulon axon guidance pathway in ens development.  we report the comprehensive gene expression profiles of a lineage-specific population of enteric progenitors, their derivatives, and their microenvironment during normal enteric nervous system development. our results confirm previously implicated genes and pathways required for ens development, and also identify scores of novel candidate genes and pathways. thus, our dataset suggests various potential mechanisms that drive ens development facilitating characterization and discovery of novel therapeutic strategies to improve gastrointestinal disorders. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords enteric nervous systemneural cresttranscriptomerna-sequencingzebrafishhirschsprungsphox2bhttp://dx.doi.org/ /100001024roy j. carver charitable trustid0etabg19824kuhlman julie a. http://dx.doi.org/ /100009227iowa state universityid0e2hbg19826kuhlman julie a. http://dx.doi.org/ /100005280american association of university womenid0epgbg19825roy-carson sweta issue-copyright-statement© the author 2017 