key: cord-0255803-anm6eb27 authors: Li, Tao; Morselli, Marco; Su, Trent; Mulugeta, Million; Larauche, Muriel; Pellegrini, Matteo; Taché, Yvette; Yuan, Pu-Qing title: Comparative transcriptomics reveal highly conserved regional programs between porcine and human colonic enteric nervous system date: 2022-02-25 journal: bioRxiv DOI: 10.1101/2022.02.24.480770 sha: 352d600c8ec1ebfd48d1452058b80f1eaeaaa514 doc_id: 255803 cord_uid: anm6eb27 The porcine gut is increasingly regarded as a useful translational model. The enteric nerve system (ENS) in the colon coordinates diverse functions. However, knowledge of the molecular profiling of porcine ENS and its similarity to that of human is limited. We identified the distinct transcriptional programs associated with functional characteristics between inner submucosal and myenteric ganglia (ISG, MG) in porcine proximal and distal colon (p-pC, p-dC) using bulk RNA sequencing (RNA-seq) and single-cell RNA-seq. Comparative transcriptomics of MG in corresponding colonic regions of porcine and human revealed highly conserved programs existing in p-pC and p-dC, which explained >90% of their transcriptomic responses to vagal nerve stimulation (VNS), suggesting that p-pC and p-dC could serve as predictors in translational studies. The conserved programs specific for inflammatory modulation were displayed in porcine with VNS. This study provides a valuable transcriptomic resource for understanding of human colonic functions and neuromodulation using porcine model. The porcine gastrointestinal (GI) tract is increasingly regarded as a useful translational research model due to its structural and functional similarities with human. Compared to mouse, rat, dog, cat or horse, both porcine and human are omnivores and colon fermenters, and have similar metabolic and intestinal physiological processes and microbial composition [1] [2] [3] [4] . Among the important similarities relevant to colonic function, both taenia and haustra/sacculation are present in porcine and human colon, while missing in dogs and most rodents, despite the porcine proximal colon being orientated in a spiral fashion 2, 5 . Additionally, compared with human, porcine has a similar chromosomal structure and a comparably sized genome with a 60% sequence homology and 93% correspondence in relevant biomarkers 6, 7, 8 . The porcine size is also relevant closely mimicking drug dose volumes used in human and for evaluating medical devices 9 . Thus, the porcine model has untapped potential for a greater understanding of the underlying colonic functions of human in health and diseases 10 . The GI tract is highly innervated by both extrinsic (sympathetic and parasympathetic) and intrinsic (enteric nervous system, ENS) nervous systems. The ENS is endowed with complex reflex circuits that control a variety of GI physiological functions via networks of neurons and glia interacting with cells located in the mucosal and muscular layers 11 . Uniquely the ENS is the only part of the peripheral nervous system that can function independently of central nervous system innervation. Therefore, an intact ENS is essential for life and its dysfunction is often linked to digestive disorders 11, 12 . In mammalian colon, the ENS comprises two main ganglionated plexuses: myenteric plexus, located between longitudinal and circular layers of muscularis externa, and submucosal plexuses (SMP) within the submucosa. Unlike the SMP in rodents that is single-layered 13 , both porcine and human colonic SMP consists of an inner submucosal plexus , located close to muscularis mucosae, and an outer submucosal plexus, located on the luminal side of the circular muscle layer 14 . Despite the anatomical similarity, it is still unknown to what extent transcriptional programs are conserved between porcine and human colonic ENS. Such information is of importance for the use of porcine model in translational studies. Recent studies in mice have demonstrated heterogeneity in neuronal identity, morphology, projection orientation and synaptic complexity within the ENS located in different colonic segments in line with the diversity of colonic motility patterns 15 . The region-dependent molecular characterization of porcine ENS has been little investigated so far. In addition, colonic functions are modulated through parasympathetic innervation, responsible for regulating colonic secretion and motility 16 . We have recently reported that electrical vagal nerve stimulation (VNS) triggered pan-colonic contractions in the porcine 17 . However, it is unknown how such VNS influences colonic transcriptional programs and whether the impacts of VNS are region-specific. The incomplete ENS characterization is largely due to longstanding technical challenges in isolating enteric neurons 18 . Laser-capture microdissection (LCM) has been the method of choice for accurately targeting and capturing cells of interest from a heterogeneous tissue sample 19 like the colon 20 . The enriched cells can be used for bulk RNA sequencing (RNA-seq), which becomes a widely adopted method for profiling transcriptomic variations 21 and is a promising choice to characterize the region-dependent molecular profiling in ENS. More recently, single-cell RNA sequencing (scRNA-seq) has emerged as a powerful method for exploring gene expression profile at single-cell resolution 22 . However, its high cost limits its application across a large population of samples 23 . Considering this constraint, several computational methods have been exploited to deconvolve cell type composition in bulk RNA-seq data using scRNA-seq reference datasets [24] [25] [26] [27] which largely facilitate the use of scRNA-seq data from a small number of subjects in clinical setting involving a large number of subjects. In the present study, we profiled the transcriptomes of myenteric and inner submucosal ganglia (MG, ISG) from porcine proximal, transverse and distal colon (p-pC, p-tC, p-dC), and MG in parallel from human ascending, transverse and descending colon (h-aC, h-tC, h-dC) using LCM coupled with bulk RNA-seq analysis. By cross-comparing the transcriptional profiling of MG in the corresponding colonic regions between porcine and human, we aimed to characterize the conservation of gene programs derived from orthologous genes. Then, based on the high-quality orthologous genes, we evaluated discrepancy of the region-specific gene programs in MG and ISG after electrical stimulation of the celiac branch of the abdominal vagus nerve in anesthetized porcine. Additionally, we performed scRNA-seq coupled with bulk gene expression deconvolution to identify cell-type marker genes and putative interactions among neurons and glia in MG of porcine colon and further revealed the regional transcriptomic heterogeneity in the ENS of porcine with and without VNS. The obtained empirical cumulative density function (ECDF) profiles and both Kolmogorov-Smirnov and Pearson's Chi-squared test (p > 0.05) based on the averaging normalized transcription levels of each orthologous gene demonstrated that gene expression profiling in MG of human colonic segments (h-aC, h-tC, h-dC) could be predicted according to that of corresponding porcine colonic segments (p-pC, p-tC, p-dC), following the probability distribution shown in ECDF profiles (Fig. 1a, c) . In addition, pair-wise comparison addressed a significant Spearman's correlation (p < 0.001) (Fig. 1b) . The differentially expressed genes (DEG) enriched significantly (p < 0.05) by Gene Ontology terms were selected for comparisons among colonic segments, between MG and ISG, and the porcine and human. The functional linkages related to the orthologous DEG were extracted to reflect functional similarities across species. Twelve porcine orthologous DEG (9 upregulated and 3 downregulated) in comparison of h-aC-MG and h-dC-MG were involved in the regulatory networks, mediated by 483 DEG (271 upregulated and 212 downregulated), when comparing p-pC-MG and p-dC-MG (Fig. 1d) . The functional linkages related to the 12 porcine orthologous DEG which achieved 100% linkage coverage and represented all the characteristics of p-pC-MG and p-dC-MG. Moreover, these conserved functional linkages accounted for more than 90% of all functional linkages of the regulatory networks responding to VNS when comparing p-pC-MG and p-dC-MG. However, partial coverage was shown when matching the functional linkages derived from human orthologous DEG in comparison of p-pC-MG and p-dC-MG and DEG in comparison 6 of h-aC-MG and h-dC-MG (Fig. 1d) . We did not find any porcine orthologous DEG when comparing h-aC-MG or h-dC-MG with h-tC-MG. Therefore, only p-pC and p-dC were considered in the subsequent analyses. Preferential enrichment for genes associated with MG and ISG functions in porcine proximal or distal colon The functional linkages related to each DEG list in our study and those related to the DEG list screened by the high-quality human orthologous genes achieved almost 100% linkage coverage and the cell-subtype gene markers also fell into the high-quality orthologous gene list. Two databases (Gene Ontology Biological Processes and WikiPathways) were used to interpret the RNA-seq data. The interactions were defined when the networks mediated by DEG involved in WikiPathways of interest and Biological Processes (BPs) shared BPs. We mainly focused on the scenario where there were interactions between the pathways of interest unless otherwise stated. Innate similarities between p-pC and p-dC were found when comparing respective MG and ISG. The upregulated genes in MG of both p-pC and p-dC accounted for a larger percentage of identified genes involved in the top five BPs and led to the higher enrichment of the BPs (Supplementary Fig. 1 ). The BPs connected with those such as nerve development, neuronal projection guidance and neuron/glial cell differentiation. In addition, these genes in MG contributed to higher enrichment scores of the WikiPathways such as oncostatin M signaling pathway, glutamate binding, activation of AMPA receptors and synaptic plasticity, chemokine signaling pathway, nitric oxide biosynthesis, dopaminergic neurogenesis and TGFbeta receptor signaling, while the downregulated genes in MG of both p-pC and p-dC made contribution to higher enrichment scores of the WikiPathways such as activation of angiotensin pathway, acetylcholine synthesis, synaptic vesicle pathway, neurotransmitter release cycle and dopamine transport and secretion ( Fig. 2) . However, there were still some differences found in comparison of MG and ISG between p-pC and p-dC. The upregulated genes in p-pC-MG contributed to higher enrichment scores of the WikiPathways such as nicotine activity on dopaminergic neurons, while the upregulated genes in p-dC-MG exclusively resulted in neurotransmitter uptake and metabolism in glial cells (Fig. 2) . In p-pC-MG, the downregulated genes contributed to larger enrichment scores of the BPs ( Supplementary Fig. 2a, c) , which connected with ENS development, blood vessel development and glial cell differentiation, while the BPs involving the upregulated genes connected with ENS development, synaptic organization and chemical synaptic transmission. In ISG, the enrichment differences between p-pC and p-dC were not as large as those between p-pC-MG and p-dC-MG. The biggest enrichment disparity happened to blood vessel development and the downregulated genes in p-pC-ISG contributed to its higher enrichment score ( Supplementary Fig. 2b, d) . The upregulated genes in p-pC-MG led to larger enrichment scores of the WikiPathways such as acetylcholine synthesis, synaptic vesicle pathway, neurotransmitter receptors and postsynaptic signal transmission, neurotransmitter release cycle and dopamine uptake, unlike those relate to Oncostatin M signaling pathway, cytokines and inflammatory response, chemokine signaling pathway, NO/cGMP/PKG mediated neuroprotection, dopaminergic neurogenesis, TGF-beta receptor signaling and smooth muscle contraction (Fig. 3a) . All listed pathways interacted with smooth muscle contraction. The upregulated genes in p-pC-ISG contributed to most prominent enrichment scores of nicotine activity on dopaminergic neurons, neurotransmitter receptors and postsynaptic signal transmission, dopaminergic neurogenesis, and dopamine biosynthesis process. The downregulated genes in p-pC-ISG resulted in enrichment scores related to nitric oxide biological action, chemokine signaling pathway and TGF-beta receptor signaling (Fig. 3b) . Distinct cell-type-specific pathway enrichment across porcine colonic segments detected by scRNA- The scRNA-seq data showed that there were 10 cell clusters in muscularis externa of p-pC, p-tC (data not shown) and p-dC detected with similar cell markers (Fig. 4a, b) . The neuronal and glial clusters were annotated with markers GAP43 and CLDN11, respectively, which were validated by single-molecule fluorescence hybridization (smFISH). GAP43 and CLDN11 were co-expressed with ELAVL4 and GFAP, respectively ( Fig. 5A -C, M-O). Cell-cell interaction analysis revealed the strongest neuron-glia interaction in p-pC, while the interactions in p-tC (data not shown) and p-dC were almost the same (Fig. 4a, b) . Three neuronal and two glial sub-clusters were further annotated with markers ACLY for cholinergic neurons, Table 1 ). The WikiPathways with top ranking enrichment scores in comparison of p-pC-MG and p-dC-MG were selected to evaluate contribution of each cell subset to pathway enrichment. Among the pathways involving the upregulated genes in p-pC-MG, cholinergic neurons played a leading role and contributed to enrichment of neurotransmitter release cycle, which was also attributed to glutamatergic neurons, enrichment of synaptic vesicle pathway, which was also attributed to SLC41A1_glia, and enrichment of chemokine signaling pathway and TGF-beta receptor signaling (Fig. 6 ). Among the pathways involving the downregulated genes in p-pC-MG, nitrergic neurons contributed to enrichment of dopaminergic neurogenesis, TGF-beta receptor signaling and chemokine signaling pathway. However, glutamatergic neurons contributed more to enrichment of chemokine signaling pathway. Although cholinergic neurons contributed to enrichment of neurotransmitter release cycle, their contribution was much less than that to enrichment of the pathway involving the upregulated genes in p-pC-MG (Fig. 6 ). The enrichment ratio of the same pathway involving upregulated and downregulated genes in each DEG list (e.g., p-pC-MG versus p-dC-MG, or p-pC-MG versus p-pC-ISG) was calculated and compared between porcine with and without VNS to evaluate its influences at both single and bulk resolution. The DEG, derived from comparison of p-dC-MG and p-dC-ISG in porcine with or without VNS, were involved in almost the same WikiPathways, despite different numbers of DEG. Only one differentially regulated gene under VNS was found, but enrichment of the WikiPathway mediated by the gene was not significant (Supplementary Fig. 4a) , suggesting a significant overlap in the transcriptional landscape under these two situations. Similarly, comparison of p-pC-ISG and p-dC-ISG revealed that no differences were observed in porcine with and without VNS ( Supplementary Fig. 4b ). Compared with porcine without VNS, BP enrichment patterns in porcine with VNS were almost the same based on comparison of p-pC-MG and p-pC-ISG. There were higher enrichment scores of the BPs involving upregulated genes in p-pC-MG than those involving downregulated genes in p-pC-MG. However, there were some alterations of BP enrichment patterns in porcine with VNS based on comparison of p-pC-MG and p-dC-MG. Enrichment scores of the BPs such as G protein-coupled receptor signaling pathway, anatomical structure morphogenesis and system process, involving upregulated genes in p-pC-MG, were higher than those involving downregulated genes in p-pC-MG (Supplementary Fig. 1 (Fig. 9, 10 ). The cross-comparison of the regional transcriptomic profiling between porcine and human indicates that the myenteric ganglia in proximal or distal colon express a conserved core gene program of orthologous genes from porcine to human with regional heterogeneity. Interestingly, VNS on the porcine modulates the conserved transcriptional program in the myenteric ganglia of proximal or distal colon, which may be implicated in the VNS-altered colonic physiological functions as we reported before 17 , suggesting that human colonic functions and neuromodulation could be modeled in the porcine. Recent studies have identified several conserved marker genes in ENS cell subtypes between human and mouse based on orthologous gene expression 18, 28 . According to our bulk RNA-seq data from both porcine and human ENS, we uncovered 7246 and 2601 conserved genes between h-aC-MG and p-pC-MG (p-value > 0.001) and between of h-dC-MG and p-dC-MG (p-value > 0.001) respectively from a total of 12291 highquality orthologous genes using R package SCBN. The functional linkages of gene regulatory networks, where single genes are positioned, have been used as powerful indicator, reflecting mechanistic insights and functional similarities across species [29] [30] [31] [32] . Pathway analysis could identify biological pathways that are enriched in DEG lists and highlight the overlaps among pathways based on the shared genes 29 . We extracted the gene networks to indicate the differentiation between h-aC-MG and h-dC-MG and between p-pC-MG and p-dC-MG. Twelve porcine orthologous genes from comparison of h-aC-MG and h-dC-MG were found, of which only 2 genes were the conserved genes based on the outputs from R package SCBN. Nonetheless, the 12 orthologous genes served as driver genes, whose networks represented all the characteristics of p-pC-MG and p-dC-MG in porcine, in terms of 100% linkage coverage (Fig. 1) . It was worth noting that the networks derived from the 12 porcine orthologous driver genes interpreted >90% functional alterations in p-pC-MG and p-dC-MG of porcine subjected to VNS (Fig. 1 ). The networks mediated by the human orthologous genes from comparison of p-pC-MG and p-dC-MG did not realize 100% coverage of the linkages derived from the DEG lists in comparison of h-aC-MG and h-dC-MG, probably because of the lack of complete annotation and assembly of the porcine genome 33 . Thus, although there were few conserved genes between porcine and human colon, the networks of interacting orthologous genes orchestrated almost the same complex colonic functions across species, and p-pC and p-dC could serve as predictors, with important implications in understanding of human colonic functions and neuromodulation. Our work shows for the first time the ENS transcriptomic heterogeneity between MG and ISG and along the proximal-distal axis of the porcine colon. Comparisons of MG and ISG in each of p-pC and p-dC uncovered several pathways responsible for normal colonic physiological functions such as peristalsis and secretion. The inflammatory signaling pathways prevailed in MG, including interactions between immune cells and neurons (e.g., oncostatin M signaling pathway and cytokines and inflammatory response), cytokine production (e.g., chemokine signaling pathway) and anti-inflammation (e.g., TGF-beta receptor signaling) (Fig. 2) . Of them, the enrichment of anti-inflammation, involving >99% of the DEG, was the highest. The reason for the inflammatory signaling pathways predominant in MG is probably because the inflammatory bowel disease (IBD) risk genes (e.g., SMAD3 and OSMR) were significantly expressed in MG. Our data also showed that all inflammatory signaling pathways closely interacted with smooth muscle contraction, indicating that the inflammatory status influenced smooth muscle activity. In addition, effects of nitric oxide mediated by nNOS and neurotransmitter clearance mediated by SLC6A4 were dominant in MG. Both nNOS and SLC6A4 that lowers the levels of serotonin are responsible for relaxation of smooth muscle [34] [35] [36] . Furthermore, the enrichment of glutamate binding, the activation of AMPA receptors and synaptic plasticity was much higher in MG than in ISG and directly interacted with the BPs such as the regulation of blood circulation and smooth muscle contraction in MG. By contrast, ISG had higher enrichment of activation of angiotensin pathway, synaptic vesicle pathway, acetylcholine synthesis, binding and downstream events and dopamine transport and secretion (Fig. 2) , all of which interacted with each other. Angiotensinogen (AGT) that was located in the core of activation of angiotensin pathway and selected in our study has been involved in the development of GI mucosal injury due to SARS-CoV-2 and other factors 37, 38 . Therefore, we used the enrichment of this pathway to evaluate the mucosal integrity which showed that only ISG contributed to the modulation of mucosal integrity. Moreover, we found that the activation of angiotensin pathway, acetylcholine synthesis, binding and downstream events and dopamine transport and secretion shared the BPs namely signaling and cell communication, suggesting the potential roles of acetylcholine-dopamine balance in maintenance of mucosal integrity. It is worth noting that neurotransmitter uptake and metabolism in glial cells, mediated by SLC1A3, a glial high affinity glutamate transporter, was dominant in MG only in p-dC (Fig. 2) . A steep oxygen gradient along the proximal-distal axis of colon exists, which leads to formation of reactive oxygen species and the resulting glutamate excitotoxicity [39] [40] [41] . Glutamate uptake in glia guarantees the neural circuit integrity because SLC1A3 connected the BPs such as synaptic signaling and neurogenesis. When comparing MG between p-pC and p-dC, we found that the inflammatory signaling pathways such as cytokines and inflammatory response, chemokine signaling pathway and TGF-beta receptor signaling were prevalent in p-dC-MG (Fig. 3) . There is evidence that hypoxia and inflammation are intertwined 41 . In addition, p-dC-MG compared to p-pC had much higher enrichment of smooth muscle contraction, NO/cGMP/PKG mediated neuroprotection and dopaminergic neurogenesis (Fig. 3) , with all of them having interactions. It has been reported that GUCY1A1 (guanylate cyclases), which was involved in NO/cGMP/PKG mediated neuroprotection, is activated by nitric oxide and plays a key role in inhibiting neuroinflammation 42 . From this perspective, we infer that the combination of nitrergic and dopaminergic neurons contributes to the regulation of smooth muscle contraction. This was confirmed by our scRNA-seq data showing that the anti-inflammatory status, determined by more BPs shared by TGF-beta signaling, chemokine signaling pathway and anti-inflammatory signaling, was closely related to the differentiation of two neuronal populations (Fig. 6) . This is consistent with recent reports, implicating the role of TGF-beta in differentiation of nitrergic and catecholaminergic enteric neurons 43, 44 . Recent study showed that the interactions of these two neuronal populations could stabilize asynchronous contractile activity and lead to the generation of colonic peristalsis in proximal colon 45 . However, our results suggested that the contribution of dopaminergic and nitrergic neurons to colonic peristalsis is much greater in p-dC than in p-pC. By contrast, p-pC-MG had more complex neuronal wiring, highlighted by the high enrichment of nicotine activity on dopaminergic neurons (mediated by TH), acetylcholine synthesis (mediated by ChAT), synaptic vesicle pathway, neurotransmitter receptors and postsynaptic signal transmission (mediated by HTR3A), neurotransmitter release cycle and GABA receptor signaling (mediated by GABRB3 and GABRA4) (Fig. 3) . Cholinergic neurons played leading roles based on their contribution to the top five WikiPathways (Fig. 6) , which is in accordance with the findings of Li et al 15 . The comparison of ISG between p-pC and p-dC revealed that the inflammatory signaling pathways such as chemokine signaling pathway and TGF-beta receptor signaling prevailed in p-dC-ISG, while in p-pC-ISG, there was higher synaptic activity, characterized by the high enrichment of nicotine activity on dopaminergic neurons (mediated by TH) and neurotransmitter receptors and postsynaptic signal transmission (mediated by HTR3A and HTR3B) (Fig. 3) . However, we also found high enrichment of blood vessel development and effects of nitric oxide (mediated by nNOS) in p-dC-ISG. Nerve fibers may contribute to the perivascular circulation based on the distribution of nNOS and its catalytic activity with nNOS-derived NO maintaining vasodilation 46 . Unlike the comparison of MG between p-pC and p-dC, dopaminergic neurogenesis and dopamine biosynthetic process were predominant in p-pC-ISG. Evidence showed that dopamine could promote colonic mucus secretion 47 . Our previous study has indicated that the functional sphere of the influence of VNS on porcine colon was pan colonic, from p-pC to p-dC, and VNS increased contractions across the colon 17 . The evidence also showed that the vagal nerve endings often synapse onto neurons in the myenteric plexus 48 . Consistent with this report, our data showed that the major changes induced by VNS happened in MG, especially in p-pC-MG. Although VNS led to some differences of pathway enrichment patterns when comparing p-pC-MG with p-pC-ISG, the alterations were not as prominent as those in comparison to p-pC-MG and p-dC-MG characterized by the enrichment pattern reversal of more than half the investigated pathways (Fig. 7) . Compared with p-pC-MG, the obvious changes in p-dC-MG were the enhanced enrichment of smooth muscle contraction and glutamate binding, activation of AMPA receptors and synaptic plasticity, which were interactive. This suggested the crucial roles of synaptic plasticity mediated by glutamate in regulation of smooth muscle contraction. Interestingly, cholinergic neurons contributed to the enrichment of glutamate binding, activation of AMPA receptors and synaptic plasticity, implying that cholinergic neurons may corelease both acetylcholine and glutamate, and the released glutamate could activate postsynaptic neurons (glutamatergic neurons) 49 . Our cell-cell interaction data also confirmed the potentially upregulated connection between cholinergic and glutamatergic neurons under VNS (Fig. 10b) . Together with the improved enrichment of GABA receptor signaling, we speculated that VNS promoted the occurrence of rhythmic phasic contractions and fecal pellet expulsion in p-dC. However, in p-pC-ISG, there was the decreased enrichment of acetylcholine binding and downstream events, neurotransmitter release cycle and synaptic vesicle pathway. Of them, VNS led to the highest alterations in the enrichment of glutamate release cycle (GLS/GLS2). This may be related to the VNS-induced the activation of an anti-inflammatory response 40, 41 . This is in accord with the improvement of the mucosal integrity, characterized by the reduced enrichment of activation of angiotensin pathway (Fig. 7b) . Gap junctions greatly contributes to electrical transmission 50 . Here, we used the enrichment changes of gap junctions to interpret the influences of VNS on the transcriptomic profiling. When comparing p-pC-MG with p-pC-ISG, the enrichment of gap junctions, anti-inflammatory signaling and smooth muscle contraction was surprisingly reduced in p-pC-MG by VNS (Fig. 7b) and had a significant positive correlation ( Supplementary Fig. 5 ). In addition, we found that glutamate release cycle and gap junctions showed the highest similarity score, suggesting that glutamate transmission may complement the function of gap junctions. Consequently, glutamate release triggered acetylcholine binding and downstream events (CHRNA5) due to their significantly positive correlation ( Supplementary Fig. 5 ). The potentially synergistic functions of glutamate and acetylcholine may be important for functional modulation in p-pC-MG. Collectively, the alterations of the pathway enrichment under VNS may reduce disease risks in view of the decreased enrichment of human orthologous risk genes for intestinal and extra-intestinal diseases, associated with gut dysmotility and inflammation in p-pC-MG (Fig. 7b) . Compared with p-dC-MG, enrichment of gap junctions was highly increased in p-pC-MG, while enrichment of synaptic vesicle pathway, neurotransmitter receptors and postsynaptic signal transmission and neurotransmitter release cycle was decreased. We speculated that there is direct mechanistic relationship between the formations of neuronal gap junction coupling and the disappearance of chemical transmission in most neuronal cell types as reported 50 . Previous studies showed that VNS reduces intestinal inflammation through the cholinergic anti-inflammatory pathway 51 . In line with the statement, we found that the enrichment ratio of pro-to anti-inflammatory signaling (pro/anti) in porcine without VNS was 1.5 times higher than that in porcine with VNS, suggesting the anti-inflammatory effects of VNS (Fig. 7a) . The enrichment of pro-inflammatory signaling under VNS was attributed to some extent to nitrergic neurons as supported by a recent report showing that NOS1 led to an outburst of inflammatory reactions in macrophages via activator protein-1 52 . Whereas, enrichment of anti-inflammatory signaling was attributed to both cholinergic and glutamatergic neurons and SLC41A1_Glia in p-pC-MG (Fig. 8, 9 ). The enteric glia were thought to be essential to reduce the intestinal inflammatory response under VNS 53 . In the vagal circuitry, glutamate is a primary neurotransmitter involved in key gastrointestinal functions 54 . However, its anti-inflammatory effects are rarely reported except that excessive release of glutamate into the extrasynaptic space led to inflammatory responses, which could be corrected by glial cells in CNS 55 . Our result indicated that the synergies of glutamatergic neurons with enteric glia were also potentially present in ENS under VNS because of the increased enrichment of neurotransmitter uptake and metabolism in glial cells (GLUL, conversion of glutamate and ammonia to glutamine) p-pC-MG (Fig. 7a) . Duo to the direct interactions between smooth muscle contraction and anti-inflammatory signaling, the anti-inflammatory properties indeed reduced the enrichment of human orthologous risk genes for intestinal and extra-intestinal diseases, associated with gut dysmotility and inflammation (Fig. 7a) . We found that OSMR, the only risk gene for IBD was expressed significantly in p-pC-MG compared to p-dC-MG under VNS. OSMR was enriched in SLC41A1-Glia and the OSMR-mediated pathways shared BPs with pro-and anti-inflammatory signaling. The enrichment ratio of the shared BPs (pro/anti) was 1.3 in porcine with VNS, while pro/anti was 1.2 in porcine without VNS, which did not show the significant anti-inflammatory effects of VNS on IBD. Future work should involve more risk genes for IBD to convincingly identify the exact correlation between vagal function and IBD. Finally, we derived the VNS influences on p-dC-MG from the improved enrichment of gap junctions in SLC41A1-Glia of p-dC-MG, compared with p-pC-MG (Figs. 8d, 9) . Surprisingly, unlike the pro-inflammatory impacts of nitrergic neurons in p-pC-MG, they played the antiinflammatory roles in p-dC-MG, together with glutamatergic neurons and SLC41A1_Glia. Moreover, the contribution of SLC41A1_Glia to anti-inflammation was much higher in p-dC-MG. Interestingly, the effects of VNS on glutamatergic neurons and SLC41A1_Glia in p-dC-MG led to switch of their contribution from inflammation to anti-inflammation (Fig. 9) . Different from the attribution of glutamatergic neurons to anti-inflammation in p-pC-MG under VNS, the increased enrichment of glutamate binding, activation of AMPA receptors and synaptic plasticity in cholinergic neurons ensured the correction of glutamatergic neurotransmission dysfunction in p-dC-MG (Fig. 7a) . Thus, we inferred that the synergies of glutamatergic neurons with cholinergic neurons contributed to anti-inflammation in p-dC-MG. In summary, this study compared for the first time the transcriptomic profiles in the colonic ENS between porcine and human and across two layers of enteric plexus along the proximal to distal colonic regions in porcine. The regional-specific gene programs were identified at both bulk and single-cell resolution in porcine and regional-dependent highly conserved core transcriptional programs were revealed between porcine and human. VNS exerted its effects mainly on the myenteric ganglia in porcine. Regional-specific transcriptomic responses to VNS shared >90% of the conserved core transcriptional programs between porcine and human. Our study provides a fundamental resource for better understanding of the importance of porcine model in the colonic translational research. HiSeq 3000 sequencer as 50 base pair single-end reads. All of the reads were then aligned to the sus_scrofa or human genome using STAR and the DEG lists were generated using edgeR (Supplementary method 3) and used for pathway enrichment analysis. The mucosa and submucosa were peeled off from each p-pC, p-tC and p-dC of 4 naïve porcine (2 of each sex). The muscularis externa containing myenteric ganglia were processed for cell suspension preparation The scRNA-seq data were processed for cell selection, filtration, normalization, scaling, cell cluster identification and annotation of neuronal and glial clusters using R toolkit Seurat (Supplementary method 5). The neuronal and glial populations acquired from 12 scRNA-seq datasets were used to further gain the subpopulations of neurons and glia using Seurat. A total of 511, 286 and 443 neurons and 316, 318 and 389 glia were isolated from p-pC, p-tC and p-dC respectively, from which 3 neuronal (cholinergic, nitrergic and glutamatergic) and 2 glial subpopulations (marked by SLC41A1 and SPC24, respectively) were identified. Cell-cell interactions were inferred from the single-cell transcriptomic data based on the ligand-receptor (L-R)-pair list in porcine generated using R package 'LRBase.Ssc.eg.db'. Firstly, scRNA-seq data were read using R package 'DropletUtils' and then the target cells were selected using unique molecular identifiers of neurons and glia. After data normalization, the R package 'scTGIF' was used to convert Ensembl ID to NCBI gene ID and a SingleCellExperiment object was created using R package 'SingleCellExperiment'. At last, the R package 'scTensor' was used to detect and visualize cell-cell interactions. To validate the key discriminatory markers for the subsets of neuronal and glial cells clustered from scRNAseq data, we performed RNAscope, a smFISH technique with fresh-frozen sections ( Images were taken using a Zeiss LSM 710 confocal microscope (Carl Zeiss Microscopy, LLC, White Plains, NY) with a laser set of 488 (C1)/561 (C2)/647 (C3). Five to ten optical sections (Z-stack) crossing the myenteric ganglia with frame 436×436 μm and 1 μm apart (20× objective) were acquired and overlaid using Imaris 9.7 for Neuroscientists (Bitplane Inc., Concord, MA). Deconvolution of bulk RNA-seq data using scRNA-seq datasets from colonic enteric ganglia of naïve porcine The R package 'SCDC' was used to deconvolve the bulk RNA-seq data from the myenteric ganglia of naïve porcine colon with an ENSEMBLE method using default settings 24 . The 12 scRNA-seq datasets generated in this study were served as references. With marker genes specific for the specified cluster, five subpopulations of cells including cholinergic, nitrergic, glutamatergic neurons and two types of glia were identified in scRNA-seq datasets and selected to construct basis matrix in p-pC, p-tC and p-dC. For each scRNA-seq dataset with raw counts, a quality control procedure was carried out to filter the cells with the threshold of keeping a single cell from an assigned cluster of 0.7. The grid search method was used to derive the ENSEMBLE weights with the search step size of 0.01. SpearmanY from the grid search result represented the maximum of the Spearman correlation between observed cell-type gene expression in scRNA-seq datasets and predicted cell-type gene expression (probabilities) in bulk RNA-seq data. The celltype gene lists were then matched to the DEG lists (p < 0.05) that were obtained from bulk RNA-seq data analysis. In this way, the cell-type DEG lists were generated and used for pathway enrichment analysis. A gene ortholog table was first established using human genome as the reference gene list to compare transcription between human and porcine. Gene homology search was performed, using ensemble multiple species comparison tool (http://www.ensembl.org/biomart/). A high-quality orthologous gene list was extracted with whole genome alignment score above 75 56 , resulting in 12291 genes. Based on the reads of the high-quality ortholog genes, the Python/bioinfokit (v0.9.1) package was used to calculate standard RPKM (reads per kilobase of exon model per million mapped reads) expression values for the orthologous gene set, followed by log2 transformation. The cross-species normalization was completed according to the scaling procedure reported by Brawand et al. 57 Correlation, which was visualized using R/corrplot package. The DEG lists were generated based on the transcriptomic comparisons in both human and naïve porcine samples with datasets pooled from males and females. Those from the male porcine were only used to assess the effects of VNS on the transcriptomic profiling in porcine colon. In order to remove the potential contamination from intestinal muscle during LCM, unique smooth muscle marker genes derived from scRNA-seq were deleted from each DEG list based on bulk RNA-seq data. In the same way, gene network overlap was assessed based on the human orthologous DEG list from the regional comparisons in porcine and the regionally corresponding DEG list from human or the regionally corresponding DEG lists from porcine with and without VNS. The human disease risk genes relevant to Hirschsprung's disease, inflammatory bowel disease, autism spectrum disorders, and Parkinson's disease were selected based on the literatures 58-61 . Their orthologous genes were defined in the porcine colonic ENS and applied for genetic risk prediction of human diseases. Enrichment analysis of WikiPathways that regulate the expression of the disease risk genes was performed based on the above-mentioned procedure. RNA sequencing data from this study have been deposited into the National Center for Biotechnology Information Gene Expression Omnibus under accession number: GSE197106. Arrowheads in light blue and purple indicate the enrichment ratio with and without VNS, respectively. Pro, pro-inflammation. Anti, anti-inflammation. The pig as an experimental model for elucidating the mechanisms governing dietary influence on mineral absorption Comparison of the gastrointestinal anatomy, physiology, and biochemistry of humans and commonly used laboratory animals The pig as a model for human nutrition Inter-species transplantation of gut microbiota from human to pigs Comparative physiology of the mammalian colon and suggestions for animal models of human disorders Comparative analyses of multi-species sequences from targeted genomic regions Dynamics of mammalian chromosome evolution inferred from multispecies comparative maps Large animal models: The key to translational discovery in digestive disease research Importance of the pig as a human biomedical model Porcine models of digestive disease: the future of large animal translational research The bowel and beyond: the enteric nervous system in neurological disorders The enteric nervous system and neurogastroenterology Classes of enteric nerve cells in the guinea-pig small intestine Structural organization and neuropeptide distribution in the mammalian enteric nervous system, with special attention to those components involved in mucosal reflexes Regional complexity in enteric neuron wiring reflects diversity of motility patterns in the mouse large intestine Central nervous system control of gastrointestinal motility and secretion and modulation of gastrointestinal functions The effect of colonic tissue electrical stimulation and celiac branch of the abdominal vagus nerve neuromodulation on colonic motility in anesthetized pigs The human and mouse enteric nervous system at single-cell resolution Laser capture microdissection in the tissue biorepository An optimized protocol harnessing laser capture microdissection for transcriptomic analysis on matched primary and metastatic colorectal tumours limma powers differential expression analyses for rna-sequencing and microarray studies Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets Computational and analytical challenges in singlecell transcriptomics SCDC: bulk gene expression deconvolution by multiple single-cell RNA sequencing references Bulk tissue cell type deconvolution with multi-subject single-cell expression reference Determining cell type abundance and expression from bulk tissues with digital cytometry Accurate estimation of cell-type composition from gene expression data Combinatorial transcriptional profiling of mouse and human enteric neurons identifies shared and disparate subtypes in situ Pathway enrichment analysis and visualization of omics data using g:Profiler, GSEA, Cytoscape and EnrichmentMap Functional alignment of regulatory networks: a study of temperate phages Gene regulatory network inference: evaluation and application to ovarian cancer allows the prioritization of drug targets A comparative metabologenomic approach reveals mechanistic insights into Streptomyces antibiotic crypticity Advances in porcine genomics and proteomics-a toolbox for developing the pig as a model organism for molecular biomedical research Nitric oxide in gastrointestinal health and disease Serotonin transporter gene (SLC6A4) polymorphism in patients with irritable bowel syndrome and healthy controls Meta-analysis of smooth muscle relaxants in the treatment of irritable bowel syndrome Potential neuroinvasive pathways of SARS-CoV-2: Deciphering the spectrum of neurological deficit seen in coronavirus disease-2019 (COVID-19) Renin-Angiotensin system associated with risk of upper GI mucosal injury induced by low dose aspirin Reactive oxygen species, nutrition, hypoxia and diseases: Problems solved? A circuit-dependent ROS feedback loop mediates glutamate excitotoxicity to sculpt the Drosophila motor system Hypoxia and Inflammation The role of NO/cGMP signaling on neuroinflammation: A New Therapeutic Opportunity BMP2 promotes differentiation of nitrergic and catecholaminergic enteric neurons through a Smad1-dependent pathway Bone Morphogenetic Protein (BMP) signaling in development and human diseases Role of enteric dopaminergic neurons in regulating peristalsis of rat proximal colon Nitric oxide in the vasculature: where does it come from and where does it go? A quantitative perspective Dopamine promotes colonic mucus secretion through dopamine D5 receptor in rats The enteric network: Interactions between the immune and nervous systems of the gut Habenula "cholinergic" neurons co-release glutamate and acetylcholine and activate postsynaptic neurons via distinct transmission modes Electrical synapses and their functional interactions with chemical synapses The vagal innervation of the gut and immune homeostasis NOS1 mediates AP1 nuclear translocation and inflammatory response Enteric glia cells are critical to limiting the intestinal inflammatory response after injury Receptors and transmission in the brain-gut axis. II. Excitatory amino acid receptors in the brain-gut axis Inflammation, glutamate, and glia: A trio of trouble in mood disorders Cross-species single-cell analysis reveals divergence of the primate microglia program The evolution of gene expression levels in mammalian organs Mouse models of Hirschsprung disease and other developmental disorders of the enteric nervous system: Old and new players Intra-and inter-cellular rewiring of the human colon during Ulcerative Colitis Autism Sequencing Consortium. Insights into Autism Spectrum Disorder Genomic Architecture and Biology from 71 Risk Loci International Parkinson's Disease Genomics Consortium; 23andMe Research Team. A meta-analysis of genome-wide association studies identifies 17 new Parkinson's disease risk loci Comparison of cell-type specific pathway enrichment in myenteric ganglia (MG) between porcine proximal and distal colon (p-pC, p-dC). The bubble plot shows the gene percentage and enrichment of WikiPathways The The authors declare that they have no conflict of interest.