Introduction

The enteric nervous system (ENS) is a large collection of neurons, glial, and precursor cells that resides within the gastrointestinal wall and regulates gut motility and secretion along with modulating epithelial and immune cell function^{1, 2}. During fetal development, the mammalian ENS is populated by neurons and glia derived from neural crest (NC)-derived precursors^3–9. These precursors follow diverse migratory routes to colonize and innervate various parts of the gut before birth^10–12. It is not clear, however, that this lineage persists in its entirety in the adult gut, as indicated by the observed lack of expression of fluorescent reporter protein in a subpopulation of adult enteric neurons in NC-lineage-traced mice^{13, 14}. Alternative sources of enteric neurons that have been proposed in the literature include the ventral neural tube (VENT)¹⁵, or the Pdx1- expressing pancreatic endoderm¹⁴, but the interpretation of these studies has been limited by the lack of robust lineage markers for non-NC derived neurons¹⁶. In addition, while prior studies have documented cellular changes to the ageing ENS¹⁷, the developmental mechanisms behind these changes are unknown. Thus, confirmation of a second, distinct lineage of enteric neurons in adults is important for our understanding of the healthy post-natal development and aging of the ENS, as well as for the pathogenesis of acquired disorders of the ENS.

In this study, we found that while the early post-natal ENS is derived from the canonical NC- lineage, this pattern changes rapidly as the ENS matures, due to the arrival and continual expansion of a novel population of Mesoderm-derived Enteric Neurons (MENs) which represent an equal proportion of the ENS in young adulthood and with increasing age, eventually outnumber the NC-derived Enteric Neurons (NENs). We also found that, while the NEN population is regulated by glial derived neurotrophic factor (GDNF) signaling through its receptor RET, the MEN population is regulated by hepatocyte growth factor (HGF) signaling. Increasing HGF levels during maturation or by pharmacological dosing increase proportions of MENs. Similarly, decrease in GDNF with age decrease NENs; and increasing GDNF levels by pharmacological dosing increase NENs proportions in the adult ENS to impact intestinal motility.

These results indicate for the first time that the mesoderm is an important source of neurons in the second largest nervous system of the body. The increasing proportion of neurons of mesodermal lineage is a natural consequence of maturation and aging; further, this lineage can be expected to have vulnerabilities to disease that are distinct from those affecting the NEN population. These findings therefore provide a new paradigm for understanding the structure and function of the adult and aging ENS in health, age-related gut dysfunction and other acquired disorders of gastrointestinal motility.

Results

Only half of all mid-age adult enteric neurons are derived from the neural crest

We analyzed small intestinal longitudinal muscle-myenteric plexus (LM-MP) from adult (post-natal day 60; P60) Wnt1-cre:Rosa26-tdTomato mice, in which tdTomato is expressed by all derivatives of Wnt1⁺ NC-cells¹⁸. In these tissues, while GFAP, a glial marker, was always co-expressed with tdTomato (Fig 1a), tdTomato-expression was absent in many myenteric neurons (Fig 1b, c). By careful enumeration, almost half of all myenteric neurons expressing the pan-neuronal marker Hu were found to not express the reporter (percent tdTomato- neurons: 44.27 ± 2.404 SEM; enumerated from 2216 neurons from 6 mice, Fig 1d). In these lineage-traced mice, myenteric ganglia were found to contain tdTomato^high and tdTomato^lowneurons and due care was taken to image subsets (Suppl. Fig 1a). Both tdTomato^high and tdTomato^low neurons were classified as tdTomato+ and only neurons that did not show any tdTomato-expression were classified as tdTomato- neurons. Immunofluorescence staining of small intestinal LM-MP tissue from adult Wnt1-cre:Rosa26-tdTomato mice using antibodies against tdTomato protein showed a consistent lack of reporter expression in a population of myenteric neurons (Suppl. Fig 1b). This lack of reporter expression is not a function of tissue peeling or fixation, as the tdTomato- cells within the myenteric ganglia were also observed in freshly harvested unpeeled gut tissue, when imaged using a live tissue fluorescence microscope (Suppl. Fig 1c). In addition, we noticed that myenteric ganglia of the Wnt1-cre:Rosa26-tdTomato transgenic mouse show similar tdTomato aggregation in cells of the myenteric plexus from both freshly harvested and immediately fixed tissue (Suppl. Fig 1d, e). Since we previously showed ongoing neuronal loss in healthy ENS¹⁹, we reasoned that cells that show hyper-aggregation of tdTomato are in advanced stages of cell death. Cells with hyper-aggregation of tdTomato showed a lack of staining for the nuclear marker DAPI (Suppl. Fig 1e), suggesting that these are indeed cells in the advanced stages of death, and given that they are not labeled with antibodies against Hu, their presence does not alter our estimation of tdTomato-expressing neurons (Suppl. Fig 1f). Thus, the absence of reporter expression in myenteric ganglia of freshly harvested tissue and the presence of tdTomato aggregation in both freshly harvested and fixed tissues shows that our observations are not caused due to any technical issues in tissue isolation and preservation.

Fig 1:

We next sought additional confirmation that the lack of tdTomato expression in this often-used lineage fate mapping mouse model was not due to incorrect activity at the Rosa26 locus where the floxed reporter transgene is located, due to aberrant cre activity under the Wnt1 promoter, or due to issues with the antibodies used to detect the pan-neuronal protein Hu. This was provided by an analogous Wnt1-cre:Hprt-tdTomato lineage-traced mouse line, in which tdTomato was expressed from the Hprt locus in a Wnt1-cre-dependent manner (Hprt locus is X-inactivated in females, hence only adult male mice were used)²⁰; and by the Pax3-cre:Rosa26-tdTomato lineage-traced mouse line, where Rosa26-tdTomato driven by the Pax3-cre transgene labels the derivatives of the neural tube and pre-migratory NC (Suppl. Fig 1g, h)²¹. Similar lack of reporter expression was previously observed by the Pachnis Lab in adult myenteric neurons from Sox10- cre NC-lineage-traced mice, and more recently by the Heuckeroth Lab in significant numbers of myenteric neurons from Wnt1-cre:Rosa26-H2B-mCherry mice, which expresses nuclear localized reporter mCherry in NC-lineage cells, further confirming our observations of the presence of non- NC-derived neurons in the adult ENS^{13, 22}. In addition to the ANNA-1 antisera, which is known to detect all neuronally significant Hu proteins (HuB, HuC, and HuD)²³, we found that another often used antibody thought to be specific to HuC and HuD also detects HuB protein (Suppl. Fig 1l) – showing that these antibodies do not detect the expression of a specific neuronal Hu antigen.

Lineage-tracing confirms a mesodermal derivation for half of all adult myenteric neurons

Alternative sources of enteric neurons proposed previously include the ventral neural tube¹⁵, and the pancreatic endoderm¹⁴, but the interpretation of these studies was limited by the lack of robust lineage markers¹⁶. Brokhman et al.¹⁴ found evidence of labeled neurons in inducible Pdx1-cre, Foxa2-cre, and Sox17-cre lineage-traced mouse lines and inferred their derivation from pancreatic endoderm. However, in a Pdx1-cre lineage-traced mouse line, many neuroectoderm- derived neurons of the central nervous system have also been previously shown to be derived from Pdx1-expressing cells²⁴, suggesting that Pdx1 is not an exclusive endodermal marker. Foxa2 is expressed by mesoderm-derived cells in the heart²⁵ and Sox17 also labels mesoderm-derived cells, especially cells of the intestinal vasculature^26–29. We therefore hypothesized that the embryonic mesoderm may be the true developmental source of the non-NC enteric neurons.

Mesoderm-derived cells during embryogenesis express Tek²⁵ and analysis of LM-MP tissues from adult male Tek-cre:Hprt-tdTomato lineage-traced mice revealed the presence of a population of tdTomato⁺ neurons (Fig 2a). Since Tek is not expressed in the adult ENS (Suppl. Fig 2a), the presence of Tek-derived adult myenteric neurons suggested their mesodermal origin. We then used mesoderm posterior 1 or Mesp1 as a definitive developmental marker for the embryonic mesoderm^30–34 to confirm the mesodermal derivation of the non-NC derived adult enteric neurons. Reporter expression in this mouse line has been previously studied in cardiac development, which we confirmed by observing its expression in the adult heart tissue (Fig 2b). Expression of tdTomato in vascular cells of both Tek-cre and Mesp1-cre mice was more pronounced than in myenteric neurons (Fig 2a, c, d) and in cardiac myocytes of Mesp1-cre mice (Fig 2b, Suppl. Fig 2b), which reflects the variable expression of the reporter gene in both the neural crest-lineage and mesoderm-lineage specific transgenic mouse lines (Fig 2b, c; Suppl. Figs 1a, 2c). Using an analogous Mesp1-cre:Rosa26-EGFP lineage-traced mouse line, in which EGFP was expressed from the Rosa26 locus in a Mesp1-cre-dependent manner, we next confirmed the derivation of a population of myenteric neurons from Mesp1-expressing mesoderm (Fig 2e, Suppl. Fig 2d, e). Similar to our observations with the tdTomato reporter (Fig 2a-d), we observed lower expression of the reporter EGFP in myenteric neurons as compared to vascular cells (Fig 2e; Suppl. Fig 2d, e). Variable expression of CAG and other pan-cellular promoters in various cell-types and during processes of maturation have previously been reported^35–38. In addition, an earlier study by Agah et al. showed differing degrees of CAG-driven LacZ reporter activation in a cardiac-specific transgenic cre mouse line, which were unrelated to copy number, suggesting insertional and positional effects or, potentially, differential methylation³⁹. These reports are consistent with our observation and could potentially help explain the observed variable expression of CAG and CMV promoter-driven reporter genes in our study.

Fig 2:

Analysis of small intestinal LM-MP from P60 Mesp1-cre:Rosa26-tdTomato lineage-traced mice showed that tdTomato-expression was observed in about half of all myenteric neurons (Fig 2f; percent tdTomato⁺ neurons: 50.28 ± 2.89 SEM; enumerated from 484 neurons from 3 mice). Because of the variable expression of the reporter gene in the transgenic mice described above, we also labeled Mesp1-cre:Rosa26-tdTomato+ neurons with the S46 antibody, which was raised against the slow tonic myosin heavy chain protein MHCst derived from avian embryonic upper leg muscle, and is expressed exclusively by mesoderm-derived cell populations (Fig 2g)^40–43. MHCst immunostaining was exclusively observed in all Tek-cre:Rosa26-tdTomato⁺ MENs (Suppl. Fig 2f). None of the 107 MHCst⁺ intra-ganglionic cells observed (35.67 ± 6.17 SEM MHCst⁺ cells per mouse) across 3 Wnt1-cre:Rosa26-tdTomato mice expressed tdTomato, suggesting that MHCst immunostaining was exclusive to non-NC lineage neurons (Fig 2h, Suppl. Fig 3a shows the same image where the MHCst immunostained ganglia can be observed with and without the tdTomato channel for better visualization). MHCst antibody was also found to label extra-ganglionic smooth muscle cells in the LM-MP tissue (Suppl. Fig 3a). Thus, both lineage tracing and protein biomarker expression provides strong support for their mesodermal origin (Fig 2g, Suppl. Fig 3a). Along with our observations on Wnt1-cre:Rosa26-tdTomato mice (Fig 1c), our results indicate that these two distinct lineages together account for all the adult small intestinal myenteric neurons.

The proteins RET, a receptor tyrosine kinase that transduces GDNF signaling in NC-derived enteric neuronal precursors, and MET, a receptor for hepatocyte growth factor (HGF), are expressed by different subsets of adult myenteric neurons⁴⁴. MET is classically expressed by mesoderm-derived cells⁴⁵, and by using immunostaining of small intestinal LM-MP tissues from Wnt1-cre:Rosa26-tdTomato and Mesp1-cre:Rosa26-tdTomato mice, we found that the expression of MET was restricted to a sub-population of adult MENs (Fig 3a, b). By contrast, RET expression was confined to NENs (Fig 3c).

Fig 3:

We then studied whether MENs and NENs differed phenotypically. In the current paradigm regarding the neurochemical basis for the functional classification of adult enteric neurons, inhibitory enteric motor neurons express the nitric oxide-producing enzyme nitric oxide synthase 1 (NOS1), excitatory motor neurons express the acetylcholine-producing enzyme choline acetyl transferase (ChAT), and enteric sensory neurons, called intrinsic primary afferent neurons (IPANs) express the neuropeptide calcitonin gene related peptide (CGRP)^{9, 44, 46}. Immunostaining using antibodies against NOS1, ChAT, and CGRP in the LM-MP from P60 Wnt1-cre:Rosa26- tdTomato mice showed that both tdTomato+ NENs and tdTomato- MENs express these neuronal markers (Fig 3d, f, h). Quantification of proportions in these neuronal lineages shows that the majority of NOS1+ inhibitory neurons and ChAT⁺ excitatory neurons are NENs (Fig 3e; of the 616 NOS1⁺ neurons observed across 3 mice, NOS1+ tdTomato-: 27.06% ± 4.46 SEM; p < 0.0001; Fisher’s exact test); (Fig 3g; of the 912 ChAT+ neurons observed across 3 mice, ChAT⁺ tdTomato- : 26.34% ± 2.54 SEM; p < 0.0001; Fisher’s exact test). By contrast, the majority (∼75%) of CGRP+ neurons were found to be MENs (Fig 3i; of the 146 CGRP⁺ neurons observed across 3 mice, CGRP+ tdTomato-: 75.76% ± 0.63 SEM; p < 0.0001; Fisher’s exact test). These results are in keeping with a previous report by Avetisyan et al. that shows low expression of NOS1 (0% of MET⁺ neurons were NOS⁺) and ChAT (<8% of MET⁺ neurons were ChAT⁺) and abundant expression of CGRP by MET+ neurons in the adult ENS⁴⁴, previously not known to be derived from a different lineage. In addition, MENs also express Cadherin-3 (CDH3, Suppl. Fig 3b), which is known to also mark a sub-population of mechanosensory spinal neurons⁴⁷. Apart from their derivation from a distinct developmental lineage, the mean cell size of MENs was significantly larger than that of NENs (Suppl. Fig 4a; n = 143 neurons/group observed across 3 mice; Feret diameter (µm) MENs: 17.47 ± 0.50 SEM, NENs: 13.03 ± 0.36 SEM, p = 0.002; Student’s t-test).

An expanded molecular characterization of MENs using unbiased single cell RNA sequencing (scRNAseq)-based analyses

Recent studies have used single cell or single nuclei-based transcriptomic analyses to query molecular diversity of enteric neurons. These studies have used either used enrichment methods based on FACS sorting neural crest-derived cells (using reporter-expression from neural crest-specific lineage fate mapping mice)^{22, 48, 49}, Phox2b- expressing cells⁵⁰, Baf53-derived cells⁵¹ or used classical NENs markers such as Snap25 to drive identification and/or enrichment of enteric neurons⁴⁹. Instead of using these strategies, we relied upon another accepted use of scRNAseq, which allows for sequencing of diverse un-enriched cell populations from a tissue, to discover any novel cell types in that tissue⁵². We performed truly unbiased and agnostic clustering of the scRNAseq data from tissues of two 6-month-old adult male C57BL/6 wildtype mice (Fig 4a; Suppl. Fig 4b, c). We identified clusters of NC-derived cells by exclusive expression of canonical NC markers Ret and Sox10 (which include the clusters of Ret-expressing NENs and Sox10-expressing Neuroglia); and the MENs cluster by its co-expression of the genes Calcb (CGRP), Met, and Cdh3 (Fig 4a, b). With cells from both samples pooled together, we compared 1,737 NC-derived cells across the NENs (77 cells) and Neuroglia (1660 cells) clusters with 2,223 cells in the MENs cluster. A list of top marker genes for every cluster are provided in Supplementary Table 2. It is important to note here that it is well known that enzyme-based dissociation methods are unequally tolerated by diverse cell types, which is known to cause over- or under-representation of several cell types in scRNAseq^{53, 54}. The same is true for neurons, where certain neuronal subtypes may be significantly over- or under- represented in scRNAseq databases, compared to their true representation in tissue⁵⁵. Thus, the relative sizes of the NENs and MENs clusters in the scRNAseq data should not be taken as a reliable indicator of their actual proportions in tissue.

