Authors: Liam McLaughlin, Bo Zhang, Siddharth Sharma, Amanda L. Knoten, Madhurima Kaushal, Jeffrey M. Purkerson, Heidie L. Huyck, Gloria S. Pryhuber, Joseph P. Gaut, Sanjay Jain
Categories: Article, Light-sheet microscopy, Nephrons, Kidney diseases
Source: Nature Communications
Authors: Liam McLaughlin, Bo Zhang, Siddharth Sharma, Amanda L. Knoten, Madhurima Kaushal, Jeffrey M. Purkerson, Heidie L. Huyck, Gloria S. Pryhuber, Joseph P. Gaut, Sanjay Jain
The human kidney maintains homeostasis through a complex network of up to a million nephrons, its fundamental tissue units. Using innovative tissue processing and light sheet fluorescence microscopy, we mapped the 3D neurovascular connectivity of nephrons to understand how their structural organization enables coordinated functions like filtration, absorption, and blood pressure regulation. Our analysis revealed developmental changes in glomerular orientation, density, volume, and innervation from birth through aging. We discovered an extensive nerve network connecting different nephron segments and organizing glomeruli into distinct communities. These communities are linked through “mother glomeruli” that serve as control centers, creating a repeating pattern throughout the cortex. This sophisticated neural organization, which is underdeveloped in newborn kidneys and disrupted in conditions like diabetes and hydronephrosis, appears to facilitate synchronized responses to maintain fluid balance. The findings provide insights into how the kidney’s structural architecture enables coordinated function across its numerous nephrons.
The functions of organs and tissues are dependent on precise structural organization that begins during development and is maintained in adulthood. This is critical for vital organs such as the human kidney where there is a complex organization of the main functional tissue unit (FTU), the nephron, across a cortex and medulla in association with surrounding neuro and vascular networks to maintain fluid homeostasis, blood pressure, and other endocrine activities^1^. Various parts of the nephron and their associations with surrounding structures have distinct functions, including filtration by the glomerulus, regulation of blood flow and pressure by the juxtaglomerular apparatus (JGA), absorption and secretion of various solutes by the tubules, and water balance by the collecting system facilitated by a complex array of cell types^2,3^. These essential roles in fluid homeostasis and blood pressure regulation also render nephron segments as targets for medications to control fluid balance and blood pressure^4^. Model system studies have demonstrated complex positive and negative feedback loops between renal sensory and sympathetic nerves through their connections with different parts of the nephron to regulate renal blood flow and salt balance, and in kidney dysfunction, renal denervation can alleviate fibrosis^5^. Tubuloglomerular feedback and arteriolar myogenic activity at the glomerular vascular pole are ways to synchronize blood flow in an individual nephron. Elegant studies in rodent models using microscopy techniques to detect blood flow changes have shown that changes in blood flow in a single nephron are remarkably synchronized among clusters of nephrons in the same arteriolar network and also between more distant communities^6^. Developmental studies have begun to show evidence that sensory and sympathetic nerves follow the blood vessel in the kidney at a macroscopic level to their respective vascular poles of the glomeruli^7^. There at the JGA, macula densa (MD) cells of the distal nephron are sensitive to variations in sodium and chloride ion concentrations in the filtrate due to proximal tubular dysfunction. These changes can modulate the activity of the renin angiotensin aldosterone system (RAAS) by engaging several processes, including relaying calcium signal waves and other chemical mediators through the extraglomerular mesangium and the activity of different cells in the JGA, including interactions of renin-producing granular cells with renal nerves and adjust glomerular blood flow. However, it is unclear how millions of nephrons coordinate these individual responses from each nephron in the human kidney, nor what anatomical and structural patterning is necessary to enable maturity of kidney function from postnatal to adult kidneys and how this organization changes with declining kidney function with age and disease.
Here, we show the structural connectivity and organization of blood vessels and nerves (neurovasculature) with the nephron segments, and they change across the lifespan and in disease. We developed 3D multiplex immunofluorescence methods and analytical pipelines using light sheet fluorescence microscopy (LSFM) on thick human kidney slices across the cortical-medullary axis and at various ages for the human biomolecular Atlas program (HuBMAP)^8^. We report that glomerular attain a highly polarized orientation with increasing nerve density around them as development progresses. They organize into communities that are interconnected through nerves. The 3D image analyses revealed both sensory and sympathetic nerves contribute to the neural networks. The afferent/efferent nerves only distinctly diverge beneath the MD, which showed only sympathetic innervation at the JGA, but had sensory innervation along the MD’s flanks. Neural connectivity was observed within multiple segments of the same nephron, such as from the glomerulus to its thick ascending limb (TAL) and proximal convoluted tubule (PCT), between glomeruli, and between glomeruli and segments of other nephrons. Importantly, the glomerular communities showed several patterns of innervation and most notably, demonstrated adjacent glomerular communities connected via inter-communal glomeruli that we deemed “mother gloms,” based on the concept of highly connected mother tree hub points^9^. The communities were interconnected with mother glomeruli throughout the cortex in a repetitive pattern, thus providing a structural organization that is poised to sense or relay signals throughout the nephron and that we hypothesize is engineered to coordinate physiological responses by decreasing the burden on individual nephrons. The network was severely diminished in diabetic patients or obstructive nephropathy samples with advanced kidney disease, suggesting a visceral nephropathy in chronic kidney disease.
Architecturally complex solid tissue, such as kidney posit, challenges for 3D volumetric imaging, impacting tissue penetration of reagents, clearing methods, imaging, and analysis. We tested several tissue clearing methods including passive CLARITY^10^, iDISCO^11^, CUBIC^12^, and CLARITY SHIELD^13^ on fixed thick tissue slices from nephrectomies or deceased donors with CLARITY SHIELD as the method of choice here (Fig. 1). Our goal was to define anatomical organization of select FTUs along the nephron in relation to neuronal patterning, morphology at macro and microscopic scales and across various ages ranging from neonatal to aged individuals and in some cases samples with disease (Supplementary Data 1). We examined the overall organization at low power (5X) and then select field of views (FOVs) in 3D at high power (20X) with depth ranging from 1 to 3 mm. In general, kidney structures / cell types evaluated in 3D included nerves, vessels, glomeruli (gloms, podocytes immunostained for NPHS1 or PODXL), collecting duct (CD) (CD, principal cells immunostained for AQP2), proximal tubules (PT) (PT, PT cells immunostained for LRP2) and TAL (TAL, TAL cell immunostained for UMOD) of the loop of Henle (LOH)) (Supplementary Data 2). For nerves, TUJ1 (targeting TUBB3) as a pan neuronal marker, and for sensory and sympathetic nerves CGRP and TH, respectively, were used. For endothelial cells in the blood vessels CD31 or CD34 or PODXL were used. In total, 23 samples across 15 patients were imaged using LSFM (Supplementary Data 1) and 22 samples analyzed for this study. A summary of the pipeline is displayed in Fig. 1.Fig. 1End-to-end pipeline for light sheet fluorescence microscopy from tissue processing to imaging to analysis.1–4 Tissue processing and immunofluorescence staining of corticomedullary strips of tissue slices from kidney, preserving in fixative, active clearing using and performing multiplex immunofluorescence staining with antibodies targeting key functional tissue units of the kidney and / or neurovasculature. The staining was typically done for five targets across the Cy5, Cy3, 488, and 405 channels. 5 LSFM imaging. Cleared and post IF samples were imaged using Zeiss LSFM to obtain gross macroscopic view of the entire specimen at 5X followed by a series of high-resolution, 20x images to visualize subregions along cortico-medullary axis. These 20x images were typically 1 × 1 × 1 mm, although a small number were larger. 6 Data processing. Prior to any analysis, data was stitched and downsampled by 4x in both X and Y using either ZEN Blue 3.1 software or Stitchy. The purpose of this pipeline, as opposed to imaging at low resolution and skipping downsampling, was due to our access to sufficient data storage and processing resources. This enabled us to feasibly obtain high resolution data that could be stored and used for validation and display, while downsampled data could be used for analysis. 7 Data were segmented for structures of interest utilizing either Cellpose or Labkit. 8 Qualitative and quantitative data analysis for various structures and their relationships was accomplished using commercial software such as Imaris or ImageJ, in addition to in-house Python code. This included coregistration of 20x subregions onto 5x volumes, and algorithmic extraction of neural networks between glomeruli. Cartoon depictions from 1, 4, and 8 were obtained using Biorender. Created in BioRender. Jain, S. (2025) https://BioRender.com/z62w451.
Much of our knowledge about kidney organization comes from 2D sections that are limited in providing relationships between different nephron segments and extracellular material, including their connectivity. To gain insights into structure-function beyond 3D vascular networks, we first evaluated neural network in relation to human kidney FTUs. The 3D volumetric imaging using LSFM with the indicated markers from a region of the kidney lobe along the corticomedullary axis revealed dense innervation of the cortex and medulla, with nerves following the medium and small sized vessels, alongside successful visualization of all labeled structures and their 3D relationships (Fig. 2a, supplementary movie 1, 5x).Fig. 2Three-dimensional morphometric analysis by LSFM.a 3D LSFM image showing cortico-medullary organization in a 46 y/o male kidney (Sample SK3). Segmented glomeruli (magenta spheres) in cortex, vessels (green) and nerves (red filaments) entering at cortico-medullary junction (CMJ), with nerves following vessels to glomerular vascular poles or vasa recta; collecting ducts (purple) extend from outer cortex to medullary tip (Supplementary Movie 1). b LSFM of neonatal kidney showing glomeruli (magenta), nerves (red/yellow), and medullary rays (green). Neural branches from arcuate-artery bundles extend upward through cortical lobules, connecting through glomeruli at various depths (arrows). c Radial heatmaps of glomerular vascular pole orientation (n = 96 glomeruli; 1 F; Supplementary Data 3). Azimuth plane shows gaussian two-quadrant bias independent of depth; polar plane displays bimodal bias with vascular poles oriented toward kidney capsule (Supplementary Data 3). d JGA innervation (20x): TUJ1+ nerves (red) innervate entire macula densa (MD); CGRP+ sensory nerves (green) innervate flanking MD cells (Supplementary Movie 3). e–g FTU morphometry from outer cortex to medulla superficial glomeruli densely packed with corticomedullary glomeruli in columns (n = 11 samples; 4 F/3 M; Supplementary Data 4); collecting duct density increasing toward CMJ then exponentially in deep medulla (n = 8 samples; 4 F/2 M; Supplementary Data 5); peri-glomerular nerve density showing steady distribution after initial cortical rise and CMJ decrease (n = 6 samples; 1 F/3 M; Supplementary Data 6), confirmed by 20x measurements. h 20x volume showing glomeruli (magenta), collecting duct (purple), and TH+ sympathetic nerves (red). Example nerve path (1–3): interlobular vessel→afferent arteriole/glomerular pole→efferent arteriole/medullary ray. i Schematic of FTU neural connectivity. Sympathetic (red) and sensory (blue) nerves follow AA to JGA, with sympathetic targeting MD. Both types project along Bowman’s capsule to other glomeruli, cortical labyrinth, vessel-associated bundles, and medullary ray structures (predominantly PT and TAL). Post-glomerular projections extend to medullary ray, with independent cortical labyrinth innervation also observed (Supplementary Movie 3). j–m Medullary neurovascular relationships (Supplementary Movie 4): J Preferential nerve association with vasa recta versus collecting ducts. k PT and TAL coursing among vasa-recta/nerve columns with prominent innervation. l Detailed PT innervation. m Detailed TAL innervation. Source data are provided as a Source Data file.
As expected, the arcuate vessels entered the cortico-medullary junction (CMJ) and then projected interlobular arteries into the cortex that branched laterally into afferent arterioles (AA) in a step-wise manner from CMJ to outer cortex and ending at the vascular pole of each glomerulus and often branched to connect vascular poles of several closely associated glomeruli (Fig. 2a, supplementary movie 2). The nerves also followed these AA that connected several closely located glomerular communities and a common intralobular artery, consistent with recent observations^14^. The AA entered the glomerulus as glomerular capillaries and exited as efferent arterioles (EA) that further branched into capillaries surrounding the tubules or transitioned into the vasa recta in the medulla. The nerve trunks, consisting of both sensory and sympathetic nerves, follow this vascular pattern to the vascular pole, consistent with previous reports of the vascular pole as a heavily innervated site. However, we note that the distinct neural trees emerge from the arcuate vessels, tracking along the interlobular vessels, form horizontal connections between adjacent neural trees along rows of glomeruli laterally interspersed between the medullary rays (Fig. 2b), suggesting secondary organization of neural activity on a macroscale. We next examined the anatomy, structure and cell type relationships of the vascular pole of the glomeruli as it harbors JGA^15^, an important FTU in the kidney consisting of terminal AA, EA, extramesangial glomerular cells and MD; MD apposes the AA-EA with intervening extraglomerular mesangium^16^. The JGA regulates blood flow to the glomerulus and hence filtration and operationalizes highly synchronized regulation of tubuloglomerular feedback, RAAS and positive and negative feedback between the sensory and sympathetic nerves. We first examined whether the relationship of the topology of arteriole entry into the glomerulus was random or instead followed a specific orientation that fixed the JGA in a particular orientation. We developed an algorithm to measure spherical coordinates of angular orientations of glomeruli with respect to their vascular poles from outer cortex to the CMJ in azimuth and polar axis. We found a two-quadrant, Gaussian directional bias in the azimuth axis, and a two-quadrant, bimodal directional bias in the polar axis where vascular poles were usually oriented upwards towards the kidney capsule at shallow angles throughout the cortex (Fig. 2c, Supplementary Data 3). The bimodal polar orientations demonstrated that glomeruli finish maturity in one of two general orientations that may be hypothetically explained by them forming on opposite sides of a ureteric bud, resulting in two sets of glomeruli that point in opposite directions^17^. Likewise, the Gaussian azimuth bias for the third and fourth quadrants may be a result of directionally-biased vascular development in the kidney. Since vessel structures arising from the renal artery must enter kidney regions from the same initial direction^18^, glomeruli may be oriented towards these vascular inputs, inducing an angular bias.
We next examined 3D relationship of nerves with the MD and glomerulus. The nerves, upon reaching the terminal end of the AA at the vascular pole, bifurcated into two or more fibers, none of which follow the AA into the glomerulus itself. In addition to connections with the vasculature, the nerves pass beneath the MD at the JGA where Synapsin I staining revealed a high-density of synaptic varicosities indicating contact points with the MD (Supplementary Fig. S1j). Another branch of these fibers moved on from the vascular pole and traveled avascularly along the outer layer of the Bowman’s capsule (BC) surface, wrapping the glomerulus like a net, consistent with recent observations^14^, before moving onto additional targets (see later) and avoiding the urinary pole (Supplementary movie 2). While both sensory (CGRP + ) and sympathetic (TH + ) nerves were present at the vascular pole, we noted a distinct pattern of MD innervation by each (Fig. 2d, Supplementary movie 2, z-axis 2D views). The TH+ fibers innervated all MD cells spanning the MD, while CGRP nerves only targeted the flanking MD cells of the MD region, indicating a basis for some of the sensory-sympathetic feedback mechanisms in response to alterations in tubular function and blood flow and regulation of the RAAS16,19. These may be signaling through a number of synaptic proteins recently identified in the MD^20^.
Quantitative analysis revealed a decrease in glomerular density from outer cortex to corticomedullary junction where converging medullary rays occupied more space^21^ (Fig. 2e, Supplementary Data 4). CD converging radially gradually occupied greater density (Fig. 2f, Supplementary Data 5) before spiking in density upon reaching the CMJ. This resulted in an overall broccoli-like glomerular distribution, with distinct glomerular columns in the juxtamedullary cortex where converging rays took more space, but dense clumps in the superficial cortex (Supplementary Fig. S1a). This glomerular pattern in the adult kidney is consistent with a rapid increase in nephron number later in human kidney development, where they grow as arcades emanating from the branching ureteric buds^22^. The CD density continues to rise throughout the medulla, consistent with the merging of CDs from multiple medullary rays and renal lobes. In contrast, the nerve densities at different depths of glomeruli did not vary outside of an initial rise in the very superficial cortex (Fig. 2g, Supplementary Data 6). 20x innervation quantification across cortical depths likewise supported the idea that glomerular innervation was largely invariable across cortical depths (Supplementary Fig. S1d, Supplementary Data 7, 8), implicating that the slight rise observed in 5x superficial cortex was a factor of diminishing signal around thinner nerves rather than a true anatomical feature.
We next traced the path of the nerves after their tracking along the AA and innervation of JGA components and the glomerulus (Fig. 2h, Supplementary movie 3). Groups of closely apposed glomeruli often demonstrated several nerve connections to one another after leaving the outer layer of the BC, serving as communities in a neuronephron network. Concerning tubules, branches at the JGA jumped from the AA to the EA, before moving beyond the glomerulus to the medullary ray and innervating other nephron segments, including the PT, CD, and TAL for some distance, but not reaching the CMJ. In the medullary ray, the same nerve branched to innervate adjacent TAL, PT, and CD. Although innervation of TAL was confirmed to be of the same nephron, we could not conclude this about the PT or CD due to their long paths exiting the FOV, however there was some suggestion of this due to the nearby adjacency of the innervated PT and the PCT of the nephron in question, and the connecting tubule of the innervated CD being oriented towards the DCT but exiting the volume before connecting. These data demonstrate that the same nerve connected the glomerular periphery and tubular segments from the same nephron, setting precedence for close communication between these structures (Fig. 2i). We also noted that there was cross nephron tubular segment connectivity by nerves that are a part of glomerular communities described above indicating quite a sophisticated system of cross talk among these FTUs (Supplementary Movie 3). Within the ray, PTs and TALs appeared the most prominent targets for innervation, often receiving parallel-nerve branches, while the CD appeared innervated only at specific points. More nerve targets downstream of glomeruli included short-lasting projections into the cortical labyrinth, including PCTs innervated by branches from its glomeruli and other glomeruli closely linked in the neuroglomerular network, and by nerves that returned to larger nerve bundles moving with interlobular vessels. While glomeruli and their associated nephron structures were the primary target of nerves in the cortex, the nerve bundles also sent numerous short-lasting projections into nearby cortical labyrinth, unassociated with glomeruli, as they progressed through the tissue (supported by Synapsin I and TUBB3 staining). Synapsin I in particular (Supplementary Fig. S1I) emphasized previous studies that showed nerves in the kidney cortex utilized consistent en-passant neuroeffector junctions^23^, together suggesting that nerves in the kidney cortex provided some degree of innervation at all points traversed.
In the medulla, vasa recta projected downward from arcuate vessels in interspersed columns, followed by networks of nerves projecting from the arcuate-associated nerve bundle. Organization of nephron structures within the medulla differed from the medullary rays, with both the proximal tubule and LOH traveling amidst vasa recta-nerve columns, while CDs moved apart from the columns in separate, largely anervous bundles of the long loop nephrons from deep cortex (Fig. 2j, k). Consequently, innervation in the medulla appeared to primarily target the vasa recta, PTs (Fig. 2l) and LOH (Fig. 2m) far more than the CD (Supplementary Movie 4).
The extensive innervation revealed above prompted us to further examine connectivity between nerves and glomerular communities to glean insights into how kidney FTUs may work in a coordinated manner. We use the concept of a network as a collection of objects, known as nodes, and their connections to one another, known as edges. Within networks, cohesive groups called ‘‘communities’’ can be formed, which relate to a system’s functional modules^24^. Groups of nodes organized into communities may together perform distinct, semi-autonomous functions within the wider network^25^.
To discover neuro-glomerular networks from 3D LSFM imaging data, we developed an algorithm where we treated glomeruli as nodes and their neural connections as edges using both manual tracings and automated network detection after the nerves have left the vascular pole (Figs. 1, 3, see methods). We discovered repeating patterns or motifs of glomerular innervation within and across communities (Fig. 3). Network motifs represent network patterns that repeat more frequently than in a random network, which in biological systems are often associated with discrete information processing substructures of a network^26^. We identified seven motifs of neuronephro-networks that we broadly group into intra-community motifs describing basic glom-nerve relationships (Fig. 3a–c), community motifs describing patterns of community aggregation (Fig. 3d, e), and inter-community motifs, which describe the relationships between communities (Fig. 3f, g); illustrated in Fig. 3h and Supplementary Movie 5.Fig. 3Neuro-glomerular network motifs.a–c Motifs existing within a community, d–e Community motifs, and f–h represent community relationships. a A lone glomerulus (blue) receives nerves (yellow) directly projecting from an interlobular vessel (green), without direct neural projections to other glomeruli. Nonparticipating nerves and gloms are shown in red and magenta, respectively. b Paired glomeruli (blue) are closely apposed and are vascularized by an afferent arteriole (green) that forks into two vessels immediately before entering both. They share a tightly linked neural network along their vascular poles and wrapping the outer layer of the bowman’s capsules. c A pyramid motif occurs when one glomeruli (white) appears to be upstream in a network to several other glomeruli (blue). The network graph of the same glomeruli is overlayed atop this, and all following, motif images. d A keychain motif (blue) occurs in low glom density tissue when a branching neural projection from the interlobular vessel spreads between several glomeruli before looping back to its nerve bundle in a circle. Keychain motif appeared as typically five to seven nodes (representing glomeruli) of equivalent degree linked in a ring e A grape motif (solid magenta) occurs in high-glom density tissue when glomeruli oriented around an interlobular vessel surround a centralized nervous hub, whose projections link multiple glomeruli (eight or more). f An hourglass motif (solid-colored glomeruli) occurs when two distinct neuro-glomerular communities which receive distinct vascular inputs are unified by one or a few vertex glomeruli that sit between communities, referred to as mother glomeruli. g Lattices (solid-colored glomeruli) occur when multiple communities are unified by mother glomeruli. h A cartoon depicting several of the neuro-glomerular motifs and their relationships with other nephron structures. Beyond joining glomeruli into networks, nerves may establish additional connections between nephron structures in the cortical labyrinth or medullary ray the letters correspond to the images A–F. See Supplementary movie 5.
We were particularly intrigued by the inter-community motifs. The hourglass (Fig. 3h, Supplementary Movie 5—0:36–1:20) structure was observed when two grape or keychain communities had their nerve networks joined by one or more central glomeruli. Significantly, each community was observed to have a distinct interlobular blood vessel source and community-specific projections into the medullary ray; the central glomeruli had additional projections into the medullary ray from both communities in the hourglass motif. Consequently, these nerve networks and central glomeruli joined together otherwise discrete nephrovascular communities, possibly synchronizing activity between adjacent glomerular communities. In the network, hourglasses appeared as two dense groups of high-degree nodes joined by one or a few nodes that existed between these communities. Attached to the hourglass, additional interglomerular connections could also be seen to join lower-degree nodes that themselves existed as intra-community motifs such as pyramids. Furthermore, analysis of the large 3D volumes revealed the lattice motif when three or more communities joined together by central glomeruli (Fig. 3g, Supplementary Movie 5—2:24–2:50). This structure represented the largest macroscale pattern we could observe in these networks. Lattices were not restricted to linear connections between multiple communities; for example, three communities could be joined to one another in a triangle structure. In dense 5x networks, this often resembled a classic spoke-and-hub network structure. Organizations like this have been described in real world systems where interconnected networks enable information or material flow with processing at central critical nodes (hubs)^27^. This organization is advantageous as it simplifies networks of scale by optimizing hubs, a trait relevant to kidneys that must control up to a million FTUs simultaneously, with a notable weakness being a vulnerability to the disruption of hub nodes. This may suggest that the loss of these central glomeruli is more deleterious to the neurophronetworks as compared to other glomeruli. Given that the community-unifying pathways traverse along central glomeruli, we hypothesize that these nodes may be of greater importance to network function and subsequently define these objects as ‘mother glomeruli’. A summary of all these motifs and their relationships with other nephron structures with the mother glomeruli at the center of an hourglass is illustrated (Fig. 3h).
We next performed a number of analyses to determine the robustness of the neuro-glomerular network across scales at both 5x and 20x resolutions, pattern at different glomerular densities and types of nerves (Fig. 4). To reveal networks across communities at low power, the main nerve trunk to which all structures are connected were removed. One way to describe communities within a network is through its modularity. A modularity closer to 1 is indicative of more statistically ‘‘surprising’’ community structure, with greater intracommunity as opposed to intercommunity connectivity, while a value closer to 0 represents a more random arrangement^25^. Higher-glomerular density samples with abundant nerves exhibited a single, large and highly modular interconnected component, demonstrating an abundance of secondary interglomerular nerve connections beyond the trunk. This component often contained the majority of the glomeruli that were previously present if the trunk was not removed (Supplementary Fig. S2b). Lower-density samples with more diffuse nerves often had many glomeruli released into small components or removed from the network entirely if the trunk was removed, but left a few medium-sized components representing modular hub-and-spoke entities (Fig. 4a, Supplementary Fig. S2c). It subsequently appeared that samples did vary in their community synchronization, with low-glomerulus/nerve density tissue typically possessing a greater degree of more-independent modules while denser cortices had more abundant interconnections. Nonetheless, hub-and-spokes dominated amongst any larger components that remained in all trunkless networks.Fig. 4Characterization of neuro-glomerular networks at both 5x macro and 20x micro scales.a1–3 Network analysis (5X) revealing modular neuro-glomerular connectivity patterns. a1 Segmented image (Sample SK3, PPID 3778) showing spheroid glomeruli and TUJ1+ nerves (red filaments) spanning cortex; largest interconnected network component highlighted with distinct glomerular communities (solid colors) and modularity score (0–1) indicating community structure strength. a2- Connectivity map with individually numbered glomeruli; red circle highlights major component from A1. a3 Enlarged connectivity map of highlighted component with community colors corresponding to A1, demonstrating interconnected modular structure. b Log-log plot comparing relative frequencies of glomerular degree distributions (Sample SK1, PPID 3785) at dual resolutions. 20x distribution (n = 222 glomeruli, blue) shows power-law fit (α = −1.50); 5x distribution (n = 440 glomeruli, red) approximates α = 1.51. Supplementary data 10. Networks exhibit scale-free characteristics with closely matched distributions across resolutions. c Comparative analysis of modularity between reference 20x neuroglomerular networks (orange) and equivalent-sized random networks (blue). Networks display minimal modularity (<10 glomeruli) within single communities, followed by rapid increase between 10–20 glomeruli and continued gradual rise with community expansion. Largest network (222 glomeruli) achieves 0.627 modularity, demonstrating non-random, modular control structure. (n = 26 FOVs**;** 1 F/3 M; Supplementary Data 9). d 20X neuroglomerular connectivity visualization (FOV #8, Sample SK1): segmented glomeruli (magenta), nerves (red), and network connections (yellow) encompassing 370 glomeruli. Trunkless network analysis identified 222 glomeruli in single, highly-modular network. d’ Corresponding 3D-trackable network graph with numbered glomeruli. e CGRP+ sensory neuron connectivity analysis (Sample SK19, PPID 3778, 20X): segmented glomeruli (spheroids) networked via CGRP+ sensory nerves (blue filaments) demonstrating distinctive dual glomerular pyramids joined by mother glomeruli in hourglass configuration; TUJ1+ nerves shown in red. f Schematic Neural trees project from arcuate vessels along interlobular vessels, forming interconnected glomerular neuro-communities. Communities unite via mother glomeruli across cortical regions and different lobules, with connections predominantly occurring along inter-ray lateral pathways. Medullary rays contain parallel-oriented projections downstream of glomerular innervation. Source data are provided as a Source Data file.
Within 20x FOVs, similar patterns emerged; however, the modularity and the strength of community structure varied proportional to the number of glomeruli captured in the image and these observations show a non-random distribution (Fig. 4b–d, Supplementary Data 9). This modularity trend was the result of images with a few glomeruli typically containing only small intracommunity/community-level structures, while 20x images with many glomeruli captured multiple communities in one FOV. Likewise, a sharp increase in modularity can be observed between 20x images containing between ten and 20 glomeruli and also supported in comparative analysis in random networks of equivalent node/edge count. Thus, below 20 glomeruli, samples demonstrated inconsistent community structure that could be less-than-random, but that images with >20 glomeruli samples yielded networks with greater-than-expected community organization. These data therefore suggested that community-level motifs (keychain and grape) contained between 10–20 glomeruli in general, which together formed intercommunity structures (hourglass and lattice) with seemingly systemic rather than a random community organization (Welch’s t test p = 0.013). The intercommunity connectivity diminished as glomeruli became less dense, but these patterns could nonetheless be observed in both dense and diffuse glomeruli tissue so long as a sufficient number of glomeruli were incorporated in the image. Nonetheless, these samples differed in their apparent network involvement, with communities linked at more frequent points in the high-density sample, implicating multiple mother glomeruli, but the low-density sample demonstrating sparser interlinkages containing only one or two mother glomeruli.
Further investigation into the innervation of mother glomeruli (as defined as connecting multiple network communities) in representative volumes demonstrated that mother glomeruli consistently received greater average innervation compared to non-mother glomeruli (Supplementary Fig. S3, Supplementary Movie 8). Region-restricted analysis of glomeruli further showed enrichment of mother glomeruli was explained more by increased enrichment at the vascular pole than the tubular pole emphasizing the likely importance of nerve-JGA interactions.
The observation that the trendlines of glomeruli number and networks at 5X and 20X resolutions from same sample overlap suggest that these are scale-free networks (Fig. 4b). A scale-free network is one not defined by its size or scale, but rather the existence of highly important nodes or hubs^28^, consistent with the spoke and hub structure visible in our networks. A scale free network generally has the property of the fraction of nodes with degree k following a power law distribution p(k) = k^-(α), where α > 1, and 3> α > 2 being more strongly scale-free^29^. Of the reference 5x networks analyzed for this property, seven of nine had degree distributions that satisfied the baseline assumption of α > 1, with all trendlines having r^2 > = 0.8, which also applied to our largest 20x network (Fig. 4d, Supplementary Data 10). It should be noted that only one network satisfied the stronger assumption of 3 > α > 2. Such scale-free, hub-oriented networks carry the important advantage of being strongly resistant to network damage due to accidental or random failures, as the loss of non-critical nodes will scarcely impact scale-free network function as compared to a random network. However, such networks display a critical weakness as the loss of even 1.5–5% of critical nodes can splinter a network into disconnected subgroups. For example, this can be seen in cell metabolic pathways, where mutations in the majority of proteins may hardly affect a cell, but the knockout of the few key metabolic proteins will have a higher chance of causing cell death^28^. On the other hand, a heavy-tailed, inverse power distribution, where a large proportion of nodes were low-degree is suggestive of a complex network more prone to extreme behaviors, such as complete-network activation, when compared to a random arrangement^30^, a factor that would aid any simultaneous management of parallel nephrons.
Finally, we examined if sensory nerves (CGRP + ) contribute to these communities as it would support sensory-sympathetic neural cross-talk across glomeruli and communities in the kidney to modulate neurophysiology of kidney that has largely been implicated to the vascular pole at the JGA and^31^. Manual nerve tracing in these 3D volumes and CGRP network results demonstrate CGRP-stained sensory nerves themselves are sufficient to generate a neuroglomerular network with some level of community structure (Fig. 4e), consistent with our observations of CGRP nerves wrapping and extending further between gloms; TH+ nerves also contributed to the community network, appearing to almost never diverge from Tuj1 stain unlike CGRP (Supplementary Fig. S1g), suggesting a greater abundance of sympathetic nerves overall. A macro-scale drawing summarizing the above neuronephro-networks utilizing mother glomeruli is shown in Fig. 4f.
We next investigated changes in neurovascular nephron organization of the human kidney over life span, including newborn to aged kidneys (1 d neonate, 1mo infant, young adult 18-40 yo, adult 40–75 yo and aged >75 yo), using LSFM to understand how structural changes may impact kidney function (Fig. 5). One advantage of pediatric kidney tissue due to its small size was that it enabled visualization of large areas spanning the whole kidney lobe by LSFM. In the neonate kidney, glomeruli were densely packed, showing a well-organized cortex and a pattern reminiscent of radial nephron development during branching morphogenesis^32^ (Fig. 5a; Supplementary movie 6). In the cortico-medullary axis, glomerular columns were interspersed in columns between medullary rays that merged into clumps as the medullary rays blended in deeper cortex. This is consistent with glomerular arcades forming in deep cortex, where a single ureteric bud branch induces multiple nephrons to rapidly increase nephron number^22^. Interestingly the neonatal glomeruli were bean-shaped compared to spheroidal in the adults (Fig. 5b) and their polar angles lacked the apparent bimodality of the young adult (neonates fail to reject Kolmoglorov-Smirnov normality test; p = 0.08 vs. young adults reject K-S normality test; p = 1.8e-4), possibly explained by the compact packing of neonatal glomeruli that later settle into shallower angles by radial enlargement of the kidney (Fig. 5c). While the vascular networks in the cortex and medulla were well established, the nerve network was visually and quantitively immature relative to the adult kidney. At 1 d, the nerve bundles followed the arcuate artery and branched into the cortex alongside interlobular vessels, but had not yet developed connections to glomeruli or other nephron structures (Supplementary Movie 7). In the medulla, nerves followed vasa recta distal to CDs and did appear to interact with the PT/TAL but with less overall complexity than the adult tissue. The medulla was visually more compact with structures in the outer and inner medulla less defined^33,34^. (Supplementary Movie 4)Fig. 53D Neuroglomerular morphometry across life span.a Neonatal renal lobe (Sample SK14, PPID-3961; 1 d) at 5X (left) and 20X (right). Underdeveloped innervation of dense glomeruli and medullary rays (purple). b Bean-shaped neonatal glomeruli. c Vascular pole polar angles in neonate. (n = 51 glomeruli; 1 M; Supplementary Data 3). d Nephron innervation in infant (Sample SK16, 1 m): renal lobe (5X, left); cortex (20X, right) showing inter-glomerular innervation (white lines). e Young adult kidney (Sample SK2) displaying abundant connections at 5X (left) and 20X (right). f Adult kidney (Sample SK3) showing decreased renal lobe (5X, left) and cortex (20X, right). g Aged kidney (Sample SK12) demonstrating variable glomerular sizes, increased glomerulosclerosis, reduced cortico-medullary slice (5X, left) and cortex (20X, right). h Glomerular density from neonate through adulthood (n = 19 samples; 5 F/7 M; Supplementary Data 1, 11). i Glomerular volumes (20X) across ages. Two-sided Welch’s t p < 0.01 between all groups except adult vs aged. (Neonate n = 1 sample, 166 gloms, Min: 2.94e5, Q1: 4.14e5, Med: 4.69e5, Q3: 5.24e5, Max: 6.31e5; Infant n = 1 sample, 79 gloms, Min:4.85e5, Q1: 8.71e5, Med: 9.57e5, Q3: 1.17e6, Max: 1.53e6; Young Adult n = 2 samples, 179 gloms, Min: 1.13e6, Q1: 1.60e6, Med: 1.77e6, Q3: 1.97e6, Max: 2.36e6; Adult n = 9 samples, 183 gloms, Min: 2.84e6, Q1: 6.51e6, Med: 8.00e6, Q3: 9.10e6, Max: 1.27e7; Aged n = 2 samples, 36 gloms, Min: 2.70e6, Q1: 3.65e6, Med: 5.49e6, Q3: 7.93e6, Max: 1.39e7). (5 F/5 M; Supplementary Data 12). j Glomerular innervation (5X). Two-sided Welch’s t-test: p < 0.01 between all groups except infant vs young adult (p < 0.05), no significance infant vs adult/aged. (Neonate n = 1 sample, 21014 gloms, Min: 0.0, Q1: 4.20e-3, Med: 5.45e-3, Q3: 7.70e-3, Max: 4.20e-2; Infant n = 1 sample, 7101 gloms, Min: 0.0, Q1: 5.28e-3, Med: 7.06e-3, Q3: 9.08e-3, Max: 6.83e-1; Young Adult n = 2 samples, 1751 gloms, Min: 0.0, Q1: 1.42e-2, Med: 2.35e-2, Q3: 3.73e-2, Max: 6.99e-2; Adult n = 6 samples, 3015 gloms, Min: 0.0, Q1: 1.23e-2, Med: 1.73e-2, Q3: 2.24e-2, Max: 1.01e-1; Aged n = 1 sample, 774 gloms, Min: 0.0, Q1: 4.59e-3, Med: 1.54e-2, Q3: 2.05e-2, Max: 6.06e-2). (2 F/4 M; Supplementary Data 6). Source data are provided as a Source Data file.
As kidney development progressed, the glomerular density decreased, volume increased steadily till young adulthood before exponentially increasing in the adult age group (Supplementary Data 11). In the aged kidneys, there was increased as well as decreased glomerular volume observed likely due to hypertrophy of resident glomeruli to compensate for progressive nephron loss supported by previous studies^35^ (Supplementary Data 12). Importantly, the neuro-glomerular networks and community connectivity were noted at 1mo and are well-established into connected communities and mother glomeruli in the young adult (Fig. 5d–g, Supplementary Fig. S2b). Interestingly, neuroglomerular networks appear more diffuse in adult kidneys, with far lower secondary glomerular connectivity beyond the nerve trunk (Supplementary Fig. S2c, Supplementary Data 13). Large modular structures could be detected but were demonstrably smaller in their glomerular inclusivity compared to similar structures in the young adult. Patterns of innervation and organization were otherwise comparable. Adult glomeruli displayed some signs of early obsolescent sclerosis, where the tuft peels away from the tubular pole^36^. In aged individuals there was a notable heterogeneity in glomerular shapes including small-shrunken sclerotic and oblong that may explain the increased density observed in this age group^36^ (Supplementary Fig. S1f). In addition, the aged group demonstrated a 5x nerve network of high modularity and comparable glomerular incorporation to the young adult (Supplementary Fig. S3), as opposed to the more diffuse adult networks, despite having glomerular density comparable to adults. These observations are supported by quantitative analysis of trends across life span in glomerular density, volume and innervation (Fig. 5h–j, Supplementary Movie 7).
We next investigated if neuroglomerular networks are altered in patients with kidney disease or pathology that may provide insights into structural basis for associated pathophysiology. We examined neuroglomerular networks in sclerotic glomeruli in aged kidneys, glomeruli in hydronephrosis near and away from the hydronephrotic area, and in a patient with diabetes. We discovered three types of alterations in neuroglomerular networks that may be maladaptive (Fig. 6).Fig. 6Altered neuroglomerular networks and morphology with age and disease using light sheet fluorescence microscopy.a–b Sclerotic glomeruli (>75 yo) maintain periglomerular and vascular pole innervation. Global glomerulosclerosis shows only rim nephrin staining (magenta) with persistent nerve fibers (red); nerve path around macula densa (MD) marked by orange arrow; DAPI (blue). c–d 3D connectivity of sclerotic (cyan) versus nonsclerotic (magenta) glomeruli; arrow indicates MD from B. Sclerotic glomeruli maintain network innervation (white lines, D). e–i Hydronephrotic kidney analysis (Sample SK8, PPID3909, 55 yo): Segmented cortex/medulla (e) near dilated calyx (f, yellow line) shows neural connections (yellow/white bubbles) with network graph overlay. High-power views (h,i) of outer cortex and medullary junction (green boxes, E) show reduced connectivity versus references; sclerotic glomeruli (cyan) remain connected. g Glomerular innervation reference samples show constant C_M depth density; hydronephrosis decreases with cortical depth (reference: n = 12 FOV, 612 glomeruli; n = 5 FOV, 42 glomeruli; 1 F/1 M; Supplementary Data 8). j Adult reference (SK18) showing greater network participation. k–m Diabetic kidney (Sample SK9, 63yo) shows reduced connectivity (k), abnormal punctate nerve pattern (L), and some preserved vascular pole innervation (m, yellow trace). n Nerve density comparisons. Reference Gaussian distribution; similar median with bimodality from low-innervation glomeruli; significantly reduced (Two-sided Welch’s t test). (Adult n = 6 samples, 3015 gloms, Min: 0.0, Q1: 1.23e-2, Med: 1.74e-2, Q3: 2.24e-2, Max: 1.01e-1; Aged n = 1 sample, 774 gloms, Min: 0.0, Q1: 4.59e-3, Med: 1.54e-2, Q3: 2.05e-2, Max: 6.06e-2; Diseased n = 3 samples, 1053 gloms, Min: 0.0, Q1: 6.64e-4, Med: 1.9-e-3, Q3: 3.62e-3, Max: 2.70e-2). (2 F/3 M, 6). O Nuclear density versus glomerular volume across lifespan/disease (n = 101 glomeruli, ten samples, 5 F/4 M): age-related decline with compensatory hypertrophy; diabetic glomeruli larger with lower cell-density; hydronephrotic smaller but normal cell density; sclerotic show reduced volume/cell-density (Supplementary Data 12, 15). P Glomerular density and connectivity neonate to young adult shows decreased density but increased connectivity; adults show age-related density/connectivity decrease; aged show paradoxical connectivity increase from sclerosed glomeruli (Neonate n = 1 patient, 1 sample; Infant n = 1 patient; 1 sample; Young Adult n = 1 patient, 2 samples; Adult n = 2 patients, 4 samples; Aged n = 1 patient, 1 sample) (2 F/4 M; Supplementary Data 13). Source data are provided as a Source Data file.
Sclerosed glomeruli do not function but the aged patient kidneys surprisingly maintained substantial innervation of sclerotic glomeruli, incorporating them into their neuroglomerular networks in low density regions (Fig. 6a–d). Even MD associated with the sclerotic glomeruli had nerve connectivity (Fig. 6b, c). Keeping nonfunctional glomeruli in the network could be maladaptive by feeding the apparently-intact network bad information and affect kidney homeostasis. One explanation for this pattern would be that as sclerotic glomeruli accumulate with age, the network responds with additional connectivity as an attempt to over-compensate and also not disrupt the highly connected network by loss of innervation to the sclerotic glomeruli.
We examined neuroglomerular network in outer cortex and inner cortex closer to region of dilated calyx in a kidney with hydronephrosis in comparison to a reference kidney (Fig. 6e–j, Supplementary Data 7). The hydronephrotic network was much reduced compared to the reference kidney, and quantification of glomerular nerve density showed progressive decline from outer cortex to deeper areas closer to the hydronephrotic segment, whereas the neuroglomerular density remained constant throughout the cortex in reference kidneys of similar density. There was increased number of sclerotic glomeruli in deeper cortex but they still retain connectivity to the glomerular network.
We also examined neuroglomerular networks in patients with diabetic kidney disease to gauge if there are innervation abnormalities in the kidney similar to neurodegeneration in neuropathy (Fig. 6k–m). Remarkably, small fibers were observed at the JGA but remaining areas displayed punctate pattern of Tuj1 (TUBB3) staining and no distinct fibers reminiscent of axonal neurodegeneration^37^. To further increase confidence in reduced innervation in diabetic kidneys we performed confocal 2D immunofluorescence on 15 additional patients (five^5^ diabetics with mild CKD with no kidney dysfunction and, five diabetic patients with severe CKD). The patients with severe CKD showed marked reduction in glomerular innervation compared to reference and mild CKD diabetes samples while the mild CKD patients showed an intermediate phenotype (Supplementary Fig. S4a, Supplementary Data 14). An acute condition such as AKI did not exhibit any reduction in neuroglomerular connectivity, suggesting that the reduction in network is correlated with chronic kidney conditions and thus may explain fluctuations in vital and metabolic profiles (Supplementary Fig. S4b).
Manual and automated quantification of the above trends in disease and also across lifespan are depicted in Fig. 6n and o and morphological variations in glomerular density (Supplementary Data 15) and networks in Fig. 6p (Supplementary Data 13).
3D multiplexed molecular imaging poses unique challenges in human solid organs due to complex architecture and matrix. Even if the methods of clearing and immunofluorescence are overcome, the amount of data generated pose significant challenge in analyses to understand structure-function relationships. The application of 3D LSFM to human kidney samples and analytical approaches presented here are key advances technologically and conceptually that are applicable to other organs. Our approach provides insights into how 3D structural organization permits various functions and explain kidney pathophysiology. In addition to confirming overall organization of the main neurovasculature branches and nephron components in previous anatomical studies, we provide findings on how these relationships and anatomical organization change from neonatal to aged time periods, such as the finding that avascular neural connections to cortical nephrons develop postnatally. We further report neuronephron patterning in diseased samples or with abnormal pathology discovering several changes in cortico-medullary anatomy and neural networks. Here, the observed changes in glomerular volumes, density, innervation and neuroconnectivity in cortico-medullary axis of key FTUs, across ages and conditions uncovers how mature function and homeostasis may be maintained in the kidney (see below).
One discovery from our work is the observation that nerves connect glomeruli in a community and across communities through secondary connections encompassing the glomeruli and between glomeruli by both sympathetic and sensory nerves. These further connect to other nephron segments, including various tubular FTUs indicating that there may be unrealized complex feedback mechanisms among nephrons beyond the sensory-sympathetic regulation at the vascular pole^19^. The secondary connections caused by nerves meandering between glomeruli may invoke a coordinated response in a connected community by amplifying or dampening activity compared to a lone glomerulus.
Our discovery and introduction of the concept of mother glomerulus connecting glomerular communities suggest that they may serve as important local control centers, similar to the hubs of spoke-and-hub networks^27^. For example, an increased density of sympathetic fibers may make a mother glom more fine-tuned in responding to and modulating changes in glomerular blood flow and blood pressure. This neural connectivity across the kidney of the glomeruli and the associated JGA may be a way synchronization of blood flow occurs in the kidney by close coordination between the nerves and the vessels and their target FTUs. The MD cells harbor several neuronal features, including maculoopodia^38^ and remarkably synchronize calcium spikes^20^. The mother gloms may control synchronicity between glomeruli of different communities via MD calcium signaling and explain basis for blood flow synchronization across the kidney. This idea is supported by our analyses demonstrating greater average neural enrichment of mother glomeruli that was more a result of enrichment along the vascular pole neighboring the MD, as compared to the tubular pole. Serving as an intersection point for sensory fibers from multiple communities may allow mother glom to serve as a highly sensitive and efficient control centers than if every glom was equally involved in the network^27^. Our data supports a scale-free network as mentioned above indicating that this design in kidney makes it less prone to random failure of nodes in the network where the mother glomeruli serve as key hub points. Being a scale-free network would importantly make the neuronephro-network resistant to the accidental failure^28^ of FTUs due to accumulating sclerosis through age or minor damage, as the loss of less neurally-involved gloms would prove largely inconsequential to network behavior.
These observations lead us to provide interpretations for several clinicopathological observations. For example, the effectiveness of denervation relieving symptoms in diseased kidneys^5^ indicates these networks may frequently reach a point of maladaptive failure in age and disease. Our images demonstrate perturbation of the neural network in both age and disease, through incorporating sclerotic components within unusually well-connected networks, or under-connected networks where nerves have retreated from diseased regions. A case may be made that the abundance of spokes allows the networks to sufficiently compensate for glomerular loss until random processes have led to the disincorporation, or instead decoherence through over incorporation, of mother glom hub points. When this threshold is crossed, the network would cease to function as intended, with further neural signaling confounded by dysfunctional network components to the point of maladaptation, worsening disease states and creating a vicious cycle that further damages the network. By completely removing the network through denervation, the kidney is reduced to a uninnvervated, auto-regulated state, potentially improving symptoms through the elimination of a bad network. The persistence of symptomatic relief after reinnervation^5^ may consequently be attributed to the new network resetting itself from bad components, reassembling itself around remaining FTUs in a structure of proper functioning.
Our observation of diminishing neuroglomerular connectivity in aged, disease or neonatal samples may provide an explanation for changes in eGFR and relation to hypertension due to variations in nephron number. With evidence of fewer glomeruli at birth being associated with diseases like hypertension^39,40^, this does raise a question of whether a similarly diffuse neural network is a factor for long-term kidney outcomes. It may also be the case that low-density glomeruli simply do not require as extensive a network to regulate due to their gross reduction of parallelized components, where diffuse networks would be hypothetically appropriate as glomeruli become less dense with age^36^.
We recognize several limitations of our study. 3D imaging methods for solid organs are still being improved and clearing methods in conjunction with LSFM are limited in multiplexity, have long experimental time exceeding 3 weeks and show variability in antibody penetration and staining. Note that 3D LSFM imaging of human solid organs also require meticulous manual validation of each labeled structure for segmentation and masking of artifacts to train machine learning algorithms for automated analysis and reconfirm labeling patterns and segmentation; an example of our manual-automated iterative learning results are shown (Supplementary Movie 9). Further, uniform deeper tissue imaging is challenging with increasing depth in LS microscopes. While we present several observations based on thousands of nephrons, often from two or more samples per patient, in future expanding this to more diverse group of patients is needed. Nevertheless, we were able to discover important trends across age trajectory for several of the observations several consistent with previous works.
Overall, the 3D maps of the human kidney in micro and macroscales of several key FTUs and their relationship to neurovasculature and organization into neighborhoods and communities and transitions at key postnatal time point in health and disease provide previously unappreciated insights into how structure informs on function and its impact on homeostasis, disease and therapy.
For 3D LSFM imaging experiments, antibodies selected were commercially available, validated by confocal microscopy, several FOVs were imaged from each sample and in several cases more than one sample was used for same patient (Supplementary Data 1). For adult kidneys, two or more patients were in each age group with >5 FTUs analyzed.
Sample sizes were not predetermined by statistical methods due to nature of this study. The strength lies in the number of individuals analyzed for 3D volumetric imaging, technologies represented for orthogonal validation, and number of features analyzed.
The following statistical tests were used to establish Two sample Welch’s t test for all comparisons with n > 5 was used to establish significance between cohorts. Welch’s t test is an alternative to Student’s t test that offers greater power when population variances cannot be assumed to be equal^41^.Two sample Wilcoxon rank-sums test was used when n = 5 (only the confocal datasets)^42^, due to offering greater power when t test assumptions of normality and variance cannot be met due to low n.Kolmoglorov–Smirnov test for normality was used to assess possible non-normality of distributions n ≥ 50^43^.
The research complies with all ethical regulation. All experiments on human samples followed all relevant guidelines and regulations. For pediatric samples, human tissue was obtained through the non-profit United Network for Organ Sharing facilitated by the International Institute for the Advancement of Medicine and processed into the University of Rochester LungMAP.BRINDL biorepository before transfer to the Kidney Translational Research Center (KTRC) for use in the Pediatric Center for Excellence in Nephrology at Washington University. Next of kin authorization (parent or guardian) for consent for research was obtained for all donations. The project is overseen by the University of Rochester Research Subjects Review Board as non-human subjects research protocol authorized for the use of the de-identified tissue from the deceased donor (RSRB00056775). Adult reference kidney samples as part of the HuBMAP consortium and disease kidneys were collected by the KTRC under a protocol approved by the Washington University Institutional Review Board (IRB #201102312). Informed consent was obtained for the use of data and samples for all participants at Washington University, including living patients undergoing partial or total nephrectomy or from discarded deceased kidney donors. The study was not designed for studying demographic differences including male and female sex thus no analysis was done to this end. Sex was not used as selection criteria. However, both male and females sources of tissue are included as they became available for the study and shown in Supplemental data Table1. No sex-based analysis was done due to lack of numbers across age-span. Sex of participants was determined based on self-report or from medical records.
The process for procuring pediatric tissue has been reported. Briefly, the kidney was handled as transplant quality, stabilized in transplant buffer solution, and referred for research as unable to place for transplantation. Adult kidney samples were obtained from deceased or surgical nephrectomy donors as reported {10.17504/protocols.io.8epv5r64jg1b/v1; 10.17504/protocols.io.261ge56zdg47/v1}. An ~2 cm ×3–5 mm × 1–3 mm slice of kidney was dissected from a renal lobe, typically containing both cortex and medulla and fixed with 4% paraformaldehyde (PFA) overnight at 4 °C and stored in phosphate buffered saline (PBS) with preservative until use. Briefly, kidney was handled as transplant quality, stabilized in transplant buffer solution, and referred for research as unable to place for transplantation. Adult kidney samples were obtained from deceased or surgical nephrectomy donors as reported {10.17504/protocols.io.8epv5r64jg1b/v1; 10.17504/protocols.io.261ge56zdg47/v1}. An ~2 cm × 3–5 mm × 1–3 mm slice of kidney was dissected from a renal lobe, typically containing both cortex and medulla and fixed with 4% PFA overnight at 4 °C and stored in PBS with preservative until use. All imaging, QC, and data analysis were performed without grouping, and in the same manner regardless of age or disease. Due to this, comparative analysis required knowing tissue disease or age status, and thus blinding was not used. Samples were used as soon as they became available, rather than being grouped into cohorts beforehand.
A detailed, step-by-step procedure for clearing through imaging is available 10.17504/protocols.io.eq2lyjxmqlx9/v1
After fixation, samples were cleared using the CLARITY SHIELD Active clearing protocol^13^. First, samples were immersed in 20 mL of SHIELD-OFF solution (10 mL LifeCanvas SHIELD-epoxy, 5 mL LifeCanvas SHIELD-buffer, 5 mL de-I water) for 3 days at 4 °C. Next, samples were immersed in 20 mL of LifeCanvas SHIELD-ON reagent with shaking for one day. Finally, samples were immersed in LifeCanvas delipidation buffer and underwent active clearing using the LifeCanvas Smartbatch+ delipidation chamber, enhancing delipdation effectiveness and speed. The samples were washed twice with in 1x PBS with rocking for three hours each, then once more overnight. Washed samples were blocked with 1%BSA/0.2%skim milk/0.3% Triton X-100 in 1x PBS (PBS-BB)/0.1% Sodium Azide at room temperature (RT) for 2–3 days with rotation.
Details of the antibodies used are provided in Supplementary Data 1 and 2. Glomeruli were stained with the anti-NPHS1, an antibody marker for nephrin that labels podocytes in the glomerular tuft^44^, or anti-PODXL, which labels podocytes and glomerular capillaries among other vessels. CDs were stained with anti-AQP2, labels principal cells in the CD. Nerves were stained with anti-Tuj1, an antibody that labels TUBB3. Endothelial cells of the blood vessels were stained with the anti-CD31 or anti-CD34 antibodies. Sympathetic nerves were stained with the anti-TH labeling tyrosine hydroxylase. Sensory nerves were stained with anti-CGRP. Synaptic termini for the nerves were stained with anti-Synapsin 1. TAL was stained with anti-UMOD. Proximal tubule was stained with anti-LRP2. Primary antibodies were used at 100 or 200 dilution. Secondary antibodies used were anti-mouse alexa-488, anti-rabbit-cy3 and anti-goat-cy5 and anti-sheep-cy5 and at 150 dilution. Nuclei were stained with DAPI at 1000 dilution.
Samples were incubated with primary antibodies in PBS-BB with rotation either at RT for 7 days at RT or at 37 °C for 2 days, then at RT overnight. Samples were removed from primary and washed with 0.3% Triton X-100 in 1x PBS/0.1% Sodium Azide at RT with rotation overnight with several changes followed by incubation with fluorescence-conjugated secondary antibodies in 0.3% Triton X-100 in 1x PBS/0.1% Sodium Azide at RT with rotation for 2–4 days. Next, samples were washed with 0.3% Triton X-100 in 1X PBS/0.1% Sodium Azide at RT with rotation with several changes overnight at RT. Samples were incubated with DAPI 1000 in 0.3% Triton X-100 in 1X PBS/0.1% Sodium Azide at RT for one day to label nuclei. Finally, samples were washed with 0.3% Triton X-100 in 1x PBS/0.1% Sodium Azide at RT overnight at RT; DAPI could also be included during incubation with the secondary antibody. After staining, the sample were stored in the dark, wrapped in aluminum foil in 1x PBS/0.1% Sodium Azide prior to imaging. We typically previewed the sample using confocal microscopy before LSFM. When ready for LSFM, the sample was refractive index matched with LifeCanvas EasyIndex. The sample was transferred into a new 5 mL Eppendorf tube and incubated sample at RT in a 50/50 mixture of EasyIndex solution for one day. Next, the sample was incubated in EasyIndex at RT for one day.
All samples were imaged on the ZEISS Lightsheet 7. A small amount of super glue was used to affix the sample to a glass coverslip, which was mounted in an imaging chamber filled with 1.52 refractive-index Cargille Oil. First, the entirety of the sample was imaged with a 5x objective with a zoom ranging between 0.4 and 1.0x to capture the entire sample, depending on sample size. The 5x image was used to identify several regions of interest that were imaged with a 20x resolution lens and dual 10x illuminators at 1.0 zoom, using a 2 × 2 tile grid and imaging 1 mm deep in Z-axis, resulting in 20x FOVs around 1 mm^3^ on average. All images thus specified were taken using Cy5, Cy3, 488, and 405 laser lines. Additional large 20x FOVs up to 8 mm^3^ were imaged on select samples to evaluate the relationships of nerves and nephrons at mesoscale at high resolutions; here, data were instead imaged on two or three of the four laser lines to facilitate data processing and make these large FOVs conducive for analysis. Key imaging parameters are summarized in Supplementary Data 17 and individual image-associated metadata are available in Zenodo (see data availability section).
The general pipeline to transform raw image data from the Lightsheet 7 into 3D isosurfaces involved the following steps.StitchingThe Lightsheet 7 microscope captures large volumes by creating overlapping Z-stacks, requiring stitching into a singular image. For 5x volumes and 20x volumes below 500 GB, the stitching was with the stitching function in ZEN Blue 3.3. For 5x volumes and 20x volumes above 500 GB, the stitching software Stitchy was used to stitch the volumes.Down samplingTo enable efficient data processing and analysis, both 5x and 20x volumes were down-sampled 4x in both the X and Y dimensions, resulting in a 16x down-sampling. For images stitched in ZEN Blue, the ZEN downsampling algorithm was applied on stitched files by setting pixel sizes in X and Y to 0.25. For images stitched with Stitchy, down-sampling occurred concurrently, generating a secondary file with X and Y pixel sizes expanded four times.MaskingSegmentation through the creation of binary masks reduces an image to solely positive or negative signal, enabling efficient computation. For this study, binary masks were generated utilizing either Cellpose2 or Cellpose 3^45,46^ or the ImageJ extension Labkit^47^.Cellpose is a trainable machine learning software that allows for custom model segmentation. Models were constructed by manually segmenting 2D slices from various LSFM datasets. A model to segment anti-NPHS1 labeled Glomeruli was constructed from 179 images across six 20x datasets, and 16 images across four 5x datasets. A model to segment anti-AQP2 labeled CDs was constructed from 46 images across three 20x datasets. The model to segment DAPI-labeled nuclei was constructed from three images across two datasets. These models were trained from scratch based on the dataset with 10,000 training epochs each. Importantly, these models were retrained as needed whenever they inadequately segmented an image, receiving additional manually segmented datasets from the new image to ensure accurate segmentation.Labkit, an extension for ImageJ, allowed trainable machine learning segmentation of heterogenous objects in lightsheet volumes. A custom model was trained from scratch for each FOV segmented. Labkit was trained to segment channels containing filamental objects including Anti-Tuj1 labeled nerves for 5x and 20x volumes, anti-CD31 or anti-CD34 labeled vessels for 5x and 20x, and labeled tubule structures, such as CDs, PTs, and TALs for 5x.Files created by both cellpose and labkit were saved as grayscale TIFF files. These masked TIFFs were converted to 8 bit binary in FIJI or Python, setting all nonzero signal to the maximum 8 bit value of 255.Thresholding and isosurface generation
To convert masked volumes into isosurfaces in Imaris for both visualization and analysis, binary TIFF masks were first converted to.ims files using the Imaris File Convertor Program, inputting the appropriate X, Y, and Z voxel dimensions. Binary channels were then imported into Imaris overtop the stitched, down-sampled volumes. Since binary channels separate objects into signal or background, all binary channels had the least possible intensity threshold of 0.125 applied to them, removing no pixels from the mask to the isosurface. Surface smoothing was also disabled for all isosurface generation. However, voxel-size based thresholding of objects was applied as needed to remove any small, mislabeled masks from the isosurface. All isosurfaces were manually validated before any analysis.
Glomeruli in 5x were segmented with their 5x Cellpose 2 model and turned into an isosurface in Imaris. An additional step during isosurface generation involved using the Imaris surfaces seed separation method, which allowed the separation of adjacent fused objects by morphological size. Due to the thickness of the light sheet, often empty space regions just adjacent to the tissue yielded a low intensity signal where there should have been none. These glomeruli shadows would sometimes be selected during segmentation as positive signal. To obtain accurate glom counts and volume information, these incorrectly selected objects were manually removed from the resulting surface within the Imaris isosurface. The resulting glomeruli isosurface was used to mask the original binary glom channel. This newly masked channel was exported for analysis.
Tubule structures in 5x were segmented using a Labkit model and converted into an isosurface in Imaris. Since we took advantage of region specificity of glomeruli and CD and stained them with the same secondary antibody, for CDs, NPHS1 glom signal was removed using the 5X glomeruli isosurface to mask the tubule segmentation if necessary. The resulting masked tubule channel was exported as a TIFF image for analysis.
Nerves were segmented using a Labkit model, then traced in Imaris using the Filaments tool, which automatically converted the segmented nerve mask into a more cohesive filament network. The resulting filament network was converted into a binary channel using the Imaris Convert Filaments to Channel Xtension and exported as a TIFF. This represented a refined nerve mask for 5x that was used for network analysis.
Vessels in 5x were segmented using a Labkit model. Due to the flexibility of labkit segmentation, vessels in 5x were either segmented for CD31/CD34 signal, which returned a mask displaying the vessel wall stain, or they were segmented for empty space within the vessels, which returned a mask displaying the entire volumes of larger vessels. The latter approach was used to both better visualize vessels as 3D objects, and to account for inconsistent CD31/CD34 staining. Masks segmented for empty space within vessels were converted into an Imaris isosurface and had any erroneously selected regions, often at the edge of the tissue, manually deleted. Although vessel masks in 5x were only used for visualization and not quantitative purposes due to the inconsistency of CD31/CD34, later PODXL was used as it consistently performed better than CD31/CD34; note PODXL also stains podocytes in the glomeruli.
Basic sample measurements were obtained using Imaris. The volume of an entire sample was obtained using the Imaris surfaces tool. First, a temporary, down-sampled image was created to increase processing speed. Next, a new surface was autogenerated on the 488 channel of the 5x image. With smoothing disabled, this surface was thresholded until the surface enveloped the entirety of the tissue volume but none of the background. The volume of the resulting surface was obtained from Imaris’s surface statistics tool (Supplementary Data 11). The surface was also used to generate a binary mask to represent the spatial properties of the entire sample.
The length and width of the sample were obtained using the Imaris measurement points tool. The cortical and medullary lengths of the sample were also obtained with the measurement points tool, measuring from the top of the cortex to the corticomedullary junction, and from the corticomedullary junction to the base of the medulla.
Glomeruli counts of an entire sample was estimated by 3D watershedding the 5x analysis glomeruli masks. Watershed^48^ is an image processing method for separating fused binary objects. We used the watershedding algorithm “Distance Transform Watershed 3D” from the FIJI plugin MorphoLibJ^49^ a function specialized for separating binary objects in 3D. This method estimated glom count with 97% accuracy in sample SK5 compared to manual count, for example.
Glomeruli density in the cortex specifically was estimated using a script written in Python. This was accomplished by dividing the total number of glomeruli by the volume of their convex hull. Since glomeruli appeared in cortical regions but not medullary ones, a volumetric boundary that surrounded only regions with gloms was necessary to calculate density. A 3D convex hull^50^ was generated around the analysis glom mask using the scipy spatial module ConvexHull. This represented the smallest convex 3D polyhedron that enveloped all binary signal, which was used to represent the volume of all glomeruli-occupied space. Dividing the number of glomeruli by the volume of the convex hull thus yielded the number of glomeruli per unit of 3D space in the cortex (Supplementary Data 11).
The neural innervation of glomeruli was procedurally quantified by skeletonizing^51^ the Tuj1 nerve mask in FIJI, followed by using a Python script to calculate an average, normalized density of nerve skeletons within 10 μm of glomeruli. The Tuj1 binary mask was used over the filaments nerve mask since it represented a more accurate description of the size and shape of the nerve fibers, and responded better to skeletonization, while the nerve mask was better utilized for establishing connectivity specifically. First, Tuj1 masks were skeletonized in FIJI using the plugin Skeletonize (2D/3D). This tiff file was converted to a 3D numpy array alongside the analysis glom masks.
To enable accurate comparative analysis across images taken at differing resolutions, all arrays underwent a resolution normalization step where the scipy zoom() function re-interpolated the glomerulus and nerve-skeleton arrays to all have the same voxel sizes as one another. In essence, low-resolution arrays were stretched until their voxel sizes were the same as the smallest reference array in the dataset. The interpolated nerve skeletons were then skeletonized once more with the skimage skletonize_3d function. This was especially important for 5x images taken at separate zooms, as normalization allowed all nerve skeletons to represent the same magnitude of nerve across all volumes being compared.
The glomerulus masks were then dilated by 10 µm in all dimensions using the OpenCV dilate function. Next, the non-dilated glom array had its indices inverted, turning 255–0 and 0–255, which was multiplied against the dilated glomerulus array to return an array of all regions within 10 µm of a glomerulus. The total pixel count of this search-region array was likewise obtained. In 5x (Supplementary Data 6, 11), the dilated binary glomerulus array was multiplied against the nerve skeleton array to mask out all nerves not within 10 µm of the glomerulus. In 20x (Supplementary Data 7), rather than using dilated glomerulus to mask nerves like 5x, instead the search region itself was multiplied against the nerve skeleton array to mask out both all nerves not within 10 µm of gloms and all signal within the tuft. This differing strategy between 5x and 20x was utilized because the majority of any nerve signal within gloms tufts was noise at 20x, but not at 5x, due to the resolution differences. The number of remaining white pixels in the masked nerve array was then calculated.
The density of nerve skeletons around a glomerulus in the volume was finally calculated by dividing the pixel count obtained for the nerves by the pixel count obtained by the search region. Ultimately, this density was the value used for comparative analysis. This script also returned an estimated true search volume and total glom search volume by multiplying the search volume and glomerulus masks by the normalized 3D voxel sizes. For simple comparison between broad sample categories (such as reference vs diseased), the automated neural data was also aggregated across FOVs by calculating for each an average weighted by the total number of glomeruli analyzed in a sample. First, the total number of glomeruli analyzed in a sample was estimated by dividing the total volume of gloms analyzed in that sample by the average glomeruli volume of that sample calculated from 20x volumes (Supplementary Data 12). This value was used to obtain an aggregate mean of a category by first finding the total number of estimated glomeruli in that category (such as total number of all reference gloms across all samples), then for each sample in said category, dividing that sample’s estimated glomeruli analysis count by the total estimated glom analysis count and multiplying this ratio by that sample’s average glomular nerve density score. Adding up every value for each sample in a category obtained this way produced a single aggregate estimate for the glom-neural density of that category.
This nerve-skeleton strategy was later reapplied with slight modifications to produce the data shown in Supplementary Fig. 4/Supplementary Table S17. Namely, a parallel method was used to isolate individual glomeruli and their neighborhoods in a small subarray that could be dilated in parallel to isolate the nerve skeletons around individual glomeruli, producing a more glomeruli-specific distribution that was used to analyze the differences in innervation magnitudes between separate glomeruli within a volume.
The neural innervation of glomeruli was manually validated by tracing the filaments within 10 µm of gloms using Imaris (Supplementary Data 8). First, only gloms fully contained within an FOV were selected and used to create a binary mask, which was exported as a TIFF. These gloms were dilated by 10 µm in all dimensions using Python (OpenCV) to establish a boundary region. This dilated TIFF was imported back into Imaris. All nerves in the dilated regions were manually traced from the Tuj1 channel in Imaris using the Filaments function. The lengths of nerves around a particular glom were then determined by selecting only that glom’s nearby traced filaments in Imaris and finding the summed segment length in the Imaris statistics tool. The volume of dilated gloms in the FOV was calculated by converting them into an isosurface in Imaris and exporting their volumes through the statistics tool. Innervation density in an FOV was then determined with the following (Total dilated glomerulus vol—glomerulus vol)/Total nerve length of glomeulus.
The density of structures across a tissue sample was obtained by using Python. This was accomplished by slicing through the masked 5x data in the direction of the Y axis and dividing the total voxel count returned by np.count_nonzero() by the voxel count of the binary mask. This outputs for each plane were saved in excel as two columns representing depth in the tissue and the density of that plane, respectively. Glomeruli and CD density distributions sometimes contained extreme density outliers obtained from regions towards the top and the bottom of the tissue totally occupied by a glomerulus (Supplementary Data 4) or duct (Supplementary Data 5). These outliers were removed using the equation Outlier = Q1-(1.5(IQR)); Q3 + (1.5(IQR)), with Q1 representing the data’s first quartile, Q3 its third quartile, and IQR the value of Q3–Q1. Nerve densities were grouped into buckets of n = 5. In most instances, samples were imaged from cortex to medulla in the y-axis. Consequently, these distributions were used to analyze how densities changed through cortical and medulla depths. However, differently oriented images were not used for this specific analysis purpose. Planar distributions for nerve densities (Supplementary Data 6) were used as a surrogate measure general innervation distribution across a sample, as shown in the innervation violin plots of Fig. 5.
20x FOVs were coregistered onto 5x volumes using Imaris. Due to the Lightsheet 7 requiring switching immersion chambers between 5x and 20x imaging sessions, all 20x FOV regions had to be coregistered post-imaging, rather than in-tandem. While imaging 20x volumes, the gross location of the 20x image was first labeled on a 2D image of the 5x volume to act as a rough guide for realigning the 20x. After stitching and downsampling, the 20x FOVs had their 647 channels extracted, converted into an ims file, and imported into Imaris atop the 5x image. Using the labeled 2D 5x image as a reference, Imaris was used to rotate and translate the 20x FOVs into their correct locations within the 5x ims file, using glomeruli as fiducial markers. Since this coregistration was accomplished manually, it did not result in exact pixel-to-pixel correspondence, but did approximately describe values such as the cortical depth of a 20x FOV. After each coregistration, the cortical depth of a 20x FOV was obtained using Imaris’s measurements points tool to measure the depth from the top of the cortex of the tissue slice in 5x to the center of the 20x FOV (Supplementary Data 16).
Segmenting glomeruli 20x glomeruli were roughly segmented in Imaris with the spots tool, then finely segmented in Cellpose, then manually refined using the Imaris surfaces tool. The stitched, downsampled 20x volumes were first converted into Imaris images. Within Imaris, regions containing glomeruli were first manually labeled using the Imaris spots feature, which allowed the masking out of larger spherical regions containing gloms, which represented a rough segmentation to aid in computationally locating gloms. The use of spots proved to increase segmentation accuracy by acting as a queue to the model that a glom would be inside a spot, however this was skipped for certain high-density images that had sufficiently strong signal. The masked glomerulus file was exported from Imaris as a TIFF and segmented using its Cellpose model, with the mask then converted into an isosurface within Imaris. In a few instances where NPHS1 stain was weak, nearby areas of tubules or ducts would occasionally be masked adjacent to a glom mask, causing the two to fuse in the isosurface. These incorrectly segmented regions were manually cut off and deleted using the Imaris surface cut tool. Within the isosurface, any nearby glomeruli that erroneously fused were also cut apart manually along one clipping plane with the surface cut tool. This represented the analysis glom isosurface used to quantify glom morphologies.To accurately characterize glom-nerve networks, a second segmented glomerulus mask was created that represented both an accurate segmentation of gloms, but also had no touching regions between adjacent gloms. To accomplish this, two clipping planes were used to split the junctions of any fused gloms within the analysis isosurface with the Imaris surface cut tool, and the resulting, small in-between region was deleted. This edited isosurface was used to mask over the original 20x binary gloms. The resulting masked glomeruli were exported for analysis. In certain high-density glom volumes, this was instead accomplished algorithmically by first downsampling the 20x volume, a step that all glom files underwent anyway while creating 20x networks, followed by watershedding using MorphoLibJ’s “Distance Transform Watershed 3D”^49^. Segmenting tubule structures CDs were segmented by their Cellpose model and converted into an isosurface in Imaris. PTs and TALs, which did not have a Cellpose model trained, were instead manually segmented within Imaris using the filament tracer tool to trace the inside of the tubule, for visualization and qualitative analysis when relevant. Segmenting nerves Nerves in 20x volumes were segmented using a Labkit model, then traced in Imaris using the Filaments tool, which automatically converted the segmented nerve mask into a more cohesive filament network. The resulting filament network was converted into a binary channel using the Imaris Convert Filaments to Channel Xtension and exported as a TIFF. This represented a refined nerve mask for 20x that was used for network analysis. Segmenting vessels Vessels in 20x volumes were segmented using a Labkit model. These could be converted into either a filament network or isosurface in Imaris; however, due to its inconsistent stain, the vessel segmentation was used only for visualization and not quantification. Segmenting nuclei
Nuclei in 20x volumes were segmented using their Cellpose model. Only binary glomeruli fully contained within the FOV were used to mask the 405 channel in Imaris. The resulting masked channel was exported as a TIFF, which was used to segment for nuclei. As the high-frequency 405 laser did not penetrate deeper tissue as well as the other laser lines, oftentimes only gloms at the surface of an FOV were included in this analysis, while those deeper in the tissue were excluded due to their reduced signal.
Glomerulus volume and sphericity were obtained from the glomerulus isosurface in Imaris. Only glomeruli contained entirely within the 20x volume were selected within the glom analysis isosurface. The volume and sphericity for those gloms were exported from the Imaris statistics tool. (Supplementary Data 12).
Cell counts in a glomerulus were estimated using Cellpose 2 segmentation. The segmented mask obtained from cellpose maintained identities of individual cells by labeling them with separate grayscale values. The number of cells within an individual glom could therefore be obtained by first using Python to label the analysis glom masks with separate grayscale values as well. By iterating through all labeled gloms, each could be converted into a boolean mask that when multiplied against the cell array would return an array only containing all cells within one glom. These cells were then counted by flattening their 3D array, removing any zero values, and converting the resulting list into a unique set whose length was equivalent to the number of unique objects. Cell densities could likewise be obtained by dividing this value by the volume of its corresponding labeled glomerulus object. (Supplementary Data 15).
Glomerulus orientation was categorized by the angular direction in which its vascular pole was oriented in 3D space using Imaris to label the location of vascular poles, followed by using a Python script to measure angles. First, the gloms separated for network analysis had the locations of their vascular poles labeled in 3D space using Imaris’s spots feature. A python script was used to generate a completely white volume of the same dimensions as the channels in Imaris. Next, Imaris’s surface mask feature was used to create binary masks of only the analysis glom isosurfaces fully contained in the FOV. The spots at the vascular poles were also used to mask over the white volume to create a TIFF containing only binary spots denoting the locations of vascular poles. Both of these binary channels were exported as TIFF files and downsampled in FIJI to speed up processing.
Next a Python script was used to obtain angular information about the glomeruli from the exported channels. The script first calculated the centroids of all vascular pole spot markers from the pole tiff and matched them to the nearest glomerulus object in the separated glom mask. The centroid of this glomerulus was calculated, and both pairs of centroids were matched. The line between these centroids was treated as a surrogate for the 3D orientation of these gloms, with the vascular pole spot centroid as the end point and the glom centroid as the base point of a vector. First, the azimuth angle was calculated with the numpy method np.arctan2(dy/dx), while the polar angle was calculated with np.arctan2(dy/dz). Note that np.arctan2() function is programmed to specially accommodate arctan calculation in all four angular quadrants. The resulting angles for each glom/pole pair were exported in an excel file (Supplementary Data 3).
To visualize the orientations of gloms, verify correct centroid matching, and search for common fates, the centroid pairs were further used to draw rays from the glom centroid, through the pole centroid, to the end of the volume. This was done using a Python script that first used the slope of the vector to find the furthest allowable point in 3D space away from the glom centroid, before using the command np. np.linspace () to estimate a set of linear coordinates of the vector within the array and set these indices to white. A TIFF file was exported with all the drawn rays, which could be opened over the original channels in Imaris to easily view exactly how gloms are oriented in 3D (Supplementary Fig. 1Sh).
Manual and automated analyses demonstrated similar trends between FOVs. To examine this, same FOVs were used for both analyses. For example, manual and automated analysis of the same FOVs throughout samples SK1 and SK2 [PPID: 3785] yielded a shallow linear trendline, with the automated trendline a little less positive than the manual one. For further validation, this comparison was normalized by adding the difference between the manual and automated SK1 FOV2 calculations to all the automated measurements. This was necessary to bring automated measurements into a generally similar magnitude range as manual measurements, since automated measurements relied on a different calculation. A two-sample Welch’s t test between these manual and normalized automatic distributions failed to reject the hypothesis that these distributions are different (p = 0.761). While not entirely perfect, the automated pipeline allowed efficient comparison of innervation magnitudes that were representative of true biological structure, trading a small amount of accuracy for a massively increased sample size **(**Supplementary Fig. 1Sd).
To enable this analysis, we developed “NetTracer3D”, a Python tool able to generate large, undirected networks from 3D images, in addition to providing analytical support. While detailed application will be published elsewhere, we describe its operation and outputs below. With highly customizable parameters, NetTracer3D offers several derivatives of its core algorithm to allow network generation for a variety of circumstances, accuracies, and speeds. Two versions of NetTracer3D were developed, utilizing the same core logic, with V2 offering superior customizability and computational performance on larger datasets. The majority of analysis in this study was accomplished using V1, except for the large 20x dataset shown in Fig. 4d, AKI 5x dataset in Supplementary Fig. S3b, and the results in Supplementary Fig. 4 S/Movie 8. Note that for the purposes of this study, V1 and V2 returned largely similar results due to the parameters used
Both versions of NetTracer3D used the same core logic. Glomerulus and nerve masks were designated as nodes and edges in a network, respectively. In general, nerve masks existed as one or a few interconnected objects in the volume that if utilized directly to establish a network would always result in erroneous connections between the majority of nodes. To establish discrete node-node connections, the node masks were first used to discretize the network. This was accomplished by allowing the nodes to extend outward in a 3D-search space and break any edges passing through their space into separately recognized connections. Essentially, a nerve passing over a glomerulus would be split by the glom’s search space into a before and after segment. By differentiating edges in this manner, distinct connections between objects could be recognized within the neural network. Connections could also be made through edges retained within the search spaces of two neighboring nodes. Any connected objects were then saved in the network as a node pair, which could be analyzed using the Python NetworkX package^52^. As this algorithm is designed to function for any generic 3D image containing potential node and edge objects, it additionally provides future analytical use beyond this study.
The processed glom and nerve binary masks were converted to 3D binary numpy arrays in Python. Since the nerves wrap the outer boundary of the BC rather than entering inside the glomerulus, first, the glom array was dilated by 10 µm in all directions to allow each glom to envelop its innervating nerves. After enveloping nearby nerves, these gloms were used to mask out that region, effectively splitting up the nerve network into a series of glom-to-glom neural connections, rather than one single interconnected bundle. This was accomplished by inverting the indices of the dilated glom arrays, marking the dilated gloms as black null regions and everything else as white regions. The inverted array was multiplied against the nerve array, removing any nerve segmentations that came within ~10 µm of gloms. This allowed the nerve network to be broken up into subconnections between gloms, whereby the separated nerves served as edges in a network, and the gloms as the nodes. The resulting array contained the “outer edges”. Due to the fact that nearby gloms sometimes had nerve connections between one another that were completely enveloped by the dilated regions and therefore eliminated from the edge network, a second set of edges was created by multiplying the non-inverted dilated glom array by the nerve array. The resulting array contained the “inner edges”.
To examine subcommunities of gloms in absence of the nerve trunk that provided extensive network innervation, at this point the nerve trunk could be removed. This was accomplished by identifying the largest object in the outer edges network and setting its indices to black. For consistency, all volumes analyzed for their glom-nerve networks were processed twice, both with and without the trunk, to compare the analytical results.
Due to being an earlier implementation, NetTracer3D V1 utilized a greedier algorithm that required the user to perform an intermediate downsampling step at 20x resolution. The purpose of this step was to prepare the glomerulus array for an ID-specific dilation into each other’s search spaces, but also served as necessary to allow nerves whose segmentation was not a complete filament to become reconnected. In V2, this reconnection is accomplished via a dilation algorithm within the edge array, but in V1, downsampling accomplished an equivalent goal, essentially squishing punctated nerve masks back together.
In only the first implementation of the algorithm, V1, both sets of edges were saved as 8-bit TIFF files. Each edge TIFF was opened in FIJI and had a maximally brightening look-up table applied. The resulting image was downsampled in FIJI by five in all dimensions using the scale function with default settings. The downsampled edge images had a maximally brightening LUT applied before they were saved. The reason for downsampling in FIJI instead of Python is that this specific set of steps in FIJI will downsample the binary network, yet prioritize retaining positive indices that represent nerve signal. By comparison, the scipy downsample method will maintain the ratio of positive to negative elements while downsampling by eliminating positive indices, losing network coherence. This FIJI downsampling method effectively allowed regions of the network that had broken up into punctates due to imaging or segmentation artifacts to reconnect to one another without the need of a large dilation/closing step. Both sets of downsampled edges were then opened in Python as 3D binary numpy arrays. Both edge files were labeled, a scipy method that assigns unique numerical identities to all separated objects within the array. Finally, the labeled array containing the outer edges was dilated by a single pixel in all dimensions, which helped further reconnect any broken regions of the network while also allowing it to overlap the dilated gloms that had previously been used to mask over the nerve network. This overlap was crucial in identifying which edge interacted with which gloms, while the single pixel dilation in the downsampled file proved relevant in reconstituting network loss due to imaging artifacts. Several downsampling and dilation parameters were tested, with the above steps best demonstrating a strong ability to accurately identify connections based on comparative observations with the 20x images themselves.
The glomerular array was labeled to give each glom a unique numerical identifier. The glom array was downsampled in all dimensions by 5, resulting in an array with the same dimensions as the downsampled edges. Next, the dilation space around each glom was determined before a subarray around that glom was obained. By dilating the glom within the subarray and utilizing boolean operations, labeled edges that interacted with that glom specifically could be isolated and exported, used to search for which glomeruli interacted with the same segment of nerves.
Setting up for Network Search – 5x V1 Due to the lower resolution of 5x volumes, downsampling and dilating the nerve masks was not used, as nerve masks did not typically display artifact-induced gaps. Samples were first dilated by 10 µm to establish detection neighborhoods using OpenCV. Next, overlapping gloms were separated using MorphoLibJ’s “Distance Transform Watershed 3D”^49^. This method preserved the correct identities of dilated, separated gloms in 5x with 95.4% accuracy versus manual separation in sample SK5. Note that this dilation followed by re-watershedding was no longer necessary in V2 due to the implementation of smart-dilation (see NetTracer3D V2 improvements). In V1, the dilated, watershedded glom masks were separated and labeled with separate grayscale values index by the plugin, and therefore could be used directly to create networks. The dilated glom mask was used to break the network into discrete units in the same manner as for 20x volumes.
For identifying pairwise connections, a computationally efficient method was developed to allow fast processing in large arrays with hundreds, or even thousands, of gloms. This was accomplished by using array transforms to create a 3D mask for both the glom and edge arrays that only contained labeled components in each that overlapped with one another. Using these masks, any components of the arrays that were not overlapping could be set to 0 with one computational step. Next, these arrays were flattened into 1-dimensional lists and had any indices containing a 0 removed, significantly trimming the number of indices to search through. Finally, the shortened edge list could be compared to the glom list. All values in the edge list of the same numerical label represented the same nerve spatially, which when compared to the glom list, could be used to rapidly decipher which gloms were interconnected by a single nerve structure. This algorithm produced the same results as simply comparing all indices of the labeled edges to each corresponding index of the glom array to establish overlap, but with a processing time taking under a minute even in large arrays, as compared to hours with the rote method.
Each sublist within the master list had its pairwise connections extracted into a new list containing every connected glom pair within the network. For example, an edge yielding the connections [10, 11, 12] to the master list would result in the pairs [10, 11] and [11, 12] being added to the ultimate pairwise list. This pairwise list was exported as a two column excel file.
To allow easy identification between glom nodes described by the 2D network and the original 3D image, the numerical identity of each labeled glom was drawn into a new TIFF file. This was accomplished by identifying the centroid coordinate of each glom in the downsampled array, then drawing the labeled numerical identifier of that glom at the centroid location into a new 3D array of the same dimensions using the Python Imaging Library module. Similarly, the centroids of all glomeruli that contained pairs in the network excel export had white lines drawn between them within a separate TIFF, for connectivity visualization and validation. The resulting arrays were saved as a tiff and imported in Imaris overtop the original isosurface to visualize the relationships between the network and the 3D image. For ease of visualization, a smoothed-isosurface was applied overtop the network-connections overlay when showing large 5x networks in any figures.
A more flexible and user-friendly search algorithm was developed in tandem with the progression of the study to be used in place of the sub-array search used in V1 20x network generation for future analysis; one that allowed similar results without any downsampling step, providing a huge increase in computational performance and a more flexible algorithm. This method instead computed the distance transform of an inverted version of the gloms array to create a map of what regions were closest to which glom. A dilation could then be computed in binary, before the binary dilated regions were assigned a glom label based on their closest glom in the distance transform. In essence, this allowed multi-Gigabyte node arrays to rapidly dilate without losing their labeled IDs. This dilated space could then be searched for connections with nerves all at once, rather than having to use individual subarrays per glomerulus. To compensate for the lack of any downsampled steps, this alternate algorithm allowed the edges to instead be dilated by some factor, resulting in similar logic and outcomes as the V1 algorithm, but with enhanced user-parameter control. Networks generated by V1 and V2 for the neuroglomerular networks were overall very similar, but not necessarily identical due to V1 dilating nodes in a lower-fidelity downsampled array compared to V2, which permitted more pixel-specific dilation in the full-sized array, V1 utilizing downsampling to reconnect nerve masks compared to V2 using dilation in a full-sized array, again providing V2 greater specificity, and finally V1 utilizing dilation followed by watershedding to create 5x networks, leading to the loss of a small number of node identities, compared to V2.
The Python module^52^ was used to create a network from the pairwise Excel file containing all glom-glom connections. Community detection within the network was accomplished via the networkX Louvain algorithm^53^. Louvain was chosen as it is an efficient community detection algorithm and is well-suited for defining communities to use as clustering variables in analysis^54^. When generating Louvain networks, any repeated edges connecting the same pair of gloms gave that connection an additional weight in the network, with a weight of one meaning that only one neural connection between gloms was found, a weight of two representing two connections, etc. All community-related quantifications utilized this weighted Louvain detection, but in instances of large, dense 5x networks, an unweighted label propagation community detection algorithm^55^ was used solely for visualization due to its high speed.
Statistical analysis of networks was largely using the Python NetworkX package. Louvain networks were quantitatively described by the networkX modularity function^25^, obtained from Louvain community detection. Other statistics obtained directly from NetworkX methods included counts of gloms in a network, the number of connections per glom in the network, and the total number of communities detected by the Louvain algorithm. For quantitative analysis at 20x, only the largest connected component within an FOV was used to obtain a modularity score, as we were primarily interested in micro-network structure between gloms at high resolutions, meaning smaller communities at the edges of networks were excluded from this portion of analysis. High density FOVs (SK1 and SK2) were prioritized for network analysis since they provided by far the most information in a single image, but a small number of low-density FOVs were also included to demonstrate the gradient in community structure between networks containing a low and high number of gloms.
To establish significance in the community structures, all analyzed network FOVs had a random network generated of equivalent size to the network’s largest component, and containing the same number of nodes and edges. Nodes in the random network were iterated through serially and randomly assigned a partner until their network contained an equivalent number of edges as the volume from which they were derived The random network for each volume was generated 100 times, and its modularity was averaged. This averaged, random distribution of modularities was finally compared to the distribution obtained from samples using a two-sample t test. This method was also used to create random networks for 5x comparison. (Supplementary Data 9)
Creating a network out of a 5x sample with the nerve trunk present often revealed highly aggregated networks with low modularity and a non-obvious presence of communities. This was due to the majority of gloms being along a primary nerve tract that stretched through the entire volume, while the reduced resolution of 5x made assessing finer glom-nerve-glom connections unfeasible. This tract included the large nerve trunk that followed the arcuate vessel, and its upward projections alongside interlobular vessels into the cortex to seek out glomerulus vascular poles.
Since, we were interested in analyzing the nature of secondary glom-nerve-glom connections, it was pertinent to also develop an analysis method that disregarded this primary trunk to reveal communities in 5x. This was accomplished by removing the trunk after breaking the nerve network into discrete glom-nerve-glom units. After the trunk was removed, all connections that remained were either those that traveled directly from glom-to-glom, resulting in a massive increase in overall network modularity. Primary trunk removal was indicated by the modularity spiking by several factors after this algorithmic step, after which analysis on the network could proceed. In two instances, the algorithmic step to remove the largest nerve structure after discretization had to be repeated as the arcuate-associated nerve was outside the FOV, meaning there were multiple trunks arising from large interlobular vessels. In one instance (SK18), the main nerve trunk was not removed due to being outside the volume, causing the network to have an initially high modularity. It must be noted that due to the heterogeneity of the samples, the degree to which the trunk was removed likely varied, and therefore these networks are best treated as predictions at secondary structure rather than true descriptions.
5x networks were also analyzed with Louvain detection and other associated functions from networkX. To establish the successful removal of a trunk, modularity was first calculated as an aggregate for all disconnected components, rather than only on the largest connected component. This meant treating each component as its own community in the overall network, rather than analyzing each component separately. A high spike in modularity meant secondary groups had properly split from the main trunk. Further analysis was then performed on large connected components of the resulting 5x networks. The proportion of glomerular interconnectivity was likewise obtained by dividing the number of glomeruli in the trunk network by the number in the largest component of the trunkless network, acting as an approximate descriptor of what size modules emerged with secondary connectivity after trunk removal (Supplementary Data 13).
Degree distribution analysis of both 5x and 20x networks was accomplished in Python using NetworkX. First, the node-pair Excel files were used to populate a network with repeated connections assigned an additional weight. The degree of each node, representing the number of connections that node possessed, was obtained for each node from networkX and saved in Excel. The Excel FREQUENCY() function was used to obtain a count of the number of nodes possessing some degree k, which was converted into a proportion of nodes possessing some degree k from 1 to the maximum degree in the network. These values were used to generate a scatterplot with degree on the x axis and the proportion of nodes with that degree on y. As any degree corresponding to 0 nodes made power-curve modeling impossible, all proportions were incremented by 0.001 prior to modeling (Supplementary Data 10).
After removing the nerve trunk in applicable volumes, 5x networks often revealed a nerve network connecting the majority of gloms with a strong community structure. This interconnected network and its associated nerves were isolated with Python. NetworkX was used to generate a list of all glom IDs in the target connected network component, as identified by the networks script. The numpy function isin() was then used to create a binary mask from all interconnected gloms contained within this network, which was saved as a tiff and exported on top of the Imaris isosurface. Nerve connections between these gloms were isolated by converting the interconnected glom network mask to a boolean array and multiplying it against the labeled, dilated edge array. Edge IDs that remained were compiled in a list. The numpy function isin() could then be used to remove all nerves with IDs not present in this list, and the isolated nerve array was saved as a TIFF to be imported into Imaris.
Potential mother glomeruli were defined as nodes in the network that possessed connections between multiple communities as detected by the Louvain algorithm. These were detected from the network by first extracting the community partition from Python that grouped the nodes into communities. Each edge (containing a pair of nodes) that thereafter had two glomeruli in a separate communities/partitions represented a connection between communities and therefore a mother glomeruli.
By isolating subnetworks of gloms, we used the 5–20x coregistration to compare network prediction ability of gloms between 5x and 20x volumes.
Due to the lower resolution of 5x volumes, networks generated were expectedly less complex than the true nerve structure. For example, one hourglass structure rated by the networks script and easily visible at 20x in sample SK2 FOV 13 was seen to be broken into two disconnected communities in the corresponding SK2 5x network. As demonstrated by this observation, 5x networks lost connectivity both within and between communities, and therefore were best suited for finding high-throughput mother gloms and appraising a simplified, macro community structure.
Analysis of confocal microscopy datasets was done with the same strategy as the 3D automated 20x innervation, where a skeletonized nerve segmentation was quantified within a designated region surrounding a segmented glomerulus mask. A labkit model for nerves was trained on the confocal datasets. To eliminate signal-intensity-based bias, the diabetic and reference confocal datasets were first stacked together and used to simultaneously train one model. This model was used to segment all nerves, with the resulting mask being skeletonized and appraised for skeleton density around glomeruli in Python (Supplementary Data 14).
All dataset visualization was done with Imaris, while graphical visualization was done with matplotlib or Microsoft Excel, and figures and movies were assembled in adobe illustrator and adobe premiere pro, respectively. All statistical testing was done using the scipy stats module. Data were tested as shown, with the exception of radial angles, for which all angles were transposed away from radial arrangement by adding 6.28 to any angles below two radians during script execution. This was because of the circular nature of radial distributions, where an item of two radians in actuality is only two radians removed from one of 6.28, rather than 4. For the purposes of these data, moving all low-value radian datapoints to the top of a now-linear distribution allowed the entire dataset to be statistically tested without being incorrectly appraised due to the radial nature of the distribution.
FIJI (ImageJ) version 1.54j.
The FIJI MorphoLibJ plugin version 1.6.1
Imaris version 10.1
Python version 3.11, with the following packages and
certifi = =2024.7.4
charset-normalizer = =3.3.2
cmake = =3.30.1
colorama = =0.4.6
contourpy = =1.2.0
cycler = =0.12.1
et-xmlfile = =1.1.0
fonttools = =4.45.1
idna = =3.7
imagecodecs = =2024.6.1
imageio = =2.33.0
itk = =5.4.0
itk-core = =5.4.0
itk-filtering = =5.4.0
itk-io = =5.4.0
itk-numerics = =5.4.0
itk-registration = =5.4.0
itk-segmentation = =5.4.0
joblib = =1.4.2
kiwisolver = =1.4.5
lazy_loader = =0.3
matplotlib = =3.8.2
mpmath = =1.3.0
networkx = =3.2.1
numpy = =1.26.2
opencv-python = =4.8.1.78
openpyxl = =3.1.2
packaging = =23.2
pandas = =2.1.3
Pillow = =10.1.0
platformdirs = =4.2.2
pooch = =1.8.2
powerlaw = =1.5
pyparsing = =3.1.1
python-dateutil = =2.8.2
Python-Louvain = =0.16
pytz = =2023.3.post1
requests = =2.32.3
scikit-image = =0.22.0
scipy = =1.11.4
six = =1.16.0
tifffile = =2023.12.9
tqdm = =4.66.4
tzdata = =2023.3
urllib3 = =2.2.2
validators = =0.33.0
xmltodict = =0.13.0
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Supplementary Information Description of Additional Supplementary Files Supplementary Movie 1 Supplementary Movie 2 Supplementary Movie 3 Supplementary Movie 4 Supplementary Movie 5 Supplementary Movie 6 Supplementary Movie 7 Supplementary Movie 8 Supplementary Movie 9 Supplementary Data 1–18 Reporting Summary Transparent Peer Review file
Source Data