Automated computational analysis reveals structural changes in the enteric nervous system of nNOS deficient mice

Neuronal nitric oxide synthase (nNOS) neurons play a fundamental role in inhibitory neurotransmission, within the enteric nervous system (ENS), and in the establishment of gut motility patterns. Clinically, loss or disruption of nNOS neurons has been shown in a range of enteric neuropathies. However, the effects of nNOS loss on the composition and structure of the ENS remain poorly understood. The aim of this study was to assess the structural and transcriptional consequences of loss of nNOS neurons within the murine ENS. Expression analysis demonstrated compensatory transcriptional upregulation of pan neuronal and inhibitory neuronal subtype targets within the Nos1−/− colon, compared to control C57BL/6J mice. Conventional confocal imaging; combined with novel machine learning approaches, and automated computational analysis, revealed increased interconnectivity within the Nos1−/− ENS, compared to age-matched control mice, with increases in network density, neural projections and neuronal branching. These findings provide the first direct evidence of structural and molecular remodelling of the ENS, upon loss of nNOS signalling. Further, we demonstrate the utility of machine learning approaches, and automated computational image analysis, in revealing previously undetected; yet potentially clinically relevant, changes in ENS structure which could provide improved understanding of pathological mechanisms across a host of enteric neuropathies.


Machine learning analysis reveals occult changes in the Nos1 −/− ENS.
In order to assess if changes in ENS composition, in terms of neuronal subtype gene expression, leads to changes in ENS structure which are visually indistinguishable, we used readily available machine learning tools; combined with in-house scripting (Supplementary Information), to quantitatively compare ENS network characteristics (e.g., network density, interganglionic area, network branching, junctioning, directionality and coherence) in multiple data sets (i.e., 105 optical sections/distal colon) in both control C57BL/6J and Nos1 −/− distal colon. To ensure reproducibility and capture the 3D-structure of the distal colonic ENS, we devised an imaging protocol whereby 21 optical sections (1 μm) were obtained in a cross-shaped array (Fig. 3a) using confocal microscopy. Following acquisition, imaged stacks were analysed after pre-processing and segmentation (Fig. 3b-e). Interestingly, while visually comparable, machine learning digital analysis revealed that TuJ1 + network density is significantly greater in the Nos1 −/− distal colon (39.2 ± 1.1%; 30 z-projected hpf, n = 6) compared to C57BL/6J controls (29.3 ± 1.9%; 30 z-projected hpf, n = 6; P < 0.0001; Fig. 4a-c). To validate our machine learning segmentation and network coverage analysis, we analysed equivalent images taken from the colon (ganglionic, transitional zone and aganglionic regions) of an Ednrb −/− (Endrb tm1Ywa /J) mouse model which displays variable but graded loss of the ENS along the length of the colon. Here, as expected, our digital analysis revealed that network density was reduced in a graded fashion when comparing ganglionic (18.6%), transition zone (16.7%) and aganglionic (10.1%) regions confirming the validity of our approach (Supplementary Figure 1).

Loss of Nos1 expression causes changes in ENS branching and junctions.
To determine how the observed alterations in overall network coverage reflect specific changes in network composition, we further quantified branch number, branch length and the number of junctions present within the ENS.
Notably, TuJ1 + networks in the Nos1 −/− distal colon demonstrated a higher degree of branching when compared to C57BL/6J control tissues (    Fig. 5g) than that of control C57BL/6J mice. However, despite differences in junction number, total branches and total branch length, mean branch length was not found to be statistically different between the C57BL/6J (24.8 ± 1.2 pxls; 30 z-projected hpf, n = 6) and Nos1 −/− colon (21.5 ± 1.4 pxls; 30 z-projected hpf, n = 6; p = 0.08; Fig. 5h). Together, these data suggest that global loss of nNOS results in the development of a ENS with greater interconnectivity and increased neuronal projections. Similarly, we again sought to confirm the validity of these branching analyses in the Ednrb −/− colon: where progressive loss of myenteric ganglia is known to occur, from ganglionic to aganglionic regions, alongside overgrowth of extrinsic neuronal fibres. Importantly, in this well characterised model, our digital analysis was able to positively identify and quantify this graded loss of the ENS. Here   www.nature.com/scientificreports/ Loss of Nos1 has no effect on network directionality or coherency. Having demonstrated that network coverage, branching and junctioning are altered in the Nos1 −/− distal colon we next sought to assess if these changes led to alterations in network directionality (Fig. 6a-d) and coherency (i.e., the degree to which local features are orientated) 21 . Interestingly, neither network orientation (C57BL/6J: 19.6 ± 12.1° vs Nos1 −/− : − 13.9 ± 11.2°; P = 0.05, 30 z-projected hpf, n = 6) nor coherency (C57BL/6J: 0.10 ± 0.01% vs Nos1 −/− : − 0.09 ± 0.01%; P = 0.40, 30 z-projected hpf, n = 6;) were found to be different (Fig. 6e,f). Although, it appeared that network directionality was different in C57BL/6J vs Nos1 −/− mice, this was due to an analytical artefact in which 'up and down directionality' were registered as ' ± 90°' , respectively, in relation to a fixed point. Taken together these data suggest that loss of Nos1 has no discernible effect on the directionality of the ENS which retains 'random' coverage (i.e., < 1% coherency).