Fig 4:

The neuronal nature of MENs was first established by testing our scRNAseq dataset for expression of genes known to be expressed by enteric or other neurons. We found detectable expression of genes Elavl2 (HuB), Hand2, Pde10a, Vsnl1, Tubb2b, Stmn2, Stx3, and Gpr88 in the MENs as well as in the NC-derived cell clusters (Fig 4b). While expression of Elavl2, Stmn2 and Hand2 has been previously observed in enteric neurons^{56, 57}, the expression of neuronal marker genes Pde10a (expressed by medium spiny neurons in the striatum⁵⁸), Vsnl1 (expressed by hippocampal neurons⁵⁹), Gpr88 (expressed by striatal neurons⁶⁰), Stx3 (expressed by hippocampal neurons⁶¹), and Tubb2b (expressed by adult retinal neurons⁶²) has not been studied. Using antibodies against these markers along with co-staining with anti-Hu ANNA1 antisera to label neurons, we visualized the expression of these neuronal markers in both tdTomato⁺ NENs and tdTomato- MENs of the small intestinal myenteric plexus tissue from adult Wnt1-cre:Rosa26- tdTomato mice (Fig 4c, Suppl. Fig 3l-q). In addition, the MENs cluster was also enriched in previously characterized ENS markers, such as Ntf3 and Il18^{2, 63} (Fig 4b). We validated the MEN- specificity of these markers using immunochemical analyses (Fig 4c; Suppl. Fig 3b - k). In addition, we also detected >40 neuronally significant genes whose expression was unique to or enriched in MENs scRNAseq cluster, when compared with other larger cell clusters of macrophages, smooth muscle cells, neuroglia, and vascular endothelium (Suppl. Fig 5). These genes include neurotransmitter receptor genes (such as Gabra1, Gabra3, Gria3, Grik5, Grind2d, Npy1r); ion channels genes (such as Cacna1a, Cacna1g, Cacnb3, Clcn3); gap junction genes (Gjc1, Gjb5); transient receptor potential channel genes (Trpv4, Trpc1, Trpc6); potassium channel genes (such as Kcnn1, Kcnj8, Kcnd3, Kcnh3, Kcns3); hormone encoding genes (Ghrh, Gnrh1, Nucb2) along with other neuronal genes such as, Prss12 (encoding for Neurotrypsin)⁶⁴, Uchl1 (encoding for pan-neuronal marker PGP9.5)⁶⁵, Cplx2 (encoding for Complexin 2)⁶⁶, Gpm6a (encoding for neuronal membrane glycoprotein M6-A)⁶⁷, and Vamp2 (encoding for Synaptobrevin 2)⁶⁸ (Suppl. Fig 5). Owing to the small size of the NENs cluster, a comparative analysis between these two neuronal lineages for these neuronally significant genes was not possible. However, the discovery of these neuronally significant genes in the MENs cluster, and the subsequent validation of a set of them provides evidence for the neuronal nature of MENs. Further examination of the MENs cluster yielded additional MENs-specific marker genes Slpi, Aebp1, Clic3, Fmo2, Smo, Myl7, and Slc17a9, whose expression by adult enteric neurons has not been previously described and which we also validated (Fig 4b).

Next, we examined the MENs cluster for expression of a mesenchymal gene that would underscore the mesodermal nature of MENs. Decorin (Dcn) is a gene coding for a chondroitin- dermatan sulphate proteoglycan, which is expressed by mesoderm-derived cells, including fibroblasts and smooth muscle cells⁶⁹. Decorin was found to be highly expressed in the MENs cluster in our data (Fig 5a), and not by the neurons in the Zeisel et al. database (Fig 5b)⁴⁸. We used a validated antibody against DCN⁷⁰ and found that this protein, known to be specifically expressed by mesoderm-derived cells, is expressed specifically by tdTomato-negative small intestinal myenteric neurons of the Wnt1-cre:Rosa26-tdTomato lineage fate mapping mouse (Fig 5c).

Fig 5:

These data show that the MENs scRNAseq cluster, which can be identified by its co-expression of neuronal markers (Elavl2, Pde10a, Hand2, etc), validated cell surface marker (Cdh3), and the mesodermal marker Dcn is the same as the population of Hu-expressing neurons that do not express tdTomato in the Wnt1-cre:Rosa26-tdTomato neural crest lineage fate mapping mouse. We thus establish the distinct nature of MENs based on the co-expression of many neuronally significant genes and of mesodermal genes by this cell population of Wnt1-cre:Rosa26-tdTomato- negative neuron. Our experimental strategy of performing scRNAseq on LM-MP cells without any pre-enrichment helps detects this newly characterized cell type. However, it does not allow us to compare and contrast the transcriptomic profiles of MENs and NENs – especially given that we do not know the true size of NENs represented in our dataset. This is exemplified by the fact that cells within the Neuroglia cluster showed expression of glial genes such as Sox10 as well as neuronal genes such as Ncam1, Hand2 and Stmn2, suggesting that NC-derived neurons may be present within both the smaller NENs cluster as well as the larger Neuroglia cluster (Fig 4b).

Detection of MENs in other murine and human datasets

Since Phox2b, the transcription factor expressed by all ENS cells, shows significant higher expression in some neurons than in all ENS glial cells⁷¹, May-Zhang et al.⁵⁰ reasoned that flow sorted nuclei of Phox2b-CFP^high cells from a Phox2b-CFP reporter mouse would provide for an enriched population of ENS neurons that can be used for single nuclei RNA sequencing (snRNAseq). They first performed bulk- RNAseq on Phox2b-CFP⁺ nuclei that were flow sorted based on CFP-expression level into CFP^lowand CFP^high fractions, of which the CFP^lowfraction was found enriched with glial-specific gene expression profile. May-Zhang et al.⁵⁰ used this to provide a rationale for performing snRNAseq on Phox2b-CFP^high nuclei, which they presumed to consist of all enteric neurons. Since, MENs make up roughly half of all adult ENS neurons, we hypothesized that they should also express Phox2b. By using a PHOX2b-specific antibody, we validated the expression of this important transcription factor in MENs (Fig 5d). Next, we tested whether we could detect some of the genes that show enriched expression in MENs (namely, Smo, Aebp1, Cdh3, Fmo2, Il18, Slpi, Upk3b, Msln, Hand2) in the Phox2b⁺ bulkRNAseq data from May-Zhang et al.⁵⁰, and found that not only were these genes expressed by their Phox2b⁺ cells, but expression of these genes was also significantly enriched in the Phox2b-CFP^low fraction (Fig 5e) that was not used for subsequent snRNAseq experiments. These data provide evidence that MENs express Phox2b, both transcriptomically, as well as at the protein level. We next queried whether any MENs were represented in the May-Zhang et al.’s snRNAseq data generated on Phox2b-CFP^high nuclei, given that they reported on the presence of a cluster of adult enteric neurons that expressed mesenchymal markers. We again tested whether the MENs-expressed genes, such as Decorin, etc were also expressed by the ‘mesenchymal’ neurons in the May-Zhang et al.⁵⁰ dataset. We found a high degree of similarity between the top genes expressed by the MENs cluster in our scRNAseq dataset and the ‘mesenchymal’ neurons in the May-Zhang et al.’s snRNAseq dataset (Table 1, Suppl. Fig 4d). Furthermore, examination of the top genes shows that this Decorin- expressing neuronal sub-population is significantly different from the other Calcb-expressing neuronal populations, which we presume to be the CGRP-expressing NENs (Table 1).

Table 1:

To confirm that the ‘mesenchymal’ neurons identified by May-Zhang et al.⁵⁰ are MENs, we used an agnostic bioinformatics projection-based approach which we have previously published⁷². This approach allows us to learn latent space representations of gene expression patterns using single cell transcriptomic data, generating patterns that correspond to cell-type-specific gene expression signatures⁷². We can then project other single cell and bulk RNA transcriptomic datasets into these learned patterns to accurately quantify the differential use of these transcriptional signatures in target data across platforms, tissues, and species⁷². In this instance, gene expression patterns were learned on our murine single cell transcriptomic data using non-negative matrix factorization (NMF) dimensionality reduction and the data was decomposed into 50 distinct NMF-patterns^{72, 73} (Fig 6a). Four NMF patterns were identified that were specific to MENs (Fig 6a, b; Suppl. Fig 4e). We next projected the snRNAseq dataset from May-Zhang et al.⁵⁰ into the four MEN-specific patterns and found that all these patterns specifically labeled the cluster of ‘mesenchymal’ neurons in the May-Zhang et al. dataset of the adult ileal ENS (Fig 6b). Further, upon clustering the nuclei based on their projection weights from the four MEN-specific patterns, we identified a cluster of nuclei that showed higher usage of all four patterns, when compared to all the other nuclei sequenced in the dataset (Suppl. Fig 4f). The genes that define this cluster of cells are similar to the top-expressed genes that define our MENs cluster (Suppl. Fig 4g, Table 1). This approach, together with the high similarities in the top-expressed genes, establishes that cluster of neurons annotated as ‘mesenchymal’ neurons by May-Zhang et al. are indeed MENs. Further, one of the markers used by May-Zhang et al. to describe their mesenchymal neuronal cluster was the gene Myh11⁵⁰, which is known to be expressed by mesoderm-derived smooth muscle cells in the adult gut⁷⁴. To test whether Myh11 is expressed by MENs, we immunolabeled the small intestinal LM-MP tissue from an adult Wnt1-cre:Rosa26-tdTomato mouse with a monoclonal antibody against MYH11 (Invitrogen) and found that tdTomato-negative and Hu⁺ neurons exclusively immunolabeled for MYH11 (Suppl. Fig. 6a). We confirmed that a subpopulation of adult enteric neurons expressed Myh11 by immunolabeling small intestinal LM-MP tissues from a tamoxifen treated Myh11-cre^ERT2:Rosa26-YFP transgenic mouse line with anti-YFP/GFP and anti-Hu antibodies (Suppl. Fig. 6b). Anti-MYH11 antibodies labeled circular smooth muscle in the human duodenum, where it did not immunolabel nerve fibers and some submucosal neurons immunostained with antibodies against the pan-neuronal marker PGP9.5 but showed immunolabeling against PGP9.5-expressing myenteric neurons in the adult healthy small intestinal tissue (Supp. Fig 6c, d). Thus, we confirm that the MENs identified by us correspond to the cluster of ‘mesenchymal neurons’ observed by May-Zhang et al. As compared to our scRNAseq dataset, where MENs comprise of a significantly large numbers of cells, the MENs cluster in the May-Zhang et al. dataset is significantly smaller, possibly due to the fact that a large number of MENs were left unsequenced in their Phox2b-CFP^low fraction.

Fig 6:

Similarly, Drokhlyansky et al.⁴⁹ using MIRACL-seq found a small cluster of cells that express the MENs-signature genes (Cdh3, Dcn, Slpi, Aebp1, Wt1, Msln, Fmo2). However, the neuronal nature of these cells was not correctly identified as the study used a gene signature profile specific for neural crest-derived neurons (derived from their experiment with enriched Wnt1-cre2:Rosa26- tdTomato-expressing neural crest-derived cells). As a result, this cluster of cells – which shows similarity to the MENs transcriptomic profile, was instead annotated as Mesothelial. This ‘mesothelial’ cluster in the Drokhlyansky et al.⁴⁹ dataset shows similar gene expression profile as that of the cluster of ‘mesenchymal neurons’ from the May-Zhang et al.⁵⁰ dataset (Table 1), especially given that many of the marker genes for both clusters are similar (Lrrn4, Rspo1, Msln, Upk3b, Upk1b, Gpm6a, Wnt5a, Gpc3, Sulf1, Muc16). The ‘mesothelial’ cluster also expresses neuronally significant genes such as Gabra3 (GABA receptor A subunit 3⁷⁵), Prss12 (Neurotrypsin^{64, 76}), Synpr (Synaptoporin⁷⁷), Trpv4 (Transient receptor potential cation channel subfamily V member 4⁷⁸), and Nxph1 (Neurexophillin-1⁷⁹). Unfortunately, a deeper interrogation of this dataset was not possible due to the size and the manner in which the original data was processed⁴⁹. In addition, MENs were not represented in the snRNAseq carried out in Wright et al.²², which was based solely on neural crest-derived cells.

In contrast to these snRNAseq based studies, the Marklund Lab in two studies performed single cell RNA-sequencing (scRNAseq) on juvenile ENS cells, first on tdTomato⁺ cells from P21 Wnt1- cre:Rosa26:tdTomato lineage fate mapping mice in Zeisel et al.⁴⁸ and then from the newly characterized Baf53b-cre:Rosa26-tdTomato BAC-transgenic mouse line⁸⁰ in Morarach et al.⁵¹. Barring 2 small clusters (Clusters 5 and 11), the other neuronal clusters identified in Morarach et al.⁵¹ mapped to the neural crest-derived clusters identified earlier by them in Zeisel et al.⁴⁸. The gene expression signature of these two clusters (Cluster 5: Sst, Calb2; Cluster 11: Npy, Th, Dbh) did not match the signature of our MENs cluster or that of the cluster of mesenchymal neurons in the May-Zhang et al.⁵⁰ study (Table 1). We reasoned that the lack of MENs in the Morarach et al.⁵¹ data may be driven by the non-comprehensive nature of Baf53b-cre in this BAC transgenic mouse line labeling all neurons. On immunolabeling LM-MP tissues from two P21 Baf53b- cre:Rosa26-tdTomato mice with antibodies against Hu, we enumerated 1312 Hu⁺ cells in the ileum (759 from mouse 1 and 553 from mouse 2), 5758 Hu⁺ cells in the proximal colon (3411 from mouse 1 and 2347 from mouse 2), and 3352 Hu⁺ cells in the distal colon (1894 from mouse 1 and 1458 from mouse 2), and found that a significant population of myenteric neurons did not express tdTomato (Suppl. Fig. 6e, f). Given this non-comprehensive nature of the Baf53b-cre BAC transgenic line in labeling all myenteric neurons, and the congruence between the transcriptomic profiles of cells analyzed using the Baf53b-cre and the neural crest-specific Wnt1-cre lines, we infer that molecular taxonomy performed on enteric neurons in the Morarach et al.⁵¹ dataset is restricted to NENs.

In a recent report, Elmentaite et al.⁸¹ performed scRNAseq on dissociated single cells from the human gut at fetal, juvenile, and adult ages and catalogued multiple intestinal cell-types to generate a comprehensive gut cell atlas. Canonical enteric neurons, while adequately represented in the fetal ages, were not found in the juvenile and adult ages which were dominated by mesenchymal and other cell types⁸¹. Given that our scRNAseq approach of sequencing diverse cell types from the murine gut wall without any marker-based pre-selection was similar to this approach in human tissue, we reasoned that MENs would also be represented in their description of mesenchymal cell clusters, especially in the juvenile and adult ages. To detect putative MENs in this human gut cell atlas data, we projected the entire human mesenchymal cell scRNAseq dataset into the four murine MEN-specific NMF patterns and detected putative human MENs clusters (Suppl. Fig 4h). Since our identification of MENs and the generation of MEN-specific NMF patterns was based on post-natal healthy murine tissue, we next projected data from only post-natal healthy mesenchymal clusters into the four MEN-specific NMF patterns (Fig 6b). The four patterns identified a population of cells within the post-natal intestinal mesenchymal cell clusters suggesting that these cells are putative MENs in the human gut (Fig 6b). These putative MENs, which map to Cluster 1 in our representation of the intestinal mesenchymal cells from the post-natal gut cell atlas data (Fig 7a) were annotated previously by Elmentaite et al. as transitional stromal cells (Fig 7b). We next tested and found that the MENs marker Dcn (Fig 7c) along with pan-neuronal Uchl1 among other neuronal markers Hand2, Stx3, Slc17a9, Pde10a and Tubb2b (Fig 7d) that we earlier showed to be expressed by adult murine MENs, were expressed by the putative MENs in the adult human gut.

