Regional cytoarchitecture of the adult and developing mouse enteric nervous system

The enteric nervous system (ENS) populates the gastrointestinal (GI) tract and controls GI function. In contrast to the central nervous system, macrostructure of the ENS has been largely overlooked. Here, we visually and computationally demonstrate that the ENS is organized in circumferential stripes that regionally differ in development and neuronal composition. This characterization provides a blueprint for future understanding of region-specific GI function and identifying ENS structural correlates of GI disorders.


Introduction
The gastrointestinal (GI) tract is the only peripheral organ with a dedicated nervous system, which is sufficient to autonomously control GI function. Neurons of the enteric nervous system (ENS) are located within either the myenteric plexus (MP), which controls motility, or the submucosal plexus (SMP), which regulates luminal secretion and absorption. These plexuses develop from migrating neural crest cells (NCCs), which differentiate and cluster into ganglia.
Disrupting ganglia formation, colonization, or neuronal tract development, such as through genetic knockout of cell adhesion genes, can have profound impacts on ENS organization and GI function [1][2][3] .
Particular digestive functions occur in specific GI regions, thus local structure and ENS composition may differ correspondingly to achieve these functions 4,5 . Nevertheless, unlike precisely defined pathways and structures within the brain or spinal cord 6,7 , how neurons are organized within plexuses above the level of ganglia, and how organization differs along the length of the GI tract, are poorly characterized. Here, we present a bird's eye view of the mouse ENS, describing the previously underappreciated regional organization of the ENS, the development of this organization, and the distribution of neuronal subpopulations that comprise it.