Discussion
Understanding of the overall effects of neuronal loss in the ENS remains rudimentary. Often, in diagnosing gut motility dysfunction, damage to the ENS; or loss of neuronal subtypes, is attributed as the causative factor based on histopathological screening for a panel of common neural, glial and interstitial cell markers. While current methodology may reveal the apparent cause of disease, this approach offers little in network characterisation in terms of structure; barring detectable evidence of aganglionosis, hyper/hypoganglionosis or the presence of ectopic ganglia, and often ignores subtle compensatory factors with could be phenotypically relevant. Similarly, subtle compensatory mechanisms and/or occult alterations in network structure are often not examined in animal models of human disease. In both human and animal studies this is in part due to the sheer scale and complexity of the ENS. However, in not assessing if such changes occur, this reductionist approach may limit our understanding of developmental, and disease, processes within the gut. More recently, efforts have been made to bridge this gap using network quantification methods to highlight differences in neuronal density 22,23 , neuronal fibre numbers and branching using TuJ1 24-26 , or using artificial intelligence (AI)-driven tools in terms of unbiased quantification of neuronal numbers 27 . Here, we extend this AI-driven approach using machine learning methodologies to highlight that occult changes in network structure occur within the Nos1 −/− mouse colon. Notably, we show that, despite gross visual similarities; the Nos1 −/− ENS displayed morphological differences with increases in neuronal density, projections and branching when compared to age matched control mice. We also demonstrate that Vasoactive intestinal peptide (Vip) is upregulated, at the transcriptional level, upon global loss of Nos1; pointing to compensatory mechanisms within ENS development. Previous studies have shown that VIP displays significant colocalization within nNOS + neurons in the mature ENS [28][29][30][31] . Moreover, modulation of nNOS and VIP expression in the ENS (i.e., counterbalance of increasing VIP expression upon reduction of nNOS expression) has been noted in disease states such as Parkinson's disease 32 . Similarly, increases in VIP expression within nNOS + neurons has been shown in response to both pharmacological treatment and culture 33,34 . Our findings add further weight to the hypothesis that such compensation may be a protective mechanism whereby the ENS attempts to rebalance excitatory and inhibitory pathways upon the loss of major nitrergic relaxatory signalling. However, our findings of changes in network structure also suggest that this compensatory mechanism results in significant morphological remodelling within the ENS. Taken together, our findings highlight that, even within a relatively simple transgenic animal model, major changes in ENS network structure and composition occur which can be overlooked. Extrapolating these findings to human disease states, we propose that current methodologies are insufficient to determine the true scope of phenotypically relevant changes which occur upon disruption of the ENS. We acknowledge that a major limitation of our study is the use the Nos1 tm1Plh /J (Nos1 −/− ) knockout model for such characterisation studies. While attractive as a model, we recognise that global loss of nNOS signalling is unlikely to truly reflect disease states in the ENS; wherein loss or disruption of specific, or multiple, neuronal subtypes is likely to occur in a regionalised or sporadic manner. Nevertheless, we demonstrate that occult changes in network structure occur upon global loss of Nos1, however, the extent to which these structural changes influence physiological output is unclear. Furthermore, it is unclear if similar alterations would be observed in human disease states where additional influences such as inflammation are likely to have profound effects on ENS structure and dynamics [35][36][37][38] . Previous human studies have shown dual loss of nNOS neurons and interstitial cells of Cajal (ICC) in multiple disease states [39][40][41][42] . In keeping with these findings, our previous work has demonstrated reduction of ICC networks within the Nos1 −/− global knockout 20 ; suggesting that loss of nNOS signalling has significant effects including remodelling of the neuromusculature. Our current findings add to these data, though raise questions as to the "cause vs consequence" of many of the remodelling effects observed in diseased tissues and animal models. Importantly, global knockout models fail to address the impact of conditional loss of nNOS, in the ENS, which is characteristic of many enteric neuropathies. Therefore, it will be important to investigate how conditional loss of neuronal subtypes impacts upon network structure and the surrounding microenvironment in this context. While the imaging protocol utilised within the current study attempted to limit sampling variation and improve consistency, we appreciate that analysis of the ENS at a larger scale would be more advantageous. Unfortunately, current technical limitations such as low spatial resolution in 'z' , at low magnification, and vignetting in stitched high-power images, limited our ability to apply similar analytics to the ENS on a whole organ scale. Indeed, our failure to detect potential differences in network directionality is likely due to limitations in our imaging protocol, which limited the scope of any particular field of view, despite the ENS displaying well established directionality along the longitudinal axis 43,44 . Future improvements in imaging and acquisition will, however, provide further opportunities to utilise AI-driven approaches to reveal yet more detail in terms of network structure and pattern recognition in both development and disease settings. We conclude that this study provides critical evidence of occult changes in the ENS, including structural and molecular remodelling, upon loss of nNOS signalling. Further, we propose that novel machine learning