Fig 7.

We hypothesized that the reason these cells were not annotated as enteric neurons was due to the use of a NENs-restricted gene list. Apart from the NENs marker Ret, such gene lists often contain the pre-synaptic gene Snap25, a component of the SNARE complex, which is widely assumed to be a pan-neuronal marker in the adult ENS^{48, 49, 82, 83}. While we have previously observed the NENs-specific nature of RET, our transcriptomic data suggests that MENs either do not express Snap25 or do so in very small amounts (Suppl. Fig 4d). We tested the expression of SNAP-25 by NENs and MENs in the myenteric ganglia of an adult Wnt1-cre:Rosa26-tdTomato mouse and found that SNAP-25 expressing neurons expressed tdTomato, but not the MENs- marker MHCst, suggesting the NEN-specific expression of this canonical marker for synaptic neurons (Fig 8a). Next, we used a validated Snap25-Gcamp knock-in mouse line^84–86, where GCamp/GFP is knocked in at the Snap25 locus and hence is expressed by Snap25-expressing cells, to confirm that the expression of Snap25 is indeed restricted to a subset of myenteric neurons in the adult LM-MP layer. By using antibodies against Hu and against GFP, we found that the expression of Snap25-GFP is indeed restricted to a subset of adult myenteric neurons (Fig 8b). Finally, by immunolabeling adult small intestinal LM-MP of Snap25-Gcamp mice with antibodies against MHCst and against GFP, we confirmed the lack of Snap25 expression by a population of MENs (Fig 8c). The low or lack of expression of Ret, Snap-25, Elavl3 and Elavl4 in MENs as observed both by us and May-Zhang et al.⁵⁰ may help explain why these canonical gene lists missed the correct identification and annotation of MENs’ neuronal nature in human datasets⁸¹. Next, we examined the expression of MENs markers MHCst, MET, SLC17A9, and DECORIN in LM-MP tissues from adult humans with no known gut motility disorder and found them to be expressed by a population of myenteric neurons in normal adult human ENS (Fig 9; Suppl. Fig 7a, b, c). Thus, our bioinformatic approaches and immunofluorescence analyses of MEN-specific markers in adult human gut together shows that the human gut similarly contains a population of MENs.

Fig 8:

Fig 9:

We next tried sub-clustering the MENs cell cluster in our scRNAseq data to study whether individual subclusters could be discriminated on the basis of MENs markers or other genes. However, at the current sequencing depth, these clusters did not yield meaningful information on specific marker genes that could be used to define them (Suppl. Fig 8a, b).

The proportion of mesoderm-derived neurons expands with age to become the dominant population in the aging ENS

Since the ENS at birth consists solely of NENs¹³, we next studied the birth-date and eventual expansion of the MEN lineage in the post-natal gut. Using Wnt1- cre:Rosa26-tdTomato mouse line, we enumerated the tdTomato- MENs in LM-MP at different ages and found a significant age-associated increase in MENs proportions (Fig 10a, b, c; One- way ANOVA, F = 117.6, DFn, DFd = 5, 12; p < 0.0001). At P11, MENs were found only in few isolated myenteric ganglia (Fig 10a) and together represented only ∼5% of all myenteric neurons (tdTomato- neurons: 4.12% ± 2.98 SEM; n = 1327 Hu⁺ neurons counted from 3 mice), suggesting that they originate shortly before this age. The proportion of MENs in myenteric ganglia rises sharply thereafter: by P22, they account for ∼30% of all neurons (tdTomato- neurons: 29.63% ± 1.229 SEM; n = 742 neurons counted from 3 mice); and at P60 they represent ∼45% of all neurons (tdTomato- neurons: 46.38% ± 4.62 SEM; n = 700 neurons counted from 3 mice). During the adult ages, the proportions of MENs in the adult ENS remains relatively stable as at P180, MENs continue to represent roughly half of all myenteric neurons (tdTomato- neurons: 57.29% ± 3.62 SEM; n = 586 neurons counted from 3 mice). However, MENs dominate the proportions of the aging ENS as at the very old age of 17 months (P510), the ENS is populated almost exclusively by MENs (Fig 10b; tdTomato- neurons: 95.99% ± 1.62 SEM; n = 996 neurons counted from 3 mice). Thus, the proportions of NENs and MENs in the myenteric ganglia can be used as a biomarker for deducing ENS age, as a healthy ENS dominant in NENs would be juvenile, one with roughly equal proportions would be adult, and an aging ENS is dominated by MENs.

Fig 10:

We next tested whether the representation of MEN-specific transcriptomic signatures show a similar increase during the maturation and aging of the human gut. We tested whether the representation of the four MEN-specific NMF patterns (Fig 5a) increased with age in the human transcriptomic data from Elmentaite et al.⁸¹. We grouped the transcriptomic data from the wide- range of healthy post-natal specimens of the gut cell atlas into three groups: Juvenile (containing data on 2910 cells represented by 10 samples within the age range of 4 – 12 years), Adult (containing data on 3696 cells represented by 36 samples within the age range of 20 – 50 years), and Aging (containing data on 2848 cells represented by 37 samples within age range of 60 – 75 years) and tested whether the representation of the four MEN-specific patterns changed significantly between these age groups. Barring pattern 16 which showed a non-significant increase in projection weights between the Juvenile and Aging samples (One-way ANOVA, F = 2.10, DFn, DFd = 2, 8, p = 0.18), projection weights of MEN-specific pattern 27 (One-way ANOVA, F = 16.08, DFn, DFd = 2, 8, p = 0.0016), pattern 32 (One-way ANOVA, F = 8.13, DFn, DFd = 2, 8, p = 0.01) and pattern 41 (One-way ANOVA, F = 9.818, DFn, DFd = 2, 8, p = 0.007) all show significant age-associated increase in the human gut tissue (Fig 10d, every datapoint refers to mean projected pattern weight for cells within a defined age or age range). This suggests that analogous to the age-associated increase in MENs we observed in murine gut (Fig 10b), the proportions of MENs in the human gut also might increase with age. We also found that the representation of MENs cluster (Cluster 1) in the post-natal human gut cell atlas data expanded with maturation and age (Suppl. Fig 7d)

GDNF and HGF levels regulate the populations of the neural crest-derived and the mesoderm-derived neurons, respectively, in the maturing ENS

GDNF-RET signaling is responsible for proliferation and neurogenesis from NC-derived ENS precursor cells during development as well as for the survival of Ret-expressing enteric neurons^87–90. Similarly, HGF signaling has been shown to be essential for the proliferation of mesoderm-derived cells⁹¹. Since the expression of the HGF receptor MET and the GDNF receptor RET is exclusive to MENs and NENs respectively, we studied the correlation between age and levels of HGF and GDNF in LM- MP (Fig 10e, f). Given that both GDNF, as well as HGF expression is found in the LM-MP layer^{44, 92, 93}, we used LM-MP tissue to study how the levels of these two trophic factors change in and around the ENS tissue during maturation. We found that in agreement with a previous report⁹¹, GDNF protein levels are markedly reduced between the P10 and the P30 ages and remains reduced thereafter up to the P90 age (Fig 10g; n = 3 mice/age-group, F = 6.821, DFn, DFd = 1, 7, p = 0.0348, One way ANOVA). In addition, expression of Gdnf transcripts show significant reduction between the ages of P30 and P90 (Fig 10h; n = 6 mice/age-group, P30: 1.070 ± 0.179 SEM, P90: 0.539 ± 0.129 SEM, p = 0.037, Student’s t-test). On the other hand, HGF expression increases progressively between P10 and P90 ages (Fig 10i, j: n = 3 mice/age- group, Protein levels: F = 8.820, DFn, DFd = 1, 7, p = 0.02; mRNA levels: F = 36.98, DFn, DFd = 1, 16, p < 0.0001, One way ANOVA). The ratio of HGF to GDNF expression in ileal LM-MP tissue shows significant increase from the age of P10 to P90 (Suppl. Fig 9a, F = 48.60, DFn, DFd = 1, 7, p = 0.0002, One way ANOVA), and HGF expression is consistently higher than GDNF expression in the full thickness small intestinal tissue from adult and aging mice in the Tabula muris data⁹⁴ (Suppl. Fig 9b). We then queried the plasma proteome levels from the LonGenity cohort of 1,025 aging individuals and found that HGF levels correlated positively, while GDNF and RET levels correlated negatively with age (Suppl. Fig 9c)⁹⁵, suggesting parallels between our data from murine intestine and human plasma proteome data.

Since GDNF tissue levels are correlated with NENs proportions, we hypothesized that GDNF signaling regulates NENs proportions in maturing and adult ENS. On administration of GDNF or saline to cohorts of P10 Wnt1-cre:Rosa26-tdTomato mice over 10 days^{87, 96}, we found that GDNF treatment promoted the juvenile phenotype by enhancing the proportions of tdTomato+ NENs and correspondingly reduced the proportion of tdTomato- MENs in P20 mice to a level similar to that seen at the P10 age, while retaining the MENs proportions in saline-treated control mice remained at a level expected of its age (Fig 10k; Controls: %MENs: 25.87 ± 4.37 SEM of 1350 neurons from 3 mice; GDNF: %MENs: 3.86 ± 0.07 SEM of 1301 neurons from 3 mice; p = 0.0072; Student’s t-test). Consistent with a previous report⁸⁷, the GDNF-driven expansion of NENs and associated contraction of MENs conserved the total neuronal numbers (Fig 10k; Controls: neurons/ganglia = 20.60 ± 4.00 SEM; GDNF: neurons/ganglia = 24.39 ± 6.96 SEM; p = 0.52; Student’s t-test).

Since HGF tissue levels increase with age, we hypothesized that increasing HGF signaling drives the expansion of MENs in the maturing gut. HGF administration to cohorts of P10 Wnt1- cre:Rosa26-tdTomato mice over 10 days promoted an increase in the tdTomato- MENs population in P20 mice to levels previously observed in P60 mice, while tissues from saline - treated control mice exhibited a MENs:NENs ratio that is expected at P20 (Fig 10l; Controls: %MENs: 27.40 ± 5.49 SEM of 1970 neurons from 5 mice; HGF: %MENs: 49.37 ± 5.52 SEM of 1704 neurons from 5 mice; p = 0.02, Student’s t-test). Similar to GDNF treatment, HGF treatment also did not cause any significant change in total neuronal numbers (Fig 10l; Controls: neurons/ganglia = 20.35 ± 2.22 SEM; HGF: neurons/ganglia = 18.21 ± 0.69 SEM; p = 0.38, Student’s t-test).

Transcriptomic evidence of MENs genesis

As the proportions of MENs in the myenteric plexus rise significantly between the ages of P10 and P30, we reasoned that this phase was accompanied by significant neurogenesis of MENs. To provide transcriptomic evidence that MENs populations indeed expand during this phase, using 10X Genomics v3.1, we again performed unbiased scRNAseq and agnostic clustering, this time on unsorted cells from the myenteric plexus layer of two male mice of P21 age, when the MENs population is still expanding (Suppl. Fig 10a). Using NMF-generated MENs patterns to run projectR-based analyses (Suppl. Fig 10b) and by using the NENs and Neuroglia-specific markers, we again annotated the neuroglia, the NENs, and the MENs as before and found that at this age, we were able to sequence a similar number of MENs and NENs (∼500 cells) along with a large number (∼800) of neuroglia cells (Fig 11a). The expression of Met was again detected in the MENs cluster, as well as in a small subset of Ret-expressing NENs in the P21 cluster (Suppl. Fig 10c). MENs showed significantly higher UMI per cell, when compared to neural crest-derived cells (NENs and neuroglial cells) (Suppl. Table 4, One-way ANOVA; F = 187.4, DFn, DFd = 2, 1877, p < 0.0001), which we expected given the larger size of MENs⁹⁷ (Suppl. Fig 4a). By immunolabeling LM-MP tissue from 3 P21 Wnt1-cre:Rosa26-tdTomato male mice with antibodies against MET, we found that proportions of MET-expressing neurons at this age were significantly higher in the population of tdTomato-negative MENs than in the population of tdTomato⁺ NENs (of 152 MET- immunolabeled neurons across 3 mice, 81.51 ± 1.55 S.E.M. were tdTomato-negative MENs). The Met and Ret co-expressing NENs cluster in our data were described as populations of neural crest-derived Chat⁺ Calcb⁺ neurons by Zeisel et al.⁴⁸ at the P21 age (Suppl. Fig 10d), where they were annotated as ENT6, 8, and 9 clusters (and which correspond to clusters 1, 3, 4, 7, 6, and 7 of the Morarach et al.⁵¹ database). These Met- and Ret-co-expressing neurons were also described by Avetisyan et al.⁴⁴ as fetal-born neural crest-derived neurons that respond to both GDNF and HGF. By contrast, MENs do not express Ret and hence are the only cell population that can respond to solely to HGF. This data provides further evidence on the HGF-induced expansion of MENs during the juvenile phase.

Fig 11:

We next performed in silico analyses on the MENs cluster using our recently published Tricycle software⁹⁸, which is capable of inferring continuous cell-cycle position and can be applied to datasets with multiple cell types, across species and experimental conditions, including sparse data and shallow sequenced droplet-based dataset - thus allowing us to find evidence of cell cycling in the MENs subset. P21 MENs scRNAseq data was projected into universal cell cycle principal components (PCs) defined by the expression of 500 cell-cycle correlated genes and the continuous cell cycle position (theta) was measured as the angle from the origin (Fig 11b). Based on expression of 500 cell-cycle correlated genes, cells between 0.5π and 1.5π show hallmarks of being cycling cells (Fig 11b). Our analyses showed that a significant number of cells in the MENs cluster (213 out of 510) were present between 0.5 and 1.5π, which represents the mitotic position in the continuous cell-cycle phase (Fig 11c). We confirmed our analyses by observing that the expression of the cell cycle gene Top2a in the cells of the MENs cluster was highly correlated with this phase (Fig 11d). Next, we found that a key cell cycle regulator gene Ect2 (epithelial cell transforming 2), which encodes for the guanine nucleotide exchange factor protein ECT2, which activates RhoA in a narrow zone at the cell equator in anaphase during cell division⁹⁹, was expressed in the MENs cluster and its expression was highly correlated with the cell’s mitotic position on the Tricycle plot (Fig 11e). We then used antibodies against ECT2 to immunolabel the LM-MP from P21 mice and found expression of this key cell cycle regulator protein in DECORIN-expressing and MHCst-expressing MENs at this age (Fig 11f, g). Thus, the Tricycle-based computational analyses found evidence of cell cycling in the MENs cluster and identified a key cell cycle gene Ect2 as a marker for putative MENs precursor cells.

Reduced GDNF-RET signaling accelerates ENS aging to cause intestinal dysfunction