Results
To understand the broad layout of neurons in the ENS, we imaged large areas (up to ~50 mm 2 ) of both MP and SMP in each intestinal region: duodenum, jejunum, ileum, proximal colon and distal colon. Immunolabelling neuron somas in the MP revealed a circumferential orientation of ganglia, which in turn loosely coalesced into a noncontinuous stripe pattern perpendicular to its longitudinal axis, herein referred to as neuronal stripes (Fig. 1a). In contrast, submucosal ganglia of the small intestine (SI) did not show any such arrangement, nor did they correlate with the location or organization of epithelial crypts (Fig. 1b). Ganglia of the proximal colon, on the other hand, were organized into diagonal stripes, converging directly opposite the mesenteric border and aligning with the mucosal folds of the proximal colon (Fig. 1b). Full thickness wholemount preparations did not suggest any structural interplexus relationship, and highlighted the orientation difference between the MP and SMP in the proximal colon (Fig. 1c).
The stripe formations of the MP emerged as peaks in longitudinal axis signal intensity profiles of wholemount images, which also confirmed the lack of circumferential structure in the SMP (Fig.   1d). We confirmed this ENS organization using a separate imaging approach, synchrotron source x-ray tomography, performed on intact intestine that allows simultaneous visualization of multiple cell types of both MP and SMP (Fig. 1e,f). Taken together, the MP, but not the SMP, displays a striped neuronal architecture.
We next sought to characterize and quantify this region-and plexus-specific ENS structure.
While intensity profiles (Fig. 1d) are useful visualizations, they cannot reveal any structure that is not perpendicular to the x-axis, such as in the proximal colon SMP, and they are susceptible to distortions in tissue preparation. We therefore used conditional intensity function (CIF) plots to generate spatial probability maps of neuronal locations relative to a given neuron ( Fig. 1g; see Methods). The analysis revealed that an average neuron was always situated within a stripe and was flanked by higher order stripes, clearly displayed when the CIF is collapsed into a twodimensional (2D) graph (Fig. 1g, bottom). Quantification of these graphs revealed that the proximal colon contained the thickest stripes (Fig. 1h), while the largest interstripe distances were found in the SI (Fig. 1i). This result mirrored analysis of regional neuronal density, in which the proximal colon MP was found to be the most dense, and the duodenum the least (Fig. 1j).
The SMP was more sparsely populated than the MP, and, in contrast to the MP, the densest region of the SMP was the duodenum (Fig. 1k). Further, CIF analysis revealed that in the SMP only the proximal colon exhibited any macrostructure beyond ganglia organization (Fig. 1l).
Therefore, the organization of enteric neurons quantitatively differs between regions and plexuses of the GI tract.
To assess how neuronal stripes arise, we visualized MP neurons in intestinal wholemount preparations from late embryonic through neonatal ages ( Fig. 2a and Extended Data Fig. 1a). At embryonic day (E)14.5, when enteric neurons have populated the entire length of the mouse intestines 8 , neurons were scattered in both the SI and colon (Fig. 2a-d). Beginning at E16.5 in the jejunum and E18.5 in the distal colon, MP neurons reorganized from scattered cells into circumferential neuronal stripes, which resolved into individual ganglia at postnatal stages ( Fig.   2a,b). CIF analysis confirmed the gradual emergence of neuronal stripes in the developing ENS and the different organizational timelines between the SI and colon (Fig. 2b). To quantify the progressive organization of the developing ENS, we analyzed nearest-neighbor distances and compared them to synthetically generated data with an imposed minimum nearest-neighbor distance of 10 µm, approximately the diameter of an enteric neuron (Fig. 2c). At E14.5, the dispersion of MP neurons did not differ compared to random distributions; this dispersion transitioned to highly clustered at E16.5 in the duodenum and jejunum and at early postnatal stages in the ileum and colon (Fig. 2d). Thus, neuronal stripes arise from a gradual reorganization of neurons, which occurs embryonically in the more proximal intestine and neonatally in the more distal intestine.
We next assessed whether gut growth or neuronal cell death contribute to the emergence of neuronal stripes. Intestinal length increased ten-fold between E14.5 and postnatal day (P)21, which correlated with decreased neuron density and increased interstripe distance (Fig. 2e-j and Extended Data Fig. 1b-g). Expression of the apoptotic marker Caspase3 revealed sparse apoptotic neurons in the developing MP (<0.5% of HuC/D+ enteric neurons across regions and developmental time, data not shown). Collectively, these results suggest that gut growth, but not neuronal cell death, may influence organization of the developing MP.
We next sought to assess the regional distribution of neuronal subtypes, focusing on known determinants of subtype function, such as Ca2+ binding protein (CBP; Fig. 3a-c) neurotransmitter ( Fig. 3d-g), and neuropeptide ( Fig. 3h-k) expression [9][10][11] . We typically observed a greater proportion of total neurons expressing a given marker in the SI than the colon. There were particularly pronounced differences for CBPs calbindin and secretagogin (Fig. 3b Both PenkCre-tdT and Tac1Cre-tdT also displayed some non-neuronal tdTomato expression (Extended Data Fig. 4g,j). Interestingly, somatostatin was the only marker more prevalent in the colon than in the SI (Fig. 3k,v and Extended Data Fig. 4d). It remains unclear which neuronal subtypes predominate in the distal colon.
We also assessed whether marker expression differed between regions within the SI or colon.
Despite distribution differences in cell body localizations, staining of neuron fibers tended to be qualitatively similar between regions, for example ChATCre-tdT fibers were observed in circular muscle in both SI and colon. The exception to this was secretagogin. In the SI, secretagogin fibers were restricted to interganglionic tracts, while they additionally innervated longitudinal and circular muscle in the proximal and distal colon, respectively ( Fig. 3c and Extended Data Fig.   2c). This may suggest a different, as yet uncharacterized, function for secretagogin neurons, dependent on intestinal region. We next examined neuronal subtype representation across myenteric neuronal stripes. We selected 4 markers for this analysis, covering a range of different cell types: motor neurons (VipCre-tdT and calretinin), sensory neurons (calbindin) and interneurons (somatostatin). At the level of enteric neuronal stripes, the proportion of neurons positive for a given marker broadly reflected that of the entire region, and these distributions shifted appropriately between regions ( Fig. 3w). For instance, VipCre-tdT is in ~15% of ileum neurons overall (Fig. 3u), and each stripe was found to have 0-30% VipCre-tdT neurons, with a peak at 15% (Fig. 3w). Subtypes appear to be well distributed across stripes, with few stripes containing zero neurons of a given subtype unless that marker is expressed in very few cells, such as calbindin (Fig. 3b,w). This is further illustrated when comparing longitudinal axis intensity profiles of HuC/D signal with that of a given subtype ( Fig. 3x and Extended Data Fig. 5). Only at very low levels did the correlation between profiles vanish, as seen in the somatostatin ileum profile compared with the colon profile ( Fig. 3x). Indeed, the strength of this correlation greatly increased over a small range when plotted against the proportion of neurons positive for a given marker (Extended Data Fig.   5d). This suggests that, even with a restricted area of analysis (1800 µm 2 ), neuron subtype prevalence can be low (2-3%), yet its spatial distribution within a region can reflect that of HuC/D, providing evidence for an approximately even distribution of a given subtype across enteric neuronal stripes.

Discussion
By combining large-scale image analysis with computational methods over multiple regions and ages of the mouse intestine, we demonstrate that the mouse MP possesses a gross cytoarchitecture of circumferential neuronal stripes that differ by region in their organization, development, and neuronal composition.
Differences in ENS organization and composition may contribute to the unique functions within each GI region. Single-cell RNA sequencing has revealed differential distribution of functional neuronal classes throughout the GI tract 12,13 , observations confirmed by our results and extended to the level of neuronal stripes. This diversity differs more region-to-region than among neuronal stripes of a given region, which tend to possess similar neuronal compositions to their neighbors.
The role of individual stripes remains to be elucidated. Clonally related ENS cells inhabit overlapping domains and exhibit coordinated activity 14 ; whether these functional units map onto the macrostructure of neuronal stripes has yet to be explored. Further, stripes may serve to control rings of circular muscle, which possess the same orientation as neuronal stripes, or to coordinate longitudinal interneuron projections oriented in parallel 15 .
Another outstanding question is the mechanism that underlies the development of neuronal stripes. The spatiotemporal differences we observe in neuronal stripe development may reflect timing of NCC colonization. Development of circular muscle may also influence this patterning.
The circular muscle differentiates around the time of neuronal stripe emergence 16,17 , and inhibition of circular muscle contractions changes the anisotropy of the embryonic chick ENS 18 .
Such mechanistic insights may improve methods for generating intestinal organoids, examples of which have two ENS plexuses but lack patterning along the longitudinal axis 19,20 .
Our characterization of regional enteric organization and composition can serve as a blueprint to understand how local enteric neurocircuitry underlies region-specific motility patterns and function. Further, assessing ENS structure may help identify anatomical changes in human GI diseases with no known correlates within the ENS, such as intestinal pseudo-obstruction 21 .
Comprehensive analyses of central nervous system structure have served as platforms to dissect basic physiology and disease, and our work represents a parallel step toward an intricate understanding of the ENS.

Marker selection
For experiments assessing neuronal subtype distribution, we selected a variety of markers (via either genetically encoded reporters, above, or antibodies (Table1)) covering calcium binding proteins, neurotransmitters and neuropeptides, which are known determinants of enteric subtype function or which delineate functionally relevant populations 13,22 . Furthermore, all labels chosen marked the neuron soma, which facilitated automated counting of large areas of tissue.
Calretinin, ChAT, and Tac1 are all known to mark excitatory motor neurons. Penk also marks excitatory motor populations, in addition to some interneurons. Somatostatin has been suggested to delineate a population of interneurons 23 . Inhibitory motor neurons were labelled by VIP and nNOS. Calbindin is suggested to mark sensory neurons 24 . The functions of VGLUT2 and Gad2 neurons are not known, but were included as they are well-known neurotransmitter markers in other areas of the nervous system and had been identified in the ENS by recent scRNAseq studies as belonging to specific subpopulations 13,25 . Secretagogin, a calcium binding protein first identified in the pancreas 26,27 , was similarly identified recently, and has not yet been described in the ENS.

Immunohistochemistry
For investigations into the structure of the adult ENS, the small intestine and colon were dissected out from 2-6 months-old mice culled by CO2 and cervical dislocation. Intestines were flushed of fecal contents using cold PBS before cutting the proximal, middle and distal 2-3 cm of the small intestine to isolate the duodenum, jejunum and ileum. The colon was cut in two, and all segments were placed in cold PBS on ice. A fire-polished glass rod 2 mm in diameter was placed through each segment, and the mesentery removed. For full-thickness wholemount preparations, tissue was fixed on the glass rod in 4% PFA for 90 minutes at 4°C with shaking.
For adult single-plexus wholemount preparations, the muscularis (both muscle layers and MP) was peeled away using cotton buds, starting at the mesenteric border. For neonatal ages (postnatal days 0, 10, 15, and 21), intestines were cut open along the mesenteric border, pinned mucosa-down on Sylgard 170 in ice-cold PBS, and the muscularis was peeled away with fine forceps. Muscularis tissue was stored in ice-cold PBS and then pinned onto Sylgard 170 in a flat sheet in 35 mm glass dishes using insect pins. In some cases, the underlying mucosa containing the SMP was also kept and pinned. Tissue was immediately fixed as described above. Residual PFA was removed via PBS washes before proceeding directly to immunostaining, or stored in PBS containing 0.1% sodium azide for up to 3 weeks at 4°C.
For embryonic experiments, pregnant mice were culled by CO2 and cervical dislocation. The uterus was removed, and embryos were placed in ice-cold PBS. The intestines were dissected from each embryo, and the mesentery was carefully removed. To measure intestinal length, the cleaned intestine was laid adjacent to a ruler. In a Sylgard 170 plate, a pin was placed in the stomach and anus to keep the intestines taut and straight. Tissue was fixed in 4% PFA for 90 minutes at 4°C with shaking.
For immunohistochemistry, large intact pieces of tissue were cut into smaller pieces and transferred to WHO microtitration trays (International Scientific Supplies) containing PBS.
The following day, tissue was washed 3 times in PBT for 30 minutes each before transferring to PBT containing secondary antibodies ( Table 2) for 2 h at room temperature with shaking. Tissue was washed twice in PBT and twice in PBS, then mounted onto slides, rinsed in ddH2O and coverslipped using Fluoromount-G (Southern Biotech).
For embryonic experiments, immunohistochemistry protocol was done as in adult with the following changes. Intact intestines were placed in WHO microtitration trays (International Scientific Supplies) and treated as above. After immunohistochemistry, fixed intestines were cannulated with a cleaning wire for 33-gauge needles (Hamilton) and cut along the wire.
Embryonic intestines were mounted full-thickness.   Briefly, the tissues were extensively washed in 0.1 M sodium cacodylate buffer (EMS). This was followed by sequential staining with 2% buffered osmium tetroxide (EMS), 2.5% potassium ferrocyanide (Sigma-Aldrich) with no rinse in between followed by pyrogallol (Sigma-Aldrich), unbuffered 2% osmium tetroxide (EMS), 1% uranyl acetate (EMS), and 0.66% aspartic acid buffered lead (II) nitrate (Sigma-Aldrich) with extensive rinses in between each of the steps. The stained tissues were dehydrated in graded ethanol, propylene oxide and infiltrated with epon resin (EMS). The tissues were finally embedded in fresh epon and cured in an oven at 60°C for 48 h 29,30 .

Synchrotron imaging
The epon embedded samples were imaged at 32-ID beamline at the Advanced Photon Source (APS) in Argonne National Laboratory. The sample was placed on a rotation stage and projection images were acquired at 600 nm/pixel resolution as the sample was rotated over 180 degrees 31,32 . The acquired images were 3D-reconstructed using the TomoPy toolbox 33 .

Data Analysis
The 3D-reconstructed data were manually annotated using Knossos, an open-source software (https://knossos.app/) 34 . We identified the myenteric ganglia in between the circular and longitudinal muscle layers and the submucosal ganglia in the submucosal space. All cell bodies were identified based on the circular outline and the presence of a nucleus within the cell.

Image acquisition
Images were acquired using a 20x (NA 0.75) oil objective on a Leica SP8 confocal microscope.
Large tiled images were acquired and stitched using the Navigator mode within LASX. Z-stacks were acquired with 3 µm between each focal plane for adult tissue and with 2 µm between each focal plane for embryonic and neonatal tissue.

Cell counting
Image analysis was primarily performed using ImageJ/FIJI (NIH, Bethesda, MD). For neuronal density analysis, HuC/D images (maximum intensity projections for adults, single plane for perinatal) were blurred before thresholding and watershedding. Cells were counted using the Analyze Particles function, and density was calculated based on area of tissue measured in FIJI. To count neuronal subtypes, the Image Calculator function was used to combine thresholded HuC/D stacks (prior to maximum projection) with raw image stacks of a defined cell type. The result of this calculation was then maximally projected and counted as above.

Conditional intensity function (CIF) analysis
To visualize and evaluate spatial patterns, we calculated the conditional intensity function (CIF), which generates a spatial density map of neuron locations relative to a given neuron. Square images measuring ~1800x1800 µm for adult, ~400x400 µm for embryonic and neonatal tissue, respectively, were processed in FIJI as described above, and the XY coordinates of each neuron were obtained using Analyze Particles. We empirically estimated the CIF for a given sample by iterating over all neurons and calculating the number of neurons in a 2D grid around that neuron. Total image area was normalized to 1, divided into 100 bins per unit length. The 2D grid's width and height were both 0.7 for adult data, while they were 0.8 and 0.5, respectively, for embryonic/neonatal data. Density values were normalized to expected density based on a uniform distribution of neurons, given a value of 1. We excluded the center grid point from the resulting CIF plot, which included data from the neuron used for conditioning.
We then transformed the 2D grid into a one-dimensional line by averaging along the y-axis, using either the full y-axis (for interstripe distance calculations) or a smaller proportion (for stripe width calculations). For adult data, this smaller proportion was a length of 0.1 relative image length above and below the center, and for developmental data, we used a value of 0.2. For developmental data, which contained far fewer neurons, we smoothed with a Gaussian of 20 μm standard deviation. We then identified the first minima and first peak next to the center.
Stripe width was taken as the width at the half height from minimum to center peak. Interstripe distance was taken as the distance between the left and center peaks. For interstripe distance analysis, samples in which secondary stripes could not be unambiguously identified were not included. All analyses were implemented in Python.

Nearest-neighbor and Empirical distribution function (EDF) analysis
To assess enteric neuron organization across development, we tested each of the samples for deviation from a hypothesis of random neuron positions, known as complete spatial randomness. Our analysis was based on nearest neighbor distances. The distribution of nearest neighbor distances under complete spatial randomness can be calculated 35 . Specifically, the mean of the distribution and its variance can be approximated by a normal distribution 35 .
Embryonic and neonatal images were processed as above to generate XY coordinates, and then for each sample we calculated the mean of the data and extracted a z-score for the deviation of the data from the expected value. In addition, to better visualize deviations, we generated synthetic samples under the assumption of complete spatial randomness. We then plotted the empirical distribution function for the data as well as an envelope defined by the maximum and minimum of the empirical distribution function over 500 synthetic samples.
We observed that some deviation from randomness was driven mostly by the fact that the synthetically generated samples exhibited overlap in space while our data, taken from a singleplane image, did not overlap. To make the analysis more robust to this property, we generated synthetic samples with a minimum distance imposed. We chose this minimum distance as 10 µm, approximately the average diameter of an enteric neuron. We then generated 500 random samples and calculated their z-score and defined a metric that is the difference between the zscore of the data and the mean of the z-score of the synthetic samples. This allowed better identification of differences from other spatial properties (e.g., stripes). All analyses were implemented in Python. Microsoft Excel and peaks in the data were automatically identified. The locations of these peaks were then exported to FIJI, which created a box measuring the full y-axis and 50 µm either side of each peak location. This box would constitute a stripe, and Analyze Particles was used on thresholded images (as above) to count the number of objects within each stripe. This data was converted to frequency distributions in Prism 9.

Statistical analysis
Statistical tests and graphical representation of data were performed using Prism 9 software

Data availability
The data that support the findings of this study are available upon reasonable request.    Figure 1a. b, CIF plots of wholemount MP HuC/D immunohistochemistry in the jejunum and distal colon over developmental time. Yellow: high probability density; blue: low. Axes apply to all panels. c, 2D arrays of neurons (HuC/D+ cells) and corresponding synthetically generated data for the E16.5 and P10 distal colon. Empirical distribution function (EDF) plots for neurons (blue line) and synthetically generated maximum and minimum (gray lines). Inset values indicate z-score difference between data and mean of synthetic values for a representative sample. d, z-score difference between data and mean of synthetically generated values (mean ± SEM) derived from EDF plots in all regions of the intestine over developmental time. n = 2-11. e-f, Lengths of small intestine (e) and colon (f) (mean ± SEM). n = 3-15. gh, Neuronal density (mean ± SEM) in the jejunum (g) and distal colon (h). n = 3-11. i-j, Violin plots of interstripe distance (mean ± SEM) in the jejunum (i) and distal colon (j) as analyzed by CIFs, n = 3-11. k, Representative images of immunohistochemical labelling with HuC/D (magenta) and apoptotic marker Caspase3 (green) in wholemount preparations of the jejunum and distal colon MP. All tests one-way ANOVA. Pairwise comparisons not shown. ***p<0.001, ****p<0.0001. Scale bars as indicated. AU, arbitrary units; EDF, empirical distribution function; CIF, conditional intensity function; E, embryonic day; MP, myenteric plexus; P, postnatal day.