Methods
Animals. 6 HuC/D + neuronal quantification. HuC/D positive neuronal cell numbers were quantified, in a blinded fashion, using the "Analyze Particles" function in FIJI software 45 . Briefly, following auto-thresholding, five independent images (1 μm optical section; 512 × 512 pixels) of HuC/D + cells, at the level of the myenteric plexus, were analysed in each mouse distal colon in an automated fashion using the following settings: Size (μm 2 ) 10-Infinity; Circularity 0.00-1.00. After unblinding, average neuronal cell counts were calculated, per individual mouse, to allow inter-and intra-group comparisons.
Neuronal network image acquisition. To allow quantitative analysis of enteric neuronal networks, five Z-series image stacks (21 × 1 μm optical sections; 512 × 512 pixels) were obtained from the distal colon in a cross-shaped array using a 40× (NA1.2) water immersion objective (Zeiss, Germany). Initially, a Z-series image stack was acquired centrally, approximately 2 cm from the internal anal sphincter, followed by four subsequent Z-series stacks obtained at a distance of 1000 μm in the proximal, distal and both lateral directions using automated directional functions in Zen software (Zeiss, Germany). Confocal micrographs were digital composites of Z-series scans. Final images were constructed using FIJI software 45 .
Neuronal network analysis. Each Z-series (21 slices) underwent a pre-filtering process and conversion to a single frame Z-project. Five such Z-projects were produced for six mice of each cohort (i.e., 30 z-projected hpf). For each mouse, z-projects were stacked and segmented using ImageJ's (version 1.53c) 'Trainable Weka Segmentation' and a classifier that was previously trained on similar images. The resulting segmentation was thresholded and/or skeletonised before standard analysis operations within ImageJ were applied, and the results exported to Excel for analysis. Analysis metrics, based on key, previously reported, characteristics of neuronal networks were assessed, including TuJ1 density, interganglionic area, and orientation. Further, to enable robust branch and junction analyses, 3D Ridge Detection was applied to produce a 3D skeletonised network which allowed inclusion of Z-plane branching in subsequent analysis. A detailed description of the processes applied can be found in Supplementary Information. qRT-PCR. RNA was extracted from C57BL/6J and Nos1 −/− tissues using TRIzol reagent (ThermoFisher, UK) and treated with DNase I (Qiagen, UK). First-strand cDNA was amplified from 1 μg RNA using SuperScript VILO cDNA Synthesis Kit (ThermoFisher, UK). qRT-PCR was performed with an ABI Prism 7500 sequence detection system (Applied Biosystems, UK) using the Quantitect SYBR Green PCR kit (Qiagen, UK), according to the manufacturer's instructions. qRT-PCR was performed in triplicate, using region-specific primers designed against mouse sequences for Gapdh, Chat, Nos1, Sub P and Vip (Supplementary Table 3). Gene expression data were expressed as a proportion of Gapdh, as a reference, using ΔΔCT calculations.
Statistical analysis. Data are expressed as mean ± standard error of the mean. Statistical analysis was performed using GraphPad Prism software (GraphPad, USA). For quantitative real time PCR (qRT-PCR) analyses fold-change comparison between control and Nos1 −/− samples was performed using ΔΔCT and statistical