Since reduced GDNF or RET levels are associated with intestinal dysfunction in patients^{83, 100, 101}, we hypothesized that alterations in GDNF-RET signaling unrelated to those seen with normal aging, would cause dysfunction. To test this hypothesis, we studied lineage proportions and intestinal function in a mouse model of reduced RET signaling. Ret-null heterozygous mice, which have been previously used to study the effect of reduced RET signaling in the adult ENS, have normal ENS structure but altered gut physiology⁸⁷. A similar mouse model with a Ret^CFP allele has a CFP reporter inserted at its Ret locus rendering it null¹⁰². Ret+/CFP (or Ret+/-) mice carrying a single copy of the Ret gene showed significant reduction in Ret transcript and RET protein expression in the early post-natal murine gut¹⁰². Similarly, we found significantly reduced Ret transcript expression in the adult LM-MP of Ret^+/- mice compared to age-matched wildtype Ret^+/+ mice (Suppl. Fig 10e). In these mice, using antibodies against CFP/GFP, Hu, and MHCst, we confirmed that the NENs marker Ret-CFP, and the MENs marker MHCst were expressed by different neuronal subpopulations (Fig 12a). Using the adult Ret+/- mice, we tested the effect of partial Ret loss on ENS lineages at two adult ages: 9 weeks (∼P60) and 16 weeks (∼P110). Using antibody against CFP/GFP to detect CFP⁺ RET-expressing neurons, we found that Ret+/- mice show a significant reduction in the proportions of Ret-CFP+ NENs (Fig 12b; 9 weeks: %CFP⁺ neurons: 24.91 ± 5.42 SEM of 837 neurons from 3 mice; 16 weeks: %CFP⁺ neurons: 13.13 ± 0.98 SEM of 1227 neurons from 5 mice; p = 0.03, Student’s t-test). We observed a corresponding significant increase in the proportions of MHCst+ MENs with age in Ret+/- mice (Fig 12c; 9 weeks: %MENs: 58.74 ± 7.33 SEM of 644 neurons from 3 mice; 16 weeks: %MENs: 82.84 ± 3.58 SEM of 935 neurons from 5 mice, One-way ANOVA, p = 0.014), while control Ret+/+ mice showed no significant age-associated change in the proportions of MENs (Fig 12c; 9 weeks: %MENs: 43.27 ± 3.24 SEM of 780 neurons from 3 mice; 16 weeks: %MENs: 54.48 ± 4.07 SEM of 1276 neurons from 5 mice, One-way ANOVA, p = 0.36), which is consistent with our previous results that MENs proportions in wildtype animals are relatively stable in this time window. The expedited loss of NENs in Ret+/- mice confirms that depletion of endogenous RET signaling in the adult ENS accelerates the aging-associated loss of NENs.

Fig 12:

Having previously shown that aging mice have intestinal dysmotility¹⁰³, we tested whether the increased loss of NENs in the Ret+/- ENS, concomitant with the expansion of MENs accelerated ENS aging, causes an early onset of aging-associated intestinal dysmotility. We studied whole gut transit time (WGTT) in a cohort (n = 8) of adult Ret+/- mice and their littermate control (n = 10) Ret+/+ mice over 7 weeks, between 9 and 16 weeks of age. While 9-week adult Ret+/- mice were similar to control Ret+/+ mice, WGTT between the two genotypes diverged with age. Consistent with a prior report⁸⁷, 16-week old Ret+/- mice displayed significantly delayed intestinal transit compared to age-matched control Ret+/+ mice (Fig 12d; WGTT (in min) Ret+/+: 121.4 ± 4.01 SEM; Ret+/-: 157.3 ± 14.62 SEM, p = 0.048; Student’s t-test).

GDNF reverts aging in the ENS to normalize intestinal motility

Along with others, we have previously shown that aging is associated with slowing of intestinal motility^{104, 105}. We hypothesized that this may be a consequence of the replacement of the juvenile NENs population by MENs and therefore GDNF supplementation, by restoring a more balanced MENs:NENs ratio, may prevent age related changes in motility. We studied 17-month-old male mice (at an age where NENs constitute only ∼5% of all myenteric neurons; Fig 10c) before and after they received 10 days of intraperitoneal injection of GDNF or saline (n=5 for each cohort). While the two cohorts showed no significant difference in their intestinal transit times at baseline (WGTT (in min) Control: 192.8 ± 11.55 SEM; GDNF: 202.4 ± 7.60 SEM, p = 0.50, Student’s t-test), GDNF treatment caused significant improvement in intestinal transit (Fig 13a, WGTT (in min) Control: 175.0 ± 8.89 SEM; GDNF: 101.0 ± 8.91 SEM, p = 0.0004, Student’s t-test), reaching levels previously observed in healthy mice (Fig 12d).

Fig 13:

GDNF treatment significantly reduced the proportions of MHCst⁺ MENs (Fig 13b; Suppl. Fig 11a, b; Control: %MENs: 88.76 ± 1.48 SEM of 909 neurons from 5 mice; GDNF: %MENs: 74.12 ± 5.98 SEM of 799 neurons from 5 mice; p = 0.045; Student’s t-test) while increasing the numbers of RET⁺ NENs (Fig 13c; Suppl. Fig 11c, d; Control: RET⁺ neurons: 1.35 ± 0.05 SEM in forty-two 40X fields from 4 mice; GDNF: RET⁺ neurons: 3.19 ± 0.56 SEM in forty-one 40X fields from 4 mice, p = 0.017; Student’s t-test). Consistent with our earlier observation, GDNF treatment did not change the average neuronal numbers in the ileal myenteric ganglia (Control: 17.56 ± 1.82 SEM neurons per ganglia, GDNF: 16.01 ± 1.22 SEM neurons per ganglia, p = 0.49, Student’s t-test).

The two neuronal lineages in human disease

In addition to the detection of MENs in the adult healthy gut by bioinformatic and immunofluorescence means (Fig 4, 5, 6, 7), we next examined the relevance of the two neuronal lineages NENs and MENs in human disease. For this, we tested whether gut dysfunction associated with pathological reductions in NENs-signaling mechanisms are also associated with increased abundance of MENs. We mined previously generated and publicly available transcriptomic data from intestinal tissues from humans with normal gut motility (controls) and from patients with obstructed defecation (OD) associated with enteric neuropathy.¹⁰³. In this dataset, Kim et al. have previously shown that OD patients have significantly reduced intestinal expression of Gdnf and Ret¹⁰³. Again, using our transfer learning tool ProjectR^{72, 106}, which allowed us to query bulk-RNAseq data across species to estimate relative use of MENs-specific NMF-patterns generated earlier, we tested whether the usage of MENs-specific NMF-patterns was significantly altered between transcriptomes of OD patients and healthy humans (n = 3/group)¹⁰³. This projection analysis indicated that 2 of the 4 MENs-specific NMF-patterns showed significantly higher usage in OD samples compared to controls (Fig 14a, Suppl. Fig 12; Student’s t-test, p < 0.05), providing evidence that the enteric neuropathy in these patients is associated with a relative increase in MENs-specific transcriptional signatures. We further queried the differential gene expression analysis performed by Kim et al and found that in addition to Gdnf and Ret (as reported by them) other important NEN genes (Snap25 and Nos1) were also significantly downregulated while the MEN marker genes identified in this study (Clic3, Cdh3, Slc17a9) (Fig 14b) were significantly upregulated.

Fig 14:

Discussion

Current dogma states that the adult ENS is exclusively derived from neural crest precursors that populate the gut during fetal development¹⁰⁷. The results of this study indicate a much more complex system, one in which the fetal and early life “juvenile” ENS consisting of neural crest- derived enteric neurons (NENs) is incrementally replaced during maturation by mesoderm- derived enteric neurons (MENs). Eventually, the aging ENS consists almost exclusively of the neurons of the MEN lineage. Using a combination of floxed transgenic reporters bred into multiple lineage fate mapping transgenic mice, this study also provides the first definitive evidence of a significant mesodermal contribution to any neural tissue. Previously, the mesoderm and especially the Lateral Plate Mesoderm (LPM) was known to give rise to diverse cell-types within the viscera, including several in the gut¹⁰⁸, and our study shows that cells of this embryonic derivation also give rise to ENS neurons. A previous report on a dual origin of the ENS described adult enteric neurons from Foxa2+ or Sox17+ precursors¹⁴ and inferred that these were endodermal in origin. However, Foxa2 and Sox17 are also expressed by mesoderm-derived cells^25–28. By contrast, using two NC-lineage-traced mouse lines (Wnt1-cre and Pax3-cre), two lineage-traced mouse lines marking mesodermal derivatives (Tek-cre, and Mesp1-cre), and robust adult mesoderm markers, we identify a population of Mesp1-derived MHCst- and MET- expressing adult enteric neurons, which makes up the entire non-NC-derived population of myenteric neurons. These results provide robust evidence that the second source of adult enteric neurons is the mesoderm, and not the endoderm.

Distinct neuronal nature of MENs

While the developmental origins of neurons have been thought to be restricted to the neural tube and the neural crest, many neurons – such as placode- derived neurons, have been known to be derived from non-neural tube, non-neural crest developmental lineages¹⁰⁹. Similarly, while supposed pan-neuronal markers, such as Snap25, are expressed specifically by neurons, not all neurons express Snap25¹¹⁰. Indeed, release of neurotransmitters may occur independently of Snap25^111–116, suggesting that the expression of this protein is not central to neuronal functions. These present us with important context based on which, one can study the data generated using multiple lines of evidence to infer the neuronal nature of MENs. Firstly, we show that MENs are intra-ganglionic cells, which conditionally express reporters under mesoderm-specific lineage fate mapping models and are labeled with antibodies against the neuronal marker Hu. Secondly, MENs show expression of important ENS neurotransmitter-encoding or generating proteins CGRP, NOS1, and ChAT. Thirdly, with our single cell transcriptomic analysis, where we used the markers Calcb, Met, and Cdh3 that we had validated to identify and annotate the MEN cell cluster, we discovered an additional set of neuronal-specific marker genes (Pde10a, Vsnl1, Stmn2, Stx3, Tubb2b, Slc17a9, Hand2, Gpr88, Elavl2, Ntf3) which were found to be expressed widely or selectively within MENs. These markers include both novel as well as established enteric neuronal markers^{2, 63, 117}. Fourthly, further transcriptomic analyses of the MENs scRNAseq cluster identifies the expression of many other neuronally significant genes – including those encoding for ion channels, neurotransmitter receptors, hormones, calcium channels, etc. – which provide a further glimpse into the diverse neuronal functions of MENs. Finally, our analyses shows that MENs express the pan-ENS- expressed gene Phox2b, and that our MENs scRNAseq cluster maps to a neuronal cluster in the May-Zhang et al.⁵⁰ database, which we shows expression of ‘mesenchymal’ and ‘neuronal’ markers⁵⁰. Our data, together with our analyses of publicly available data, thus provides evidence on the neuronal nature of MENs. While the co-expression of aforementioned pan-neuronal markers establishes the neuronal nature of MENs, the expression of set of key genes in MENs suggests putative specialized functional roles of these subpopulations when they co-exist with NENs in adults. These include the expression of genes associated with putative sensory (Calcb)¹¹⁸, neuromodulatory (Ntf3)¹¹⁹, secretory (Cftr)¹²⁰, immunomodulatory (Il18)², hormone regulating (Ghrh)¹²¹, extracellular matrix (ECM)-modulating (Dcn), and motility-regulating roles (Slc17a9 or VNUT)¹²². These MEN-enriched or MEN-specific expression profiles suggest that MENs perform various roles in the regulation of diverse gastrointestinal functions. While the current scRNAseq data is restricted in providing meaningful cell clusters that can inform the true functional diversity of MENs, future work will provide detailed analyses and functional characterization of MENs.

Beyond establishing the neuronal nature of MENs, our analyses also help establish the distinct neuronal nature of MENs. Based on the lineage fate mapping experiments with multiple combinations of neural crest and mesoderm lineage fate mapping mice models and based on the expression of mesoderm-specific markers such as MHCst and Decorin, we identified that MENs are a population of mesoderm-derived cells. MHCst is a subset of myosin heavy chain isoforms that has been previously shown to be specifically expressed by a subset of mesoderm-derived muscle cells¹²³. While prior reports do not clarify the nature of the exact epitope and the gene responsible for MHCst, the expression of the smooth muscle myosin heavy chain protein MYH11 specifically in MENs suggests that the epitope detected by the anti-MHCst antibody in MENs is of the mesoderm specific protein MYH11 protein.

Further, our bioinformatic analyses of the human gut cell atlas data generated by Elmentaite et al.⁸¹ also shows the existence of a putative MENs-like cell cluster in the human gut, which shows a similar co-expression of mesenchymal and neuronal markers as observed by us in our MENs cluster and by May-Zhang et al.⁵⁰ in their ‘mesenchymal neuronal’ cluster. Indeed, Elmentaite et al.⁸¹ identified this supposed mesenchymal cluster by the expression of the mesoderm-specific gene Decorin. Further, by observing the co-expression of neuronal and mesoderm-specific markers in human gut tissue, we validate the presence of neurons that express mesoderm- specific markers or MENs in adult human tissue.

Detection and identification of MENs in other datasets

Recent studies have significantly increased our understanding of the molecular nature of various classes of enteric neurons. However, these studies differ significantly in their methods and tools used, which have precluded them from detecting or correctly identifying the developmental origins of MENs. While May-Zhang- et al.⁵⁰ observed the presence of a small neuronal cluster that expressed mesenchymal marker genes, their study did not further interrogate the developmental origins of the ‘mesenchymal’ enteric neurons. Drokhlyansky et al.⁴⁹ in their unbiased MIRACL-seq analyses and Elmentaite et al.⁸¹ in their human gut cell atlas datasets detected the MENs cluster but could not correctly identify their neuronal nature, possibly owing to the expression of mesenchymal genes in this cluster, and also due to the absent or lower than expected expression of canonical neural markers such as Ret and Snap25. In contrast to these studies that performed single cell or nuclear transcriptomic studies that were unbiased by lineage markers, other transcriptomic datasets from Drokhlyansky et al.⁴⁹, Wright et al.²², and Zeisel et al.⁴⁸ were generated using enriched neural crest-derived cells, which contributed to the absence of MENs in these datasets. An important outlier to these studies is Morarach et al.⁵¹ from the Marklund group, which utilized the newly characterized BAC-transgenic line Baf53b-cre to label and flow sort enteric neuronal cells for subsequent single cell RNA sequencing. Our analysis of ileal and colonic tissue from this line shows that it does not comprehensively label all enteric neurons. While this may seem at odds from the Morarach et al.⁵¹ report, this apparent difference in the behavior of this transgenic line can be explained by the fact that in contrast to the established and accepted practice of counting all Hu-immunolabeled neurons employed by us and others^{22, 50, 124}, the Morarach et al. study restricted the enumeration to only include neurons that showed only cytoplasmic and not nuclear expression of Hu. Neuronal Hu proteins (HuB, HuC, and HuD) show both nuclear and cytoplasmic localization¹²⁵, and since the Morarach et al.⁵¹ study only enumerated cytoplasmic Hu-expressing neurons, their characterization of this transgenic mouse line was performed using only a subset of all neurons. In addition, compared to Drokhlyansky et al.⁴⁹, Wright et al.²², and May-Zhang-et al.⁵⁰ whose protocols for neuronal enrichment were based on flow sorting neuronal nuclei, the Marklund group (in both Zeisel et al.⁴⁸ and Morarach et al.⁵¹) performed flow sorting of late post- natal enteric neuronal cells, which can cause significant stress on neuronal cells. These factors, together with the BAC transgenic nature of the Baf53b-cre line, may have contributed to the lack of MENs in the Morarach et al.⁵¹ data. Further, all neuronal cells described in the Morarach et al.⁵¹ study express Snap25 (which we show is not expressed by all adult enteric neurons), and mostly map to known neural crest-derived cell clusters identified in Zeisel et al.⁴⁸. This suggests that the rich transcriptomic atlas of enteric neurons generated by deep transcriptomic profiling of Baf53b- cre-derived neurons carried out by Morarach et al.⁵¹ is restricted only to the NENs, which comprise of half of all adult enteric neurons. Thus, prior experimental data could not establish the true developmental identity, correct annotation, and transcriptomic definition of a population of mesoderm-derived neurons or MENs in murine and human tissue.

By contrast, we were able to correctly identify the MENs cluster in unbiased scRNAseq of adult murine LM-MP cells and use it to both identify and validate the expression of MENs-specific markers in the adult murine LM-MP tissue, as well as use bioinformatic analyses to identify the MENs in publicly available murine data from May-Zhang et al.⁵⁰ and human data from Elmentaite et al.⁸¹. This allowed us to show the expression of MENs markers MHCst, MET, DCN, and SLC17A9 in many enteric neurons in adult humans, suggesting that the mesoderm-derived ENS may be a feature common to mice and humans alike. MENs in both species can readily be identified by their expression of MENs-specific markers, thus providing a convenient tool to further identify and study these neurons.

Development of MENs and MENgenesis

The presence of small population of MENs in the gut at the post-natal age day 10 (P10) and their expansion during maturation suggests that key timepoints in their development correspond with the post-natal development of the ENS. Importantly, since MENs populations are expanding at this age, performing scRNAseq-based bioinformatic analyses of the juvenile LM-MP cells not only provided us with the computational evidence of significant cell cycling in the MENs cluster at this age, but also provided us with the information on Ect2, a cell-cycle marker expressed by cycling putative MENs precursor cells. Ect2 transcripts were found expressed by cycling cells within the MENs cluster at the P21 age, which also showed significant enrichment in the expression of the cell-cycle gene Top2a. Except for a cluster of cells annotated as proliferating glial cells that show expression of the cell cycle gene Top2a, Ect2 gene expression is not detected in other murine myenteric neural crest-derived cells in the scRNAseq study performed by Zeisel et al.⁴⁸ at the P21 juvenile age. Our anatomical observation on the expansion of the MEN population in the adolescent gut along with the computational evidence of cell cycling in the MENs cluster at this age, and the detection and validation of Ect2 gene expression in a subset of MENs provide evidence of significant MENgenesis during the juvenile phase. ECT2 expression was found in MHCst- and DCN- expressing myenteric cells at this age, suggesting that many of these are proliferating MEN- precursor cells, and may not be terminally differentiated functional neurons.

The rapidly expanding MEN population is 10-11 days old in the P21 post-natal gut (as they originate in the gut just before P10 post-natal age) and thus, the age at which these cells partake in regulation of diverse gastrointestinal functions is yet unknown. This important issue precludes us from using the two MENs transcriptomic datasets to compare how MENs functions change with age. Prior studies have reported on the differing birthdates of various neural crest-derived enteric neuronal sub-populations in the developing fetal gut^{9, 51, 126}. Similar studies need to be performed in the future to understand the birthdates of the various subpopulations that comprise MENs, and to ascertain when do they become functional. A recent study by Parathan et al.¹²⁴ highlighted important neuronal and neurochemical changes in the adolescent gut. Further studies on how the ENS matures and ages thus need to be performed in the context of NENs and MENs to truly understand how shared and distinct functions of these two lineages change with maturation and age. A crucial aspect of this work will require further identifying the molecular markers for the Ect2-expressing proliferating MEN progenitor cells. Our earlier work showed that Nestin-expressing neural progenitor cells are neural crest-derived cells that contribute to steady state neurogenesis in the adult ENS¹⁹, thus suggesting that MENs are derived from Ect2- expressing non-Nestin⁺ precursor cells.

Use of scRNAseq to identify MENs

The use of single cell and single nucleus transcriptomic approaches provides an important tool that has been recently used to perform molecular taxonomy of fetal and post-natal enteric neurons and glial cells^{22, 48, 50, 51}. In addition to performing deep sequencing that allows for cluster-based molecular taxonomy, these tools can also be utilized for identification of novel cell types, where prior information on their presence and marker expression is unknown^{52, 98, 127–129}. Compared to the other studies that utilized this tool to study the ENS, we used the single cell transcriptomics to both confirm the presence of MENs and identify more MEN-specific markers; as well as to use the MEN-specific transcriptomic signature to query their presence in other murine and human datasets. In contrast to its use for performing molecular taxonomy, which requires deep sequencing of an enriched cell population, our strategy of opting for a lower read depth (shallow sequencing) is efficient, practical, and economical for discovering and identifying novel celltypes⁵². In addition, even if a gene is being expressed at a medium - low level, there is a certain probability that it will not be detected by scRNAseq in both deep and shallow sequencing approaches. This is important since the older v2.0 10X Genomics chemistry (used for our 6-month dataset) captures only 10-14% of the cellular transcriptome, and that from the newer v3.1 chemistry (for our P21 data) captures only a third of the transcriptome from a suspension of homogenous single cells^{130, 131}. While a partial solution to this issue is to increase sequencing depth, beyond a certain point, this strategy leads to diminishing returns as the fraction of PCR duplicates increases with deeper sequencing. While getting insights into cellular heterogeneity of MENs is important for furthering our understanding of how they regulate diverse gut functions, the transgenic mice that can specifically label only MENs or all enteric neurons to allow for their enrichment are currently not available^{22, 48–51, 71}. Without the availability of MENs- specific transgenic models, the scRNAseq strategy employed by us is the only currently available tool for querying the transcriptome of MENs. These reflect significant bottlenecks that have disallowed the comprehensive and an in-depth characterization of all enteric neurons in all the current studies on the adult ENS.

Expression of MET and RET by MENs and NENs

Since our scRNAseq data highlighted the lineage-specific nature of the expression of Ret and Met in the adult ENS, we studied whether these genes regulated the origin, expansion, and maintenance of the neuronal populations that express them. The expression of MET by MENs in the mature adult ENS is consistent with an earlier study by Avetisyan et al. that reported that a deletion of Met in the Wnt1-cre expressing NC-lineage did not alter the abundance of MET⁺ neurons in the adult ENS⁴⁴. If Met was expressed solely by neural crest-derived cells, then the loss of MET by enteric neurons would have been expected since the conditional deletion of Met gene in this floxed Met transgenic animal has been observed to cause a significant or comprehensive loss of MET protein expression in other organs^132–135. While at first glance, observations of the mutually exclusive nature of RET and MET provided by us and by Avetisyan et al.⁴⁴ in the adult gut might not agree with Zeisel et al.’s⁴⁸ observed co-expression of Ret and Met transcripts in NENs at the P21 developmental stage, our observation that a small proportion of NENs in the P21 gut express MET helps provide an explanation. Our P21 data and the data from Zeisel et al.⁴⁸ are consistent with one of the observations made by Avetisyan et al.⁴⁴, who showed the existence of a population of fetal neural crest-derived enteric neurons that are co-dependent on GDNF and HGF and hence would express both RET and MET. In contrast to these observations made in the fetal and juvenile ENS, the anatomical observations made by Avetisyan et al.⁴⁴ in the adult ENS (where MET immunoreactivity in the adult enteric neurons is exclusive of Ret expression) suggests that the population of fetal-derived Ret and Met co-expressing neurons is lost after maturation. Compared to the P21 scRNAseq data from Zeisel et al.⁴⁸ that shows that all Met-expressing neural crest- derived neurons co-express Calcb (CGRP), Ret, and Chat, Avetisyan et al.’s anatomical data in the adult gut shows that despite accounting for a third of all enteric neurons and expressing CGRP, MET-expressing adult enteric neurons do not co-express Ret and only a small proportion of these MET⁺ neurons co-express ChAT. Thus, Avetisyan et al.’s data agrees with our data on the adult ENS which shows exclusive expression of RET and MET by NENs and MENs, respectively, and that while MENs consist of a large proportion of CGRP-expressing MET⁺ neurons, they contain only a small population of ChAT-expressing neurons. Along with Avetisyan et al.’s data, our data suggests that the adult small intestinal myenteric plexus contains a significant proportion of neurons that do not immunolabel with antibodies against NOS1 or ChAT. This observation is consistent with earlier data from Parathan et al.¹²⁴ who showed that ∼35-40% of small intestinal myenteric neurons from 6-week-old mice did not immunolabel with antibodies against NOS1 or ChAT.

Since they do not express the neural crest-marker RET, MENs are the only neuronal population that respond to HGF alone. In this study, we show that MENs are dependent on HGF-MET signaling in a manner analogous to the requirement of GDNF-RET signaling by NENs. While our study does not discount that the exogenous HGF supplementation in juvenile mice may alter proportions of RET and MET-co-expressing neurons within the NENs subpopulation, our report focuses on HGF’s ability to drive the expansion of MENs population in the ENS. Similarly, Ret haploinsufficiency-mediated loss of the NEN lineage causes a proportional increase in the MEN lineage, while maintaining overall neuronal numbers. Conservation of neuronal number, despite significant shifts in the proportions of the two neuronal lineages suggests the presence of yet unknown quorum sensing signaling mechanisms between the two lineages that reciprocally regulate their populations to maintain the structure of the post-natal ENS. This again implies the existence of a yet undefined precursor cell responsible for the expansion of the MEN population in a manner analogous to what we have previously described for NENs¹⁹.

Relevance to aging and human disease

While it is known that Gdnf expression in the mature gut is significantly downregulated^{44, 87, 96, 136}, the functional consequences of this loss have been unclear⁸⁷. We found that reduced GDNF-RET signaling drives the age-dependent loss of NENs as this loss can be slowed or reversed by GDNF supplementation and accelerated by further reduction of Ret expression. In aging animals, GDNF-driven resuscitation of NENs was associated with a functional normalization of age-associated delay in intestinal transit. Our results identify a novel role for GDNF in maintaining or resuscitating the canonical NEN population, while preserving overall enteric neuronal numbers. In the last few years, studies have focused on identifying juvenile protective factors (JPFs), the loss of which correlates with or drives maturation- and aging-associated phenotypes¹³⁷. In this context, GDNF may therefore qualify as a JPF or a senolytic as its presence maintains the dominance of NENs in maturing ENS, corresponding to a juvenile phenotype; and its re-introduction promotes and resuscitates the genesis of NENs in adult and aging gut to restore healthy gut function. The exact nature of the cells that respond to GDNF re-introduction and generate NENs is yet unknown, but it can be hypothesized that these may include Nestin⁺ enteric neural stem cells and/or GDNF-responsive Schwann cells^{19, 138}. In addition to its effect on resuscitating NENs, acute administration of GDNF also augments peristaltic motility by increasing neuronal activity^{22, 139}. Since the aging gut has a very small number of RET-expressing GDNF-responsive NENs, we assume that GDNF-driven normalization of intestinal motility would first increase the proportions of RET-expressing NENs and then augment their activation. This further emphasizes the importance of using GDNF or similar RET agonists in the treatment of age-associated intestinal dysmotility.

While Avetisyan et al. showed that MET expression highly correlated with CGRP-expression in the adult murine myenteric plexus⁴⁴, MET co-expression by all MHCst⁺ human enteric neurons in our data and by all PGP9.5⁺ neurons in human enteric neurons by Avetisyan et al.⁴⁴ suggests that MET-expression in the human ENS may not be restricted to CGRP-expressing IPANs. In addition, recent data shows that our understanding of the nature of markers that can correctly identify IPANs in the murine and human ENS is still evolving^{50, 140}, suggesting that ascribing function to neurons using a single marker may be incorrect.

The detection of MENs cluster in the post-natal and adult scRNAseq data in the gut cell atlas generated by Elmentaite et al.⁸¹, their increasing representation with age, and the subsequent validation of the MEN-specific markers MHCst, MET and SLC17A9 in a subset of adult human ENS neurons shows that our findings may be of clinical importance. Many gastrointestinal motility functions are compromised with advancing age, resulting in clinically significant disorders¹⁴¹ Although the exact mechanisms will need to be worked out, our results indicate that a MENs- predominant ENS is associated with significant differences in gut motility. The progressive nature of the change in enteric neuronal lineages with age may have pathological implications for the elderly gut when the balance is overwhelmingly in favor of MENs. Understanding the forces that regulate parity between these two different sources of neurogenesis therefore holds the promise of arresting or reverting progressive loss of gut motility with increasing age. Our results may also have implications for the pathogenesis of disordered motility unrelated to aging, as downregulation of Gdnf and Ret expression has been associated with diverticular disease and obstructed defecation in adult patients^{100, 103, 142}. GWAS analyses further showed that Gdnf gene was in the eQTL associated with increased incidence of diverticular disease¹⁴³. It is in this context that the identification of SNAP-25 as a NENs marker gene and the transcriptomic pattern analyses of patients with chronic obstructed defecation are significant. Expression of Snap25 is upregulated by GDNF and is significantly downregulated in gut tissues of patients with diseases associated with significant reduction in GDNF-RET signaling (diverticular disease and obstructed defecation)^{83, 100, 103}. Thus, while earlier thought to be a pan-neuronal^{49, 144}, establishing the identity of SNAP-25 as a NEN lineage-restricted marker provides us not only with an important tool to query proportions of NENs in murine and human tissue but also suggests the NENs-limited nature of prior observations based on Snap25 transgenic animals^{82, 144}. Gut dysmotility disorders that present with conserved ENS neuronal numbers¹⁰⁰ have puzzled investigators on the etiology of the disease. Our results suggest an explanation based on NENs-to-MENs lineage shifts that conserve overall enteric neuronal numbers with variable effects on function.

In conclusion, using multiple lines of evidence, we show that the ENS of the juvenile and mature mammalian gut is dominated by two different classes of neurons that have distinct ontogenies and trophic factors, associated with functional differences. The shift in neuronal lineage may represent a functional adaptation to changes in nutrient handling, microbiota or other physiological factors associated with different ages. Further research in the functional differences between the neuronal lineages and factors that regulate the parity between these two nervous systems in humans during the lifetime will be important to advance our understanding of the adult ENS and the treatment of age-related and other pathological disorders of gut motility.

Acknowledgements

We would like to thank Dr. Chulan Kwon (JHU) for his kind gift of the Mesp1-cre mice, Dr. Jeremy Nathans (JHU) for his kind gift of the Tek-EGFP, Tek-cre and Hprt-tdTomato lineage-traced mice and his support, and Dr. Vanda Lennon (Mayo Clinic) for her kind gift of the ANNA1 anti-Hu antisera. We thank Dr. Akira Sawa (JHU) for his help with the Rosa26-EGFP mouse line, and Dr Meera Murgai (NCI) for her help with the Myh11-cre^ERT2:Rosa26-YFP mouse tissues. We thank Dr. Xinzhong Dong (JHU) and Dr. Mark Donowitz (JHU) for their support. The microscopy was performed on the Ross Imaging Core at the Hopkins Conte Digestive Disease Center at the Johns Hopkins University (P30DK089502) using the Olympus FV 3000rs (procured with the NIH-NIDDK S10 OD025244 grant) and we thank Dr. George McNamara for his help with confocal training and imaging. The 10X Genomics Chromium processing for scRNAseq was performed at the GRCF Core and the sequencing was performed at the CIDR core at the Johns Hopkins University. S.K. was funded through a grant from the Ludwig Foundation, a grant from the NIA (R01AG066768), a pilot award from the Hopkins Digestive Diseases Basic & Translational Research Core Center grant (P30DK089502), and a pilot award from the Diacomp initiative through Augusta University. L.A.G was funded in part by a Johns Hopkins Catalyst Award and a grant from the NIA (R01AG066768). J.S. was funded through the Maryland Genetics, Epidemiology, and Medicine training program sponsored by the Burroughs Welcome Fund. P.J.P. was funded through the Hopkins Conte Digestive Disease Center at the Johns Hopkins University (P30DK089502), NIDDK (R01DK080920), the Maryland Stem Cell Research Foundation (MSCRF130005), and a grant from the AMOS family.

Author contributions

S.K., and P.J.P. conceived the study; S.K., M.A., M.L, L.A.G. and P.J.P. designed the research studies; S.K. L.B, C.Z, A.S., A.B., and M.S. conducted the experiments; S.K., Z.H., M.K., J.L., L.B., S.B., E.V., J.S., S.N., M.S., G.L., and L.A.G. acquired and analyzed the data; S.K., and P.J.P drafted the manuscript.

Methods

Animals

Experimental protocols were approved by The Johns Hopkins University’s Animal Care and Use Committee in accordance with the guidelines provided by the National Institutes of Health. Presence of vaginal plug was ascertained as 0.5 days post-fertilization and this metric was used to calculate age of mice. Only male mice were used for the studies detailed in this report. The Wnt1-cre:Rosa26-tdTomato lineage-traced line was generated as detailed before by breeding the B6.Cg-Tg(Wnt1-cre) with the Ai14 transgenic mouse line (Jax #: 007914) containing the Rosa26- tdTomato transgene^{19, 145, 146}. The Wnt1-cre line used in this study is the well validated line for studying neural crest derivatives previously used by us and others^{18, 19, 22, 48, 147}, and not the Wnt1- cre2 line used by Drokhlyansky et al.⁴⁹. We again validated the fidelity of this line by observing whether any aberrant recombination occurred in our Wnt1-cre:Rosa26-tdTomato line by studying the intestinal mucosal tissue and found that the line behaved expectedly to label neural crest- derivatives and epithelial cells (Suppl. Fig 1i). Pax3-cre:Rosa26-tdTomato lineage-traced line was generated by breeding the Ai9 transgenic mouse line (Jax #: 007909) with the Pax3-cre transgenic mouse (Jax #: 005549). The Wnt1-cre:Hprt-tdTomato mouse was generated by breeding our aforementioned Wnt1-cre transgenic mouse line with the Hprt-tdTomato transgenic mouse line (Jax #: 021428, kind gift from Prof. Jeremy Nathans). Mesp1-cre:Rosa26-tdTomato mice were generated by breeding the Mesp1-cre transgenic mice¹⁴⁸ (gift from Dr. Chulan Kwon, JHU) with the Ai14 transgenic mice. Mesp1-cre:Rosa26-GFP mice were generated by breeding the Mesp1-cre transgenic mice with the Rosa26-EGFP transgenic mouse line (gift from Dr. Akira Sawa, JHU). Neither the Rosa26-tdTomato (Ai14) line nor the Rosa26-EGFP line showed any reporter expression in the absence of a cre driver (Suppl. Fig 1j, k). Ret+/CFP mice (MGI:3777556) were inter-bred to get a colony of adult Ret+/+ and Ret+/CFP mice. RetCFP/CFP mice died at or before term. Tek-cre:Hprt-tdTomato mice were generated by breeding Tek-cre transgenic mice (also known as Tie2-cre; Jax #: 004128) with the Hprt-tdTomato transgenic mouse line. 17-month- old male C57BL/6 mice from the aging mouse colony of the National Institute of Aging were procured for the GDNF-treatment experiment.

Human Tissues

Human tissues were obtained under IRB protocol IRB00181108 that was approved by Institutional Review Board at the Johns Hopkins University. Pathological normal specimens of human duodenum and colon were obtained post-resection. Tissues were obtained from adult donors and were de-identified such that the exact age, gender, and ethnicity of the donors was unknown.

Tissue preparation

Mice were anesthetized with isoflurane and sacrificed by cervical dislocation. A laparotomy was performed, and the ileum was removed and lavaged with PBS containing penicillin-streptomycin (PS; Invitrogen), then cut into 1-cm-long segments and placed over a sterile plastic rod. A superficial longitudinal incision was made along the serosal surface and the LM-MP was peeled off from the underlying tissue using a wet sterile cotton swab ¹⁹ and placed in Opti-MEM medium (Invitrogen) containing Pen-Strep (Invitrogen). The tissue was then laid flat and fixed within 30 minutes of isolation with freshly prepared ice cold 4% paraformaldehyde (PFA) solution for 45 minutes in the dark to preserve fluorescence intensity and prevent photo-bleaching. All LM-MP tissues post-isolation were fixed within 30 minutes of their isolation. After the fixation, the tissue was removed and stored in ice cold sterile PBS with Pen-Strep for immunofluorescence staining and subsequent microscopy.

For human tissues, duodenal tissue from adult human patients (n = 3 patients), who did not have any prior history of chronic intestinal dysmotility, that was removed by Whipple procedure was obtained. A colonic sample from a pathologically normal colonic resection from an adult donor suffering from colon carcinoma who similarly did not have prior history of chronic intestinal dysmotility was also obtained. The resected tissue was placed in ice cold Opti-MEM medium (Invitrogen) containing Pen-Strep (Invitrogen). The mucosal and sub-mucosal tissue was dissected out in the medium under light microscope and the muscularis layer containing myenteric plexus tissue was obtained. The tissue was laid out between two glass slides and fixed overnight in ice cold 4% PFA after which it was removed and stored in ice cold sterile PBS with Pen-Strep for immunofluorescence staining, optical clarification, and subsequent microscopy.

Immunohistochemistry

For murine tissue: The fixed LM-MP tissue was washed twice in ice-cold PBS in the dark at 16°C. The tissue was then incubated in blocking-permeabilizing buffer (BPB; 5% normal goat serum with 0.3% Triton-X) for 1 hour. While staining for antibodies that were mouse monoclonal, 5% normal mouse serum was added to the BPB. The tissue was then removed from the BPB and was incubated with the appropriate primary antibody at the listed concentration (Supplementary Table 1) for 48 h at 16°C in the dark with shaking at 55 rpm. Following incubation with primary antibody, the tissue was washed three times (15-min wash each) in PBS at room temperature in the dark. The tissue was then incubated in the appropriate secondary antibody at room temperature for 1 hour while on a rotary shaker (65 rpm). The tissue was again washed three times in PBS at room temperature, counterstained with DAPI to stain the nuclei, overlaid with Prolong Antifade Gold mounting medium, cover-slipped, and imaged.

Immunostaining for ChAT was performed using fixed and unpeeled murine tissues from adult male Wnt1-cre:Rosa26-tdTomato mice. The intestinal segments were obtained by perfuse fixing the mouse with ice-cold PBS and then freshly prepared 4% PFA solution. Once the intestine was harvested, the ileum was flushed with cold HBSS with 0.4 mol/L N-acetyl-l-cysteine and then PBS to remove luminal contents. After, the tissue was fixed in 4% PFA overnight. The fixed tissue was then immersed in 2% Triton-X 100 solution for 2 d at 15°C for permeabilization. The tissue was immunostained with antibodies against Hu (ANNA1) and against ChAT (rabbit, Abnova). The primary antibody was then diluted (1:100) in the dilution buffer (0.25% Triton X-100, 1% normal goat serum, and 0.02% sodium azide in PBS) and incubated for 1 day at 15°C. An Alexa Fluor 647-conjugated goat anti-rabbit secondary antibody (1:200, Invitrogen) and an Alexa Fluor 488- conjugated goat anti-human secondary antibody (1:200, Invitrogen) were then used to reveal the immunopositive structure. Finally, the labeled specimens were immersed in FocusClear solution (CelExplorer) for optical clearing before being imaged by laser confocal microscopy using a Zeiss LSM 510 META. The LSM 510 software and the Avizo 6.2 image reconstruction software (VSG) were used for the 2D and 3D projection of the confocal micrographs. The Avizo software was operated under a Dell T7500 workstation and the “Gaussian Filter” function of Avizo was used for noise reduction of the micrographs.

Colchicine treatment: For CGRP immunostaining, mice were injected with Colchicine at a concentration of 5 mg/Kg body weight 16 hours (overnight) before they were sacrificed. The mice were housed singly during this time and adequate gel packs were provided. Food and water were provided ad libitum. On the following day, the mice were sacrificed, and their LM-MP tissues were harvested as detailed above.

For human tissue: The fixed muscularis layer containing myenteric plexus tissue was removed from ice cold PBS and incubated in blocking-permeabilizing buffer (BPB; 5% normal goat serum, 5% normal mouse serum with 0.3% Triton-X) for 4 hours. The tissue was then removed from the BPB and was incubated with the appropriate primary antibody at the listed concentration (Supplementary Table 1) for 5 days at 16°C in the dark with shaking at 55 rpm. Following incubation with primary antibody, the tissue was washed five times (15-min wash each) in PBS at room temperature in the dark. The tissue was then incubated in the appropriate secondary antibody at 16°C in the dark with shaking at 55 rpm for 2 days. The tissue was again washed in dark for five times in PBS that contained DAPI at room temperature. After the final wash, the tissue was suspended in tissue clarification buffer CUBIC¹⁴⁹ for 1 hour at 4°C in the dark after which it was overlaid with Prolong Antifade Gold mounting medium, cover-slipped, and imaged. Briefly, the CUBIC optical clarification buffer was made by mixing 2.5 g of urea (25% by wt), 2.5 g of N, N, N ’, N ’-tetrakis (2-hydroxy-propyl) ethylenediamine (25% by wt), 1.5 g of Triton X-100 (15% by wt) in 35 ml of Distilled Water. The solution was shaken till the ingredients were dissolved and yielded a clear viscous solution.

Microscopy

Except for ChAT microscopy, all other imaging was done by using the oil immersion 63X objective on the Leica SP8 confocal microscope and by using the oil immersion 40X objective on the Olympus Fluoview 3000rs confocal microscope with resonance scanning mode. For thick tissues, such as human tissues, the Galvano mode of the Olympus Fluoview 3000rs microscope that enabled higher resolution imaging and averaging was used. Images obtained were then analyzed using Fiji (https://fiji.sc/).

Live tissue imaging of the mouse gut tissue was performed by harvesting small intestinal tissue from an adult Wnt1-cre:Rosa26-tdTomato mouse and immediately putting it in OptiMEM solution. The tissue was then immediately put on a chamber slide containing OptiMEM and imaged in live tissue culture conditions of the EVOS M7000 microscope under the 20X objective.

Enumeration of neurons

Enumeration of tdTomato⁺ and tdTomato- neurons was performed on tissue with native tdTomato fluorescence. The native tdTomato fluorescence in fixed tissue was intense, thus there was no need to increase fluorescent signal using anti-RFP/tdTomato antibodies.

Enumeration of enteric neurons to study alterations to ENS neuronal numbers was performed by following the well-established protocol^150–153 of counting the numbers of neurons in myenteric ganglia and comparing them between animals for statistically significant differences. Identification of myenteric ganglia was performed according to our pre-determined method published earlier¹⁹. Briefly, contiguous clusters of neurons were defined as a ganglion and the total numbers of neurons within these clusters were enumerated as numbers of myenteric neurons per ganglion. As a rule, clusters of 3 neurons or more were deemed to consist of a ganglion and our enumeration strategy did not count extra-ganglionic neurons. We imaged ∼10 ganglia per tissue for our enumeration and each group studied had n ≥ 3 mice. Identification and enumeration of neurons and detection of co-localization was performed manually by trained laboratory personnel. Fig 1d comprises of the data from 6 P60 mice, which includes data from the 3 mice that were used for the development experiment presented in Fig 8c.

Protein isolation and detection

After the LM-MP tissue was isolated, it was weighed and placed in a sterile 1.5 ml microfuge tube. 1X RIPA buffer (Cell Signaling Technology) with Halt Protease Inhibitor Cocktail (Thermo Scientific) at 5X concentration, Phosphatase Inhibitor Cocktails II and III (Sigma-Aldrich) at 2X concentrations were added to the tissue lysate buffer. Tissue was disrupted using 1.0 mm silica beads in Bullet Blender 24 (Next Advance) for 5 mins at highest setting. The lysate was incubated at 4°C with shaking for 30 mins, centrifuged at 14,000rpm for 20 mins and the supernatant was taken and stored in −80°C in aliquots. Protein concentration was estimated using Bradford assay solution (Biorad) following the manufacturer’s protocol. For immunoblotting, 40 µg of protein was loaded per well of 4%-20% gradient denaturing gel (Biorad). Protein marker used was Precision Plus Dual Color standards (Biorad). After fractionating the proteins, they were blotted onto ImmunBlot PVDF membrane (Biorad) overnight at 4°C at 17V for 12-16 hours. After blotting, membrane was blocked with Odyssey TBS blocking buffer (Li-Cor) for 1 hour at room temperature with shaking. Incubation with primary antibodies were carried out at 4°C with shaking for 24 hours. Following binding, the blot was washed 4 times with TBS-T (Tris Buffered Saline with 0.5% Tween) for 15 mins each with shaking at room temperature. Secondary antibody incubation was carried out in dark at room temperature for 1.5 hours with shaking. The blot was then washed 4 times for 15 mins each and imaged on Odyssey CLx system (Li-Cor).

Part of the aforementioned protocol was also followed for detecting HuB protein (recombinant Human ELAVL2, 1-346aa, from Biosource; Catalogue number MBS205995). The protein was loaded in two amounts (0.5 µg and 1.25 µg) in separate lanes of a 4-20% gradient denaturing gel (Biorad) along with protein marker Precision Plus Dual Color standards (Biorad). After fractionating the proteins, they were blotted onto ImmunBlot PVDF membrane (Biorad) overnight at 4°C at 17V for 12-16 hours. After blotting, membrane was blocked with Odyssey TBS blocking buffer (Li-Cor) for 1 hour at room temperature with shaking. Incubation with primary antibodies were carried out at 4°C with shaking for 24 hours. Following binding, the blot was washed 4 times with TBS-T (Tris Buffered Saline with 0.5% Tween) for 15 mins each with shak-ing at room temperature. Secondary antibody incubation was carried out in dark at room tempera-ture for 1.5 hours with shaking. The blot was then washed 4 times for 15 mins each and imaged on Odyssey CLx system (Li-Cor).

Antibody used to detect RET was a well validated antibody used previously for detecting RET in ENS¹⁵⁴ and was further validated by us using a western blot on total proteins isolated from murine small intestinal LM-MP (Suppl Figure 8e). Antibodies used are detailed in the Supplementary Table 1.

RNA isolation and quantitative detection of specific transcripts

The isolated tissue was stored in RNALater Solution (Ambion). RNA was isolated using RNeasy Mini Kit (Qiagen) following manufacturer’s protocol. RNA quantification was carried out using Epoch Microplate Spectrophotometer (BioTek). cDNA synthesis was carried by SuperScript IV VILO Master Mix (Invitrogen). Quantitative Real-time PCR was carried out using Taqman Gene Expression Master Mix (Applied Biosystems) and Roto-Gene Q (Qiagen). The probes used are listed in Supplementary Table 1.

Single Cell RNA sequencing and analyses of LM-MP tissues

Single cell preparation from murine ileal LM-MP tissues: Ileal LM-MP tissues from male littermate C57/BL6 wildtype mice were isolated by peeling as previously described¹⁹. Succinctly, mice were anesthetized with isoflurane and sacrificed by cervical dislocation. A laparotomy was performed and the small intestine was removed and lavaged with PBS containing penicillin-streptomycin (PS; Invitrogen). The ileum was then isolated which was then cut into 2-cm-long segments and placed over a sterile plastic rod. A superficial longitudinal incision was made along the serosal surface and the LM-MP was peeled off from the underlying tissue using a wet sterile cotton swab and placed in Opti-MEM medium (Invitrogen) containing PenStrep (Invitrogen). The tissues were then dissociated in Digestion Buffer containing 1 mg/ml Liberase (Sigma-Aldrich) in OptiMEM. Tissues from mouse 1 were dissociated in the Digestion buffer containing Liberase TH and tissues from mouse 2 were dissociated in the Digestion buffer containing Liberase TL. Dissociation was performed at 37°C for 30 minutes on a rotary shaker, after which the cells were centrifuged at 200g for 7 minutes, and the pellet was resuspended in ice cold sterile PBS. The cell suspension was passed through a 40 µm cell sieve and the resulting filtered cell suspension was again centrifuged at 200g for 7 minutes. This process of cell centrifugation and filtration was repeated two more times, after which the cells were resuspended in 1 ml ice cold sterile PBS. The repeated steps of serial cell washes and filtration removed clumps and debris and the viability of the resulting cell suspension was estimated to be >90% using Trypan Blue dye test. The cells were then processed through 10X Genomics Chromium V2.0 system (for the 6-month-old murine tissues) and 10X Genomics Chromium V3.1 system (for the P21 murine tissues) according to the manufacturer’s suggested workflow. The processing was done at the GRCF Core Facility at the Johns Hopkins University. The pooled libraries were sequenced on an Illumina HiSeq 2500 (for 6-month-old tissues) and Illumina Novaseq (for P21 tissues) to an average depth of 3.125×10⁸ reads per sample library. The sequencing was performed at the CIDR core facility at the Johns Hopkins University.

For processing data from scRNAseq on 6-month-old murine LM-MP tissues: Pre-processing of FASTQs to Expression Matrices: FASTQ sequence files were processed following a Kallisto Bustools workflow compatible with downstream RNA velocity¹⁵⁵. References required for pseudo- alignment of transcripts were obtained using the get_velocity_files (functionality of BUSpaRSE (https://github.com/BUStools/BUSpaRse)), with “L=98” for 10X Genomics v2.0 sequencing chemistry. Reads were pseudo-aligned to an index built from Ensembl 97 transcriptome annotation (Gencode vM22; GRCm38). Across two samples processed, a total of 578,529,125 reads were successfully pseudo-aligned. Barcodes within a Hamming distance of one to known 10X Genomics v2.0 barcodes were corrected. Reads were classified as “spliced” or “unspliced” by their complement to the target list of intronic sequences and exonic sequences, respectively, and subsequently quantified separately into expression count matrices. Spliced counts are used for all analyses.

Single cell gene expression analysis: scRNA-seq count matrices were analyzed using Monocle3. 11,123 high-quality cells were identified as meeting a 200 UMI minimum threshold with a mitochondrial read ratio of less than 20%; droplets that did not meet these criteria were excluded from the analysis. Mitochondrial counts were determined as the sum of reads mapping to 37 genes annotated to the mitochondrial genome. All genes with non-zero expression were included for downstream analysis. Raw counts were first scaled by a cell-size scaling factor and subsequently log10 transformed with a pseudo-count of 1. Normalized values are used in place of raw counts unless otherwise noted.

Prior to UMAP dimensionality reduction, batch effects between the two biological replicates were assessed and corrected via the mutual nearest neighbors (MNN) algorithm as implemented by Batchelor in Monocle3¹⁵⁶. 15 clusters of cells in the UMAP embedding were identified by Leiden community detection. 30 marker genes for each cluster were identified based on greatest pseudo R² values and used for supervised annotation of cell types by searching UniProt, Allen Cell Atlas and through literature search with Pubmed.

NC-derived cell clusters were identified by expression of NC markers Ret and Sox10. MEN cluster was identified by its expression of CGRP-coding Calcb, Met, and Cdh3. The pan-MENs protein marker MHCst was identified by labeling with an antibody S46 which labels all members of the MHCst family⁴¹. Since the antibody does not identify a single gene product, MHCst immunostaining could not be used to identify a specific gene marker for use in the annotation of the MEN cluster. For further analysis into the MENs population, the full LM-MP dataset was subset to include only the 2,223 cells annotated as MENs. These cells were re-processed as above, but with a reduced PCA dimensionality of k=20 as input for the UMAP embedding. 5 clusters of cells in the UMAP embedding were identified by Leiden community detection (k= 10, resolution = 5e - 4).

P20 LM-MP scRNA processing and analysis: Pre-processing of FASTQs to Expression Matrices: FASTQ sequence files were processed following a Kallisto Bustools workflow compatible with downstream RNA velocity. References required for pseudo-alignment of transcripts were obtained using the get_velocity_files (functionality of BUSpaRSE (https://github.com/BUStools/BUSpaRse)), with “L=91” for 10X Genomics v3.1 sequencing chemistry. Reads were processed following the same steps as the six month LMMP dataset. Spliced counts are used for all analyses.

Single cell gene expression analysis: scRNA-seq count matrices were analyzed using Monocle3. 11,264 high-quality cells were identified as meeting a 600 UMI minimum threshold with a mitochondrial read ratio of less than 20%; droplets that did not meet these criteria were excluded from the analysis. Mitochondrial counts were determined as the sum of reads mapping to 37 genes annotated to the mitochondrial genome. All genes with non-zero expression were included for downstream analysis. Raw counts were first scaled by a cell-size scaling factor and subsequently log10 transformed with a pseudo-count of 1. Normalized values are used in place of raw counts unless otherwise noted.

NC-derived and MENs cell clusters were identified by expression of established marker genes. Further, the identity of MENs was confirmed by a projection analysis using MENs-specific patterns learned in the six month dataset. To investigate differences in division potential within the ENS, the data was subset to NENs, neuroglia, and MENs for computational cell cycle analysis (NENs: n = 526; neuroglia: n = 844; MENs: n = 510). Continuous scores for cell cycle position (0-2π) were calculated via Tricycle by projecting each cell into the default reference space. The relevance of these scores was confirmed by profiling the scores’ correlation with expression of genes with known variation over the cell cycle. Lastly, each cell was unambiguously characterized as “non- cycling” (0-0.5π or 1.5-2π) or “cycling” (0.5-1.5π) by binarizing Tricycle scores.

External datasets: Spliced and unspliced count matrices for mesenchymal subset of the Gut Cell Atlas were obtained from https://www.gutcellatlas.org/81, and existing annotations were used when available. SnRNAseq counts matrices from May-Zhang et al.⁵⁰ were downloaded from GEO GSE153192. Bulk RNAseq dataset from the Human Obstructed defecation study by Kim et al.¹⁰³ was obtained from GEO GSE101968.

Pattern discovery and ProjectR analyses: Pattern discovery was utilized to identify sets of co-expressed genes that define cell-type specific transcriptional signatures. The normalized expression matrix was decomposed via non-negative matrix factorization (NMF) as implemented in the R package NNLM, with k=50 and default parameters. Cell weights for each pattern were grouped by assigned cell-type and represented by heatmap. Pattern vectors were hierarchically clustered by a Euclidian distance metric, implemented in ComplexHeatmap¹⁵⁷. These patterns were then tested on the bulk RNA-Seq expression matrix for the Human Obstructed Defecation study¹⁰³. The log2 expression (log2(rpkm + 1)) from this study was projected onto the NMF patterns using projectR⁷². Students’ t tests were performed on the projection weights from the Control and OD groups to test for differences between them. For the projection of datasets generated by May-Zhang et al.⁵⁰ and Elmentaite et al.⁸¹, log10-transformed normalized counts were projected.

Analysis of May-Zhang et al. dataset

Raw UMI count matrices for the single-nucleus RNA sequencing of enteric neurons performed by May-Zhang et al.⁵⁰ from the ileum of adult mice were obtained from GEO (GSE153192). The DropletUtils package (v1.10.3)^{158, 159} was used to compute barcode ranks, knee, and inflection point statistics for all the samples. All barcodes in the 10X genomics libraries with total UMI greater than the inflection point were considered to be cell-containing droplets and retained for further analysis. The inDrop libraries were filtered to retain all barcodes with UMI greater than 200. The filtered matrices were combined to form a SingleCellExperiment object.

Cell-specific scaling factors accounting for composition biases were computed using computeSumFactors (scran v1.18.7)¹⁶⁰. The raw UMI counts were then scaled using the scaling factors and log 2 transformed after addition of a pseudo count of 1 using logNormCounts (scater v1.18.6)¹⁶¹. Batch effects between the samples were corrected using fastMNN (batchelor v1.6.3)¹⁵⁶.

Unsupervised clustering was performed using a graph-based clustering technique. The graph was obtained using buildSNNgraph (scran v1.18.7)¹⁶⁰ and the louvain algorithm implemented in the igraph package (v1.2.7)¹⁶² was used to identify communities. Differential expression testing was performed between all pairs of clusters using findMarkers (scran v1.18.7)¹⁶⁰. One-sided Welch’s t-test was used to identify genes that were upregulated in a cluster as compared to any others in the dataset. Benjamini-Hochberg method was applied (to combined p-values) to correct for multiple testing and the genes were ranked by significance. The top ranked genes were used as marker genes for the clusters. The log fold change from the pairwise comparison with the lowest p-value was used as the summary log fold change.

The log normalized counts from this dataset were then projected into the NMF patterns learned on our scRNAseq data. This was performed using transfer learning implemented in projectR (v1.6.0). The projection weights corresponding to the four MEN-specific patterns were then clustered using graph-based clustering and visualized by UMAP (uwot v0.1.10)( https://github.com/jlmelville/uwot).

Whole gut transit time analyses

Whole-gut transit time (WGTT) for every mouse was analyzed by the method using the carmine red protocol¹⁹. Mice were placed in individual cages and deprived of food for 1 hour before receiving 0.3 mL 6% (wt/vol) carmine solution in 0.5% methylcellulose by oral gavage into the animal’s stomach. The time taken for each mouse to produce a red fecal pellet after the administration of carmine dye was recorded in minutes. The experiment was terminated at 210 minutes post-gavage and the WGTT of any mice that did not expel the red dye at the termination was marked at the value of 210 min. The mean difference in whole gut transit time (in minutes) between both the Ret+/+ and Ret+/- mice cohorts, and the GDNF and Saline-treated Control cohorts were analyzed statistically.

In vivo Injections

GDNF injection: Similar to prior report that gave sub-cutaneous injections of GDNF to post-natal mice⁹⁶, we took 6 littermate 10 day old (P10) male Wnt1-cre:Rosa26-tdTomato mice and divided into two subgroups, GDNF and Control. Each mouse in the GDNF group was injected sub-cutaneously with 50 µl of 2mg/ml of GDNF (Peprotech Catalogue #: 450-44) every other day, while the Control group was injected with 50 µl of sterile saline. The mice were given 5 doses and then sacrificed on P20, after which their LM-MP tissues were isolated as detailed above. The tissues were then immunostained with antibodies against Hu and imaged. In a separate experiment, adult (P60) mice were also injected sub-cutaneously with GDNF (100 µl of 100 µg/ml of GDNF). The mice were given 5 doses over a course of 10 days and then sacrificed on P70, after which their LM-MP tissues were isolated as detailed above. For studying the effect of GDNF on aging mice, two cohorts of 17-month-old male C57BL/6 mice (n=5 mice/cohort) were obtained from the aging colony of the National Institute of Aging. Before the start of dosing, the whole gut transit time was assayed. The animals were then injected daily sub-cutaneously either with 100 µl of saline (Control) or 100 µl of GDNF (500 µg/ml) for 10 consecutive days, after which the mice were sacrificed, and their LM-MP tissues were isolated as detailed above. The tissues were then immunostained with antibodies against Hu and imaged.

HGF injection: Similar to prior report that gave sub-cutaneous injections of HGF to post-natal mice, we took 6 littermate 10 day old (P10) male Wnt1-cre:Rosa26-tdTomato mice and divided into two subgroups, HGF and Control. Each mouse in the HGF group was injected sub-cutaneously with 100 µl of 2mg/ml of HGF (Peprotech Catalogue #:315-23) every other day, while the saline group was injected with 100 µl of sterile saline. The mice were given 5 doses and then sacrificed on P20, after which their LM-MP tissues were isolated as detailed above. The tissues were then immunostained with antibodies against Hu and imaged.

Statistics

Data was analyzed using Graphpad Prism 8.3.1 and R using Unpaired Students t-test, Simple Linear Regression, and Ordinary One-Way ANOVA.

Data

All raw data are provided in Supplementary Table 2 and 3. We imaged atleast 10 ganglia per tissue for our enumeration and each group studied had n ≥ 3 mice. Raw single cell RNA sequencing data is archived on the NCBI GEO server and can be accessed under the accession numbers GSE156146 and GSE213604.

Code

Code generated for the scRNA/snRNA analyses is available at https://github.com/jaredslosberg/6month_ENS/.

Acknowledgements

We wish to thank Prof. Lior Pachter at Caltech for his help with responding to the reviewer’s comments.

Suppl Figure 1: Presence and absence of tdTomato reporter expression in the myenteric ganglia of adult neural crest lineage-traced mice

(a) Hu-immunostained Wnt1-cre:Rosa26-tdTomato LM-MP tissue, where tdTomato-expression (red) marks neural crest-derived cells, and Hu (green) labels myenteric neurons, we observed the presence of higher intensity fluorescence (red arrows) and lower intensity fluorescence (orange arrows) in various tdTomato+ neurons. The myenteric ganglia also contained Hu+ neurons that did not express tdTomato (green arrow) and these neurons alone were deemed to be non-neural crest-derived neurons. Scale bar denotes 10 µm.

(b) Reporter non-expressing Hu-immunolabeled neurons (green, green arrows) can be observed both through the expression of the tdTomato reporter (red), as well as on antibody labeling of the tdTomato reporter (gray).

(c) Freshly harvested small intestinal tissue from Wnt1-cre:Rosa26-tdTomato mice when imaged under 20X magnification of a live tissue fluorescence microscope shows the presence of cells that express tdTomato (red) and those that don’t.

(d) Cells showing hyper-aggregation of the tdTomato reporter can be found in freshly harvested full-thickness small intestinal myenteric plexus in tissues imaged under 20X magnification of a live tissue fluorescence microscope.

(e) Myenteric ganglia from fixed LM-MP tissue of adult Wnt1-cre:Rosa26-tdTomato mice show show hyper-aggregation of tdTomato (white arrow) when imaged using a confocal microscope. Nuclei are stained with nuclear marker DAPI (blue). Scale bar denotes 10 µm.

(f) Cells showing hyper-aggregation of tdTomato (red, red arrow) show no detectable nuclear labeling with DAPI (blue) or of pan-neuronal marker Hu (green, green arrow). Scale bar denotes 10 µm.

In addition, interrogation of two other neural crest-specific lineage-traced mouse lines confirms our observations with Wnt1-cre:Rosa26-tdTomato mice. (g) In Wnt1-cre:Hprt-tdTomato adult male mice, tdTomato expression (red) is present only in a subset of Hu-immunolabeled (green) neurons (red arrows) while it is not detected in other Hu-expressing neurons (green). Similarly,

(h) only a subpopulation of Hu+ adult enteric neurons (green) is derived from the neural crest, as seen by its co-expression of tdTomato (red, red arrows) in the Pax3-cre:Rosa26-tdTomato mice while other neurons that do not express tdTomato (white arrows) are inferred to be not derived from the neural crest. Scale bar denotes 10 µm.

(i) Orthogonal views of small intestinal villus from a Wnt1-cre:Rosa26-tdTomato mouse that shows presence of tdTomato-expressing (red) neural crest-derived fibers inside the villus but a complete lack of tdTomato-expression by intestinal epithelial cells. Nuclei are labeled with DAPI (blue). Scale bar denotes 10 µm.

(j) Representational image of the myenteric plexus from a Rosa26-tdTomato mouse without any cre driver showing that in absence of cre recombinase to drive tdTomato expression, there is a lack of tdTomato expression in this layer. Nuclei are labeled with DAPI (blue). Scale bar denotes 10 µm.

(k) Representational image of the myenteric plexus from a Rosa26-EGFP mouse without any cre driver showing that in absence of cre recombinase to drive EGFP expression, there is a lack of EGFP expression in this layer. Nuclei are labeled with DAPI (blue). Scale bar denotes 10 µm.

(l) Western blot image showing detection of two different loading amounts of the recombinant HuB protein by the rabbit polyclonal antibody against HuC/D (abcam). Molecular Weight marker with defined band sizes indicated on the right.

Suppl Figure 2: Tek-expressing and Mesp1-expressing mesodermal lineage contributes to adult enteric neurons.

(a) Adult small intestinal LM-MP layer from a male Tek-EGFP mouse, where EGFP expression (green) labels Tek-expressing cells was immunostained with specific antibodies against the pan-neuronal marker PGP9.5 (red) and imaged under a confocal microscope. The Tek-EGFP expression did not label any cells within the myenteric plexus marked by PGP9.5. Nuclei are labeled with DAPI (blue). Scale bar denotes 10 µm.

(b) We observed the presence of higher intensity fluorescence (red arrows) and lower intensity fluorescence (green arrow) in various tdTomato+ cells in a section of heart tissue from a 6-month-old Mesp1-cre:Rosa26-tdTomato adult male mouse, where tdTomato (red) labels Mesp1-derived cells of the mesodermal lineage that includes cardiomyocytes and vasculature. Nuclei are stained by DAPI (blue). Scale bar denotes 10 µm.

(c) Adult small intestinal LM-MP layer from a 6-month-old male Mesp1-cre:Rosa26-tdTomato mouse, where tdTomato (red) cells are mesoderm-derived cells, was immunostained with specific antibodies against the pan-neuronal marker Hu (green) and imaged under a confocal microscope. tdTomato-expressing cells outside the ganglia (red arrows) were brighter than tdTomato+ neurons inside the ganglia (green arrows). Nuclei are labeled with DAPI (blue). Scale bar denotes 10 µm.

(d) Adult small intestinal LM-MP layer from a 2-month old male Mesp1-cre:Rosa26-EGFP mouse, where EGFP cells are mesoderm-derived, was immunostained with specific antibodies against the reporter EGFP (green) and the pan-neuronal marker Hu (red) and imaged under a confocal microscope. A population of Hu-immunolabeled neurons were also labeled with antibodies against EGFP (green arrows) while other neurons were not (red arrows), suggesting the derivation of EGFP⁺ neurons from the Mesp1-expressing mesoderm. (e) No primary control for the Anti-GFP immunostaining showed very weak EGFP expression in the fixed LM-MP tissue from these animals. Nuclei were labeled with DAPI (blue). Scale bar = 10 µm. Consistent with the observations made with the tdTomato reporter in this transgenic line, the Mesp1-cre:Rosa26- EGFP animals also show a brighter reporter signal in the vasculature and a dimmer reporter signal in the myenteric neurons.

(f) MHCst+ (cyan) labels all tdTomato+ (red) Hu+ (green) adult neurons (white arrows) but not the tdTomato- neurons (yellow arrow) in the mesoderm-specific Tek-cre:Hprt-tdTomato lineage-traced lines. Nuclei are stained with DAPI (blue). Scale bar = 10 µm.

Suppl Figure 3: Validation of MENs-specific markers by immunohistochemistry and confocal microscopy

Small intestinal LM-MP from several adult 5-month-old male Wnt1-cre:Rosa26-tdTomato mice, where tdTomato-expression (red) labels neural crest-derived cells, was used to ascertain and validate the expression of the novel MENs marker and of known neuronal marker genes.

(a) The myenteric ganglion shown in Fig 2h is presented here as a merged image (where MHCst immunostaining and tdTomato-expression are shown together) and with MHCst immunostaining without tdTomato signal to better appreciate the presence of MHCst-expressing neurons in the ganglia. MHCst-specific antibody S46 (green) was found to label tdTomato non-expressing cells within the myenteric ganglia as well as cells outside of the myenteric ganglia. Neural crest-derived cells within the myenteric ganglia expressed tdTomato. Nuclei are labeled with DAPI (blue). Scale bar denotes 10 µm.

The immunostaining of the MENs-markers elucidated through scRNAseq analyses: (b) CDH3, (c) AEBP1, (d) SLPI, (e) CFTR, (f) NT-3, (g) CLIC3, (h) FMO2, (i) IL-18, (j) MYL-7, and (k) SLC17A9 (all immunolabeled cells in these panels are green; green arrows) was found to be localized to or enriched in cells without tdTomato expression, while the tdTomato+ cells were found to show not to express these MENs-specific markers (red arrow). In addition, we also immunolabeled the myenteric ganglia from this transgenic mouse line with known neuronal markers (gray) and with antibodies against the pan-neuronal marker Hu (green) to show that the neuronal markers Elavl2 or HuB (l), GPR88 (m), HAND2 (n), PDE10A (o), TUBB2B (p) and VSNL1 (q) were expressed by both tdTomato⁺ (red arrow) and tdTomato^- Hu⁺ neurons. Orthogonal views of images where 3D image z-stacks were obtained are provided for better visualization of marker expression in the 3D space. Nuclei in these tissues were labeled with DAPI (blue). Scale bar denotes 10 µm.

Suppl Figure 4: scRNAseq metrics

(a) MENs are significantly larger in size as compared to the NENs, as the mean Feret’s Diameter of the Wnt1-cre:Rosa26-tdTomato- MENs is significantly more than that of the Wnt1- cre:tdTomato+ NENs. Student’s t-test; ** p value < 0.01.

(b, c) UMAP representation of sequenced cells generated from 10X Genomics Chromium-derived libraries from the two mouse samples, taken from separate littermate adult male mice, shows similar representation of their cells across the various clusters.

(d) Quasiviolin plot of expressed genes provide a comparison of important MEN-specific and neuronal genes between the NENs and MENs clusters in the adult ileal May-Zhang et al. dataset and the NENs, Neuroglial, and MENs clusters in our scRNAseq dataset.

(e) Top genes that are pattern drivers for the four MEN-specific NMF patterns along with their gene weights.

(f) The four MEN-specific patterns were together used to label a group of cells (True) that together showed increased usage of the MEN-specific patterns, when compared to other cell types (False) in the adult ileal snRNAseq dataset from May-Zhang et al.⁵⁰.

(g) Top-expressed genes in the group of cells (f) that showed higher utilization of the MEN-specific patterns, compared to other cells in the dataset.

(h) Representation of the four MEN-specific NMF patterns in the UMAP of our data and the projection of the May-Zhang et al.⁵⁰ murine ileal dataset and the entire gut mesenchymal single cell RNA seq dataset by Elmentaite et al.⁸¹ into the MEN-specific patterns highlights cell clusters that show higher usage of the four NMF patterns. Bottom row shows increasing usage of these patterns with age in the entire mesenchymal gut cell atlas data.

Suppl Figure 5: Expression of neuronally significant genes by MENs

Sparklines plot of gene expression in various cell clusters, as revealed by single cell transcriptomics of diverse LM-MP cells from small intestinal tissue of two 6-month-old mice, show that expression of some genes is enriched in the MENs cluster, when compared to the other major clusters of cells in the dataset. Lower to higher expression in the plot is represented by yellow to blue gradient. Subsequent analyses shows that proteins coded by these genes are neuronally significant.

Suppl Figure 6: Expression of mesodermal marker MYH11, and reporter expression under the Baf53b-cre transgenic line in post-natal myenteric neurons

(a) Orthogonal views of a myenteric ganglion of an adult Wnt1-cre:Rosa26-tdTomato mouse, where tdTomato (red) is expressed by neural crest-derived cells, shows that immunolabeling with a MYH11-specific antibody (green) labels tdTomato non-expressing cells. Nuclei are labeled with DAPI (blue). Scale bar denotes 10 µm.

(b) Representative image of the myenteric plexus layer from a Myh11-cre^ERT2:Rosa26-YFP adult male mouse, which was dosed with tamoxifen for two consecutive weeks, and which when immunostained with antibodies against YFP/GFP (green) and the pan-neuronal marker Hu (red), shows the presence of YFP-expressing neurons (green arrows) in the myenteric ganglia that also contains YFP negative neurons (red arrow). Nuclei are labeled with DAPI (blue). Scale bar denotes 10 µm.

(c) Immunostaining of healthy human adult duodenal tissue with antibodies against the pan-neuronal marker PGP9.5 (green) and MYH11 (red), shows that in the circular muscularis layer of the tissue, MYH11 labels the smooth muscle and not the PGP9.5 immunostained nerve fibers and cell bodies. Scale bar denotes 10 µm.

(d) Immunostaining of healthy human adult duodenal tissue with antibodies against the pan-neuronal marker PGP9.5 (green) and MYH11 (red), shows that in the myenteric ganglia, MYH11 immunolabels a subset of PGP9.5-expressing myenteric neurons (white arrows). Scale bar denotes 10 µm.

(e) Representative image of a myenteric ganglion from a Baf53b-cre:Rosa26-tdTomato adult male mouse, where tdTomato (red) is expressed by Baf53b expressing cells or those derived from Baf53-expressing cells, when immunostained with antibodies against the pan-neuronal marker Hu (green), shows the presence of neurons that express tdTomato (white arrows) along with neurons that do not express tdTomato (green arrows). Nuclei are labeled with DAPI (blue). Scale bar denotes 10 µm.

(f) Graphical representation of the proportions of tdTomato-expressing Hu⁺ myenteric neurons in the LM-MP tissue of two male P21 Baf53b-cre:Rosa26-tdTomato mice shows that myenteric ganglia from ileum, proximal colon, and distal colon contain significant proportions of neurons that do not express tdTomato under the Baf53b-cre transgenic line.

Suppl Figure 7: Putative mesoderm-derived enteric neurons (MENs) are present in adult human myenteric ganglia.

(a, b) The mesodermal specific markers, MET and MHCst are expressed by a subpopulation of human enteric neurons. Myenteric ganglia from normal duodenal tissue of 2 different adult human subjects express the pan-neuronal marker Hu (green), with distinct subpopulations positive and negative for MET (blue) and MHCst (red). MET and MHCst-expressing neurons (red arrows) in the human ENS are assumed to be mesoderm-derived enteric neurons (MENs) and those not expressing these MENs markers (green arrows) are presumably NENs. Nuclei are labeled with DAPI (gray). Scale bar = 10 µm.

(c) The mesodermal specific marker DECORIN is expressed by a subpopulation of human enteric neurons. Myenteric ganglia from normal human duodenal tissue when immunostained with antibodies against the pan-neuronal marker Hu (green) and DECORIN (red) shows the presence of DECORIN-expressing myenteric neurons, which are assumed to be MENs. Nuclei are labeled with DAPI (blue). Scale bar = 10 µm.

(d) Distribution of the various mesenchymal clusters of the post-natal healthy gut cell atlas data (generated by Elmentaite et al.⁸¹) by donor age. For the age of each donor, the proportion of cells annotated as belonging to each cluster is shown. Proportions of Cluster 1 are observed to expand with age.

Suppl Fig 8: Sub-clustering MENs in the scRNAseq dataset generated from 6-month-old mice.

(a) Sub-clustering the MENs cluster from our scRNAseq dataset into 5 distinct clusters (b) was found to be not based on the expression of cluster-specific marker genes.

Suppl Fig 9: HGF and GDNF ratios in maturing and adult murine gut.

(a) Ratio of HGF:GDNF values in fluorescence intensities of their western blots (normalized to their own house-keeping protein Beta Actin) from the small intestinal LM-MP of post-natal mice at different ages during their maturation.

(b) Ratio of HGF:GDNF expression values from the full thickness small intestinal tissue of adult mice taken from the Tabula muris database⁹⁴.

(c) Analyses of human plasma proteome from 1,025 individuals was performed by Sathyan et al.⁹⁵ and the correlation of various plasma proteomes with aging was calculated. Of the genes that showed significant correlation with age (p < 0.05), we found that the HGF levels in human plasma correlated positively with age (above 0), while the levels of GDNF and RET correlated negatively with age (below 0).

Suppl Fig 10: Identification of the MENs cluster in the P21 scRNAseq data and expression of Ret in NENs and in the Ret heterozygous mice.

(a) UMAP representation of all the cell clusters generated from scRNAseq performed on all cells from the P21 murine LM-MP.

(b) Projection of scRNAseq from the P21 murine LM-MP cells into the four MEN-specific NMF patterns generated using scRNAseq on 6-month-old murine LM-MP cells shows higher usage of the four NMF patterns in the putative MENs cluster in the P21 dataset (black arrow).

(c) Expression of the genes Ret and Met in the MENs and NENs clusters in the scRNAseq of P21 LM-MP cells shows that while Met-expression is detected in both MENs and a subset of NENs, the MET-expressing NENs subset shows significant expression of Ret, while the MENs do not express Ret.

(d) Sparklines plot of gene expression of various neuronally significant genes in the scRNAseq data from Zeisel et al., which was generated using flow sorted Wnt1-cre:Rosa26-tdTomato⁺ cells from a P21 mouse, shows that Met expression is observed in the clusters ENT6, 8, and 9 which also show expression of CGRP-expressing gene Calcb, the neural crest-specific gene Ret, Choline acetyltransferase-expressing gene Chat and the pan-neuronal marker PGP9.5- expressing gene Uchl1.

(e) Fold change in transcript abundance of the gene Ret in the small intestinal LM-MP tissue from Ret wildtype (Ret WT) and Ret heterozygous (Ret Het) littermate adult mice shows that Ret transcripts are significantly less abundant in the tissue of Ret Het mice, when compared to those from Ret WT mice. Ret transcript expression in every sample was normalized to the expression of the house keeping gene Hprt. Student’s t-test; * p value < 0.05.

Suppl Fig 11: Effect of GDNF treatment on MHCst-expressing MENs and RET-expressing NENs in aging mice

Representative image of MHCst (red) and Hu (green)-immunostained myenteric ganglia of (a) saline-treated control and (b). GDNF-treated 17-month-old mice show MHCst⁺ (yellow arrows) and MHCst- neurons. Note the reduction in MHCst-immunostaining on inter-ganglionic fibers in the 17-month-old mice when treated with GDNF, compared to control mice. Nuclei are stained with DAPI (blue). Scale bar indicates 10 µm.

Representative image of RET-immunostained (green) myenteric ganglia of (c). saline-treated control and (d). GDNF-treated 17-month-old mice. Nuclei are stained with DAPI (blue). Scale bar indicates 10 µm.

(e) Western blot image of a 4-20% gel performed using the RET antibody on total protein lysate from small intestinal LM-MP layer from post-natal mice shows the expected band between150 and 250 kDa of RET protein.

Suppl Fig 12: Projection of the bulk RNA sequencing data of intestinal tissue from Control and Patients with Obstructed Defecation into our murine scRNAseq-derived NMF patterns using projectR.

Transcriptomic data from intestinal tissues of patients with normal intestinal motility and those with obstructed defecation (OD) was procured from GEO (GSE101968). Non-negative matrix factorization (NMF), as implemented in the R package NNLM (https://github.com/linxihui/NNLM), was performed on the murine scRNAseq data using k=50 and default parameters; and cell weights for each pattern were grouped. Using projectR, the log2 expression (log2(rpkm + 1)) from the human bulk RNA sequencing data from control and OD patients were projected into the murine scRNAseq-derived NMF patterns. The mean projection weights from Control and OD groups were tested for statistically significant differences using Students’ t tests.