Imaging cellular activity across all organs reveals body-wide circuits

Cells across an organism must continuously coordinate with one another to sustain life and adapt to changing external environments and alterations within the body. Homeostasis is maintained through a complex network of dynamic interactions. Disease can arise from breakdowns in these intercellular feedback mechanisms. Although biomedicine has identified key interactions, such as stress responses mediated by neuroendocrine signalling or neural pathways between the brain and the gut, a vast number of mechanisms of whole-organism function remain unknown³. Our understanding is limited by the challenges of accessing body-wide cellular dynamics, driving the need for technologies to measure, analyse and model body-wide control mechanisms at the cellular level.

Modern synergies between imaging technology and protein engineering allow for time-varying molecular signals to be recorded in tissues. This has caused revolutions in fields such as neuroscience through the imaging of calcium—a fast, universal intracellular messenger involved in a wide range of cellular processes, including neuronal action potentials¹—and other signals including voltage and neuromodulators across many neurons simultaneously⁴. However, time-varying activity patterns of most cell types in the body have not yet been recorded.

For most vertebrate models, optical access to large, opaque tissues poses a currently insurmountable challenge to whole-body imaging. Transparent vertebrate animals such as young zebrafish and adult Danionella cerebrum⁵ overcome this barrier, making them uniquely suited as models for in vivo studies of cellular dynamics across the body, offering unparalleled access to the inner workings of evolutionarily conserved organs such as the liver, pancreas, gut, brain, and the immune and cardiovascular systems.

This study introduces WHOLISTIC, a platform for in vivo imaging of cellular calcium dynamics, generalizable to other molecular dynamics, across nearly all cells of transparent vertebrates, such as the young zebrafish. By extending and integrating pancellular transgenic lines expressing genetically encoded calcium indicators in almost all cells in the body, high-speed volumetric fluorescence imaging, and a suite of computational methods for registration and cell population analysis, along with whole-body expansion microscopy (WB-ExM), we captured, analysed and interpreted cellular activity across tissues and organ systems.

To image dynamic fluctuations of intracellular calcium across cell types of the body of young zebrafish (Fig. 1a), we engineered a pancellular transgenic zebrafish line that expresses the genetically encoded calcium indicator GCaMP7f under the control of the ubiquitin promoter^6,7, Tg(ubi:tTA; TRE:GCaMP7f) (Fig. 1b,c). To avoid transient or sparse expression⁸, we used the binary expression system tTA–TRE to enhance the concentration of GCaMP7f⁹, which resulted in extensive and sustained pancellular expression with no observed pathological manifestations (Fig. 1b and Extended Data Fig. 1a,b). Here we chose to focus on intracellular calcium, given its broad involvement in numerous cellular processes across all cell types¹. The approach generalizes to other sensors of voltage, membrane tension and other molecules.

The utility of pancellular imaging is contingent on the ability to discern organs and cell types within densely labelled tissue. The ubiquitous expression of GCaMP7f in tightly packed cells might have led to a conglomeration of indistinguishable cells; however, we found that major organs and tissue types could be distinguished by their distinct visual textures arising from variations in cell size, morphology, subcellular structure, baseline calcium and sensor expression level (Fig. 1d). Together, this enables the recognition of bodily tissues, including muscle, liver, kidney, brain and intestine, as well as smaller cellular populations such as the pineal gland and the endocrine pancreas (Fig. 1d), subsequently verified at a more detailed scale using WB-ExM (Methods). Furthermore, we developed a pancellular Danionella cerebrum GCaMP transgenic line to adapt WHOLISTIC to adult transparent vertebrates (Extended Data Fig. 1c,d). Thus, the expression of calcium sensors under the ubiquitin promoter has the potential to enable the monitoring of dynamic cellular signals throughout the organism, both developing and mature, while also affording the identification of tissue types, with or without the addition of cell-type-specific markers.

To acquire and process spatiotemporal WHOLISTIC data and extract interpretable cellular calcium activity time series that can be related to the tissues of origin, we developed a data collection and analysis workflow (Methods and Fig. 1e–j). Acquiring high-quality data across the body was more challenging than imaging the brain due to greater tissue inhomogeneities and complex geometry that increase scattering and diffraction, such that light-sheet microscopy^10,11 was unsuitable. Among the microscopy techniques evaluated, spinning-disk confocal microscopy proved most robust to tissue scattering and diffraction (Extended Data Fig. 1e). As most major internal organs in zebrafish are located in the anterior third of the body, we opted to primarily image this anterior region from the side (Fig. 1e and Supplementary Videos 1 and 2).

Smooth and skeletal muscle contractions throughout the body result in the displacement of cells (Supplementary Video 3), such that accurately extracting the activity of individual cells over time requires accounting for this cellular displacement. Although motion in the brain is relatively rigid, motion in the viscera is much less constrained, making standard registration algorithms ineffective or slow. To overcome this challenge, we first performed dual-colour imaging of calcium dynamics concurrently with an anatomical reference channel to aid alignment (Extended Data Fig. 1f,g). Next, we developed a registration algorithm based on iterative optical flow estimation¹² that estimates motion patch-wise (11 × 11 pixels) and combines smoothness regularization, signal-to-noise weighting of gradients and multi-resolution registration to favour convergence to an optimal tissue alignment for each time point (Methods and Fig. 1f). This achieved high-quality alignment, and most cells were accurately registered throughout the experiments (Fig. 1g, Extended Data Fig. 2a–c and Supplementary Video 4). Subsequent analysis allowed for the identification of single-cell responses to externally applied stimuli, including previously unknown responses of cartilage chondrocytes to cold temperatures (Extended Data Fig. 3 and Supplementary Video 5) and responses to ketamine by cells located at the meningeal boundary of the brain (Extended Data Fig. 4). Thus, WHOLISTIC enables the analysis of spontaneous and stimulus-driven responses across the animal at the single-cell level.

Identifying functional tissue ensembles

To extract cellular-scale calcium activity traces from registered imaging data (Fig. 1h), we used a segmentation method based on constrained non-negative matrix factorization (Voluseg)¹³ (Methods and Fig. 1i). We tuned this approach conservatively to minimize false mergers, over-parcellating the data into subcellular segments, that can later be aggregated.

We found that the majority of cells exhibit measurable, time-varying calcium dynamics even at rest (Fig. 1j–m and Supplementary Video 1). Although consistent with the known importance of calcium signalling, to our knowledge, it has never been demonstrated that most cell types across disparate organs exhibit measurable second-scale calcium fluctuations even in the absence of specific stimuli, nor that the dynamic range of GCaMP7f is sufficiently broad to capture them.

The utility of pancellular imaging increases with the ability to readily assign cells to their organs and cell types of origin. Although imaged volumes can be manually annotated with the identity of tissue types based solely on anatomical features (Fig. 1d), it is laborious, especially for sparse or highly distributed populations, and misses finer delineations within cellular populations that are defined not by their appearance but by their calcium activity patterns. We thus explored whether a more data-driven methodology, centred on the aggregation of cells with related activity patterns, could more effectively delineate functional anatomical boundaries within and across tissues. This would also provide a more compact, dimensionally reduced and interpretable representation of the data. To this end, we implemented spectral clustering¹⁴ and used a coherence-based similarity metric. Unlike correlation, coherence is invariant to phase lags and thus groups cells engaged in temporally structured activity such as travelling waves (Methods and Fig. 1j).

Examining the anatomical footprint of the spectral clusters, we found that they predominantly corresponded to sub-regions of individual tissue types, such as the kidney, gallbladder, liver, gut and muscle (Fig. 1j,k,m and Extended Data Figs. 5 and 6). This can be leveraged to organize the data into interpretable and readily identifiable clusters at the tissue level (Fig. 1j,k), which we denote ‘functional tissue ensembles’. These can then be annotated with tentative tissue-type identities when discernible. As an example, the WHOLISTIC workflow parcellates the pronephric kidney nephron into contiguous functional tissue ensembles, each exhibiting distinct dynamics, that together couple into travelling waves that occur in bursts (Fig. 1j,k). The spatial propagation and inter-burst quiescence suggest that kidney filtration and/or intraluminal flow—conceptualized as a continuous or continuously oscillating process^15,16—may in fact be intermittent and pulsatile on shorter timescales, warranting follow-up experiments. Clusters across the body exhibited a range of calcium activity, showcasing diverse temporal spectra, including oscillatory, pulsatile and bistable activity (Fig. 1l,m).

Collectively, these findings demonstrate how the ubiquitous expression of molecular sensors, in conjunction with advanced motion-corrected spinning-disk confocal microscopy and the computational workflow for cellular and tissue segmentation, together comprising WHOLISTIC, facilitate the study of cellular interactions across tissues of the body.

Body-wide skeletal muscle coordination

To exemplify the use of WHOLISTIC to study coupling of calcium activity within tissue types, we began by inspecting functional tissue ensembles corresponding to skeletal muscle, in which the rise of intracellular calcium induces contraction. Muscle ensembles were classified based on their distinctive cellular morphology (Fig. 1d), stereotyped spatial location (Fig. 2a) and characteristic activity waveform. Next, the ensembles were organized through hierarchical clustering of their calcium activity patterns and visualized by projecting their spatial footprint back into anatomical space (Methods and Fig. 2b). This approach readily enabled the identification of the primary epaxial and hypaxial trunk muscles, which exhibited correlated activity patterns (as represented by the union of the blue and green clusters in Fig. 2b), underscoring their propensity to contract synchronously, as documented in other species¹⁷. This method also identified the abdominal (also known as hypaxialis), sternohyoid, pectoral fin and branchial muscles¹⁸ (Extended Data Fig. 7a). More surprisingly, the anterodorsal portion of the epaxial muscle grouped in a distinct functional tissue ensemble (‘cervical epaxial muscle’), with calcium activity correlated with the abdominal muscle (Fig. 2c, right), a finding that, to our knowledge, has not been previously reported. These results can be visualized within the raw data by examining time points at which specific muscle groups are selectively activated (Fig. 2c, left).

To ascertain the neural underpinning of the independence of activity of the cervical epaxial muscle from the primary epaxial muscle, we traced the anterior motor nerve tracts using a transgenic line labelling motor neurons, Tg(VAChTa:eGFP) (Methods). This revealed that, although the spinal nerves predominantly innervate the majority of the epaxial and hypaxial musculature, the spino-occipital nerve—known to innervate the ventral sternohyoid and pectoral fin muscles¹⁹—extends axons dorsally to innervate the cervical portion of the epaxial muscle (Fig. 2d and Extended Data Fig. 7b). This dual innervation supports the independent regional activation, consistent with previous studies in other species demonstrating regional activation²⁰. Together, this shows how WHOLISTIC can uncover, even within extensively studied tissue dynamics—muscles—new functional units, as well as reveal previously unidentified synergies and independencies among muscle groups (Fig. 2e).

Muscle-coupled dynamics

To demonstrate the ability of WHOLISTIC to identify inter-organ interactions, we next studied the coupling between skeletal muscle and other bodily tissues. Beyond moving the animal, skeletal muscles generate internal mechanical forces detected by mechanosensitive cells body wide²¹. Connections through the nervous, vascular and other systems can additionally induce coupling between muscles and other organs. Yet, the acute effects of muscle activity on non-muscle tissue types and vice versa remain largely uncharted.

We thus asked how motor activity acutely affects organ and tissue dynamics throughout the body. We examined the temporal relationships between the activation of skeletal muscles and calcium signalling in other cells of the body (Fig. 2f–h) by fitting a lag-regression model to all cells within the imaged volume, in which cellular calcium activity is modelled as the convolution of muscle activity and a cell-specific motor-response kernel (Methods, Fig. 2f and Extended Data Fig. 7c,d; we note that this model does not purport to infer causal direction). We quantified motor-related cellular activity throughout the imaged volume using statistical parametric mapping, which highlighted responses in various tissues, including the brain, skin epithelium, notochord sheath cells and gills (Fig. 2g,h). This coupling exhibited distinct, tissue-specific responses, often enduring beyond the timescale of muscle activity, except in the brain, where a cluster of brainstem neurons showed short coupling timescales (Fig. 2g,h, row 4).

Within the brain, beyond motor-coupled neurons, a cohort of cells in the midline of the hindbrain and spinal cord exhibited temporally extended responses to muscular activity, enduring significantly longer than the neuronal responses (approximately 15 s; Fig. 2g,h, row 2). On the basis of their anatomical position and morphology, we hypothesized that these were ependymal cells, glia-like cells that line the ventricles and the central canal²². To test this hypothesis, we generated a transgenic line expressing GCaMP7f specifically in ependymal cells (Tg(foxj1a:GCaMP7f)), as well as a new pancellular transgenic line expressing a red fluorescent calcium sensor²³, Tg(ubi:tTA;TRE:jRGECO1b). Crossing these two lines enabled the simultaneous recording of definitive ependymal cell calcium activity, as well as comprehensive whole-body dynamics in independent channels (Fig. 2i) and confirmed that ependymal cells exhibit Ca²⁺ transients subsequent to motor activity, with response amplitudes that are proportional to the intensity of motor events (Fig. 2j,k and Extended Data Fig. 7e). This suggests that ependymal cells integrate body movements into their cellular state, potentially contributing to ventricular and cerebrospinal fluid homeostasis.

Revealing this coupling between ependymal calcium and muscle activity illustrates the sensitivity of WHOLISTIC for discovering functional properties even amid sparse and distributed cell populations.

Optogenetic test of brain–body coupling

To expand the analysis to include visceral organs as regressors, we considered relationships between cellular and smooth muscle activity, which exist in tissues forming the gut, the lining of blood vessels, the gallbladder and others, by incorporating smooth muscle ensembles into the lag-regression model. Among the various relationships identified, the model captured a relationship between the caudal hindbrain and smooth muscle of the gastropharyngeal sphincter. Indeed, this regressor accounted for the largest fraction of variance of a small group of cells within the caudal hindbrain, specifically within the region of the motor vagal nuclei (Fig. 2l and Extended Data Fig. 7f). As the model is agnostic to causal direction, this is consistent with these motor vagal neurons driving the sphincter calcium activation via the vagus nerve. Given the established connection between the motor vagus and the gut across species²⁴, this finding represents the identification of a conserved vagus–gut circuit from an unbiased WHOLISTIC analysis.

We sought to probe whether the relationship between this motor vagal area and the gastropharyngeal sphincter was causal. To do this, we combined WHOLISTIC with optogenetics by crossing a transgenic line expressing channelrhodopsin in the motor vagus with the pancellular genetically encoded calcium indicator line (Tg(VAChTa:gal4;UAS:CoChR-eGFP); Tg(ubi:tTA;TRE:jRGECO1b); Extended Data Fig. 7g). We found that activation of neurons within the most caudal segment of the motor vagus induces a dose-dependent calcium activation across the gastric sphincter (Fig. 2m–o), recapitulating the observed coupling identified by the model. Moreover, examination of the entire motor vagus through the successive activation of subsets of neighbouring vagal neurons revealed a topological mapping between the motor vagus and visceral tissues along the anteroposterior axis (Fig. 2p), with more anterior motor vagal neurons controlling more anterior pharyngeal muscles²⁵. We noted that neural activity exerted control of smooth muscle up to the anterior gut but not beyond. To assess whether there was an anatomical underpinning for this spatially restricted control, we visualized the distribution of motor vagal axons along the gastrointestinal tract using WB-ExM (Methods and Extended Data Fig. 7h,i), revealing a significant drop in innervation density at the level of the anterior gut, paralleling the functional results and offering a potential anatomical basis for the functional findings. This combination of WHOLISTIC, optogenetics and ExM-based anatomical mapping underscores its capacity to rapidly uncover and mechanistically examine brain–body circuits.

As a demonstration of the broad applicability of WHOLISTIC, we next considered the converse of movement-coupled modulation of tissue activity: processes that occur during extended periods of motor quiescence. In between phases of sustained motor activity, animals engage in equally important periods of muscle inactivity associated with rest, sleep and tonic immobility. These quiescent states are active phases of physiological regulation, yet although the neural substrate during them is extensively studied, the contributions of other cell types remain less characterized.

We therefore aimed to identify cells—whether neuronal or non-neuronal, within or outside the brain—that exhibited heightened activity during motor-quiescent periods. We analysed the variance in cellular calcium activity during episodes of extended muscle inactivity (more than 10 min) compared with periods of regular motor activity. This analysis highlighted a region in the central nervous system that exhibited increased activity during rest (Extended Data Fig. 8a). These signals were identifiable along the midline of the hindbrain and spinal cord, displaying high-amplitude, ultraslow oscillatory activity during motor-quiescent phases, with a periodicity of 3–7 min (Fig. 3a and Extended Data Fig. 8b,c). By contrast, during periods of motor activity, faster and abrupt muscle-locked activity was observed in these cells, and the slow oscillations were either absent or significantly attenuated (Fig. 3b–d). Imaging neuronal and astrocyte-specific lines (Tg(elavl3:H2B-GCaMP7f) and Tg(gfap:jRGECO1b)) failed to recover such dynamics, leading to the hypothesis that these slow signals may originate from ependymal cells. Ependymal cells have been associated with sleep and alterations in cerebrospinal fluid flow^26,27 and are juxtaventricular cells, with their cell bodies residing along the central canal and ventricles²². However, the oscillations were primarily observed at the base of the hindbrain (Fig. 3e), above the notochord—an area far from the ventricles or central canal, where ependymal cells are typically assumed to reside (Extended Data Fig. 8d)—bringing into question whether the oscillations originate from ependymal cells. To resolve this discrepancy, we sought to visualize the detailed morphology of ependymal cells that reside deep within the brain.

Whole-body expansion microscopy

Despite the relative transparency of the zebrafish, the accumulation of refractive index variations across the body results in light scattering, making it difficult to visualize the fine anatomy of such midline structures. Existing clearing and expansion methods are either only applicable to younger samples (5 days or less post-fertilization), achieve clearing at the expense of extensive proteolytic digestion^28,29 or require extensive wash steps to clear adequately without expansion³⁰. To overcome these limitations, we developed an advanced enzyme-free, rapid and robust WB-ExM protocol² (Fig. 3f) that enables high-quality clearing and up to 5× uniform expansion, along with the acquisition of molecular information at subcellular resolution (Extended Data Fig. 9a–i) throughout entire larval and juvenile zebrafish (Extended Data Fig. 9j), as well as mature Danionella cerebrum (Extended Data Fig. 9k). This protocol combines high-temperature disruption with a multi-round embedding strategy, gradually reinforcing and expanding the sample to ensure homogeneous expansion (Methods). WB-ExM retains immunofluorescence signals (WB-ExM-IF; Methods, Extended Data Fig. 9a,b and Supplementary Video 6) and, with modifications, is compatible with in situ hybridization (WB-ExM-FISH; Methods; Extended Data Fig. 9e–g). In addition, the use of total protein stains yields rich anatomical information (WB-ExM-Histo) with which to contextualize the immunofluorescence and FISH signals³¹ (Extended Data Fig. 9h and Supplementary Video 7). To quantitatively assess expansion quality across the sample, we developed PhotoMap, an unbiased method for extracting the deformation field of the gel (Extended Data Fig. 10), inspired by GelMap³². Expansion was found to be largely uniform (Extended Data Fig. 10f,g), except in regions of mineralizing bone, where tissue pinching was observed but remained spatially confined (Extended Data Fig. 10c,d). These tools provide a high-throughput route to comprehensive whole-body anatomical and molecular data, complementing histological, X-ray and electron microscopy approaches, and a vital basis for interpreting and validating WHOLISTIC findings.

We used WB-ExM-IF to characterize the fine morphological features of ependymal cells (Fig. 3g and Extended Data Fig. 8e–g) and found that this population of hindbrain ependymal cells, typically considered to be small cuboidal cells³³, actually extends long and dense ventral projections (approximately 120 μm), forming a sheet-like structure that runs along the entire midline of the hindbrain, reaching the floor of the brain (Fig. 3g,h and Supplementary Video 8), thereby explaining the spatial distribution of the oscillatory signal. The morphology of these ependymal cells echoes a class of ependymal cells, including tanycytes, that retain morphological features of their progenitor, radial glia cells²², and are widely present in the adult mammalian brain, suggesting an evolutionarily conserved lineage³⁴. These cells also displayed previously undescribed morphological features (Extended Data Fig. 8f–k): processes ensheathing commissural fibre crossings approximately 25 μm below the central canal, ensheathment of passing arteries and evenly spaced lateral projections surrounding neuronal tracts exiting the spinal cord echoing those recently identified in mice^26,35. Such findings, not readily observable in non-expanded samples, show how WB-ExM resolves the cellular and extracellular provenance of functional signals, supporting hypothesis generation.

Cell-type-specific imaging of ependymal cells recovered the motor quiescent-locked ultraslow oscillatory dynamics, confirming the ependymal cell origin of the ultraslow oscillations (Fig. 3i–l). To understand the spatiotemporal structure of ependymal activity, we performed cross-coherence analysis to assess the presence of coupling between the oscillatory dynamics of cells (Methods). This uncovered a spatial organization manifesting as local waves along the anteroposterior axis of the fish, encompassing both the hindbrain and the spinal cord (Fig. 3j), forming a block-like structure in the cross-coherence matrix with a length scale of approximately 75 μm. To quantitatively evaluate this, we derived a method for utilizing coherence as a metric in k-means clustering (Methods); the algorithm, blind to the anatomical location of the cells, recovered spatially continuous clusters (Fig. 3k). Ependymal cells are interconnected via gap junctions³⁶, probably accounting for the localized propagation of the signal. Calcium within ependymal cells has been associated with F-actin-mediated cellular shrinkage, resulting in increased paracellular space between ependymal cells, thereby facilitating efflux through perineuronal routes²⁶. These findings support the hypothesis that these cells form an interconnected network across the brain and spinal cord, capable of integrating and modulating neural activity, vascular signals, and cerebrospinal fluid composition. Moreover, it suggests that the observed state-dependent calcium oscillatory modes of ependymal cell activity may correspond to prolonged periods of altered cerebrospinal fluid exchange during extended motor quiescence.

Together, these analyses show that WHOLISTIC can identify state-dependent sparse cellular signals within densely labelled tissue that can subsequently be validated and followed up using cell-type-specific transgenics and WB-ExM, and highlights previously overlooked populations active during specific behavioural states.

To demonstrate the power of WHOLISTIC for studying organism-wide responses to physiological stress, in addition to spontaneous activity (Figs. 2 and 3), we examined how animals adapt to hypoxia: a condition of reduced oxygen availability that engages systemic responses across tissues and organs. Hypoxia is a fundamental physiological challenge that arises in contexts ranging from high-altitude adaptation to stroke and heart failure, and triggers genetic, metabolic, cellular and behavioural adaptations that maintain oxygen homeostasis and prioritize energy allocation to critical tissues³⁷. However, the dynamic interplay between single-cell activity changes and whole-organism physiology during hypoxia remains poorly understood, largely due to the technical challenges of simultaneously monitoring multiple tissues in real time.

To systematically evaluate the systemic effects of hypoxia across the body, animals were exposed to alternating periods of normoxia (normal O₂ levels, 21%, approximately 15 min) and hypoxia (reduced O₂ levels, 10% O₂, approximately 15 min; Fig. 4a). We observed widespread changes in baseline levels of calcium, quantified by computing the difference in average baseline calcium levels between normoxic and hypoxic phases, referred to as the ‘oxygen modulation score’ (Methods and Fig. 4b). Of note, gastrointestinal tissues showed a systematic increase in baseline calcium levels during hypoxia (Fig. 4c), potentially reflecting reduced calcium buffering capacity of the cells as a result of diminished oxygen availability³⁸. Other organs also showed changes in baseline calcium levels, although the gastrointestinal tract, at the average tissue level, showed among the most notable initial increase in baseline calcium levels (Extended Data Fig. 11a). Given that the arterial blood supply to the brain and gut would contain similar oxygen levels, we sought to understand these disparities in responses.

While considering possible causes for the differential modulation of calcium in the visceral organs compared with the brain, we noticed that the mesenteric artery, the primary artery supplying the visceral organs, appeared constricted during episodes of hypoxia (in unregistered imaging data where tissue motion is visible). This suggests a reallocation of oxygen across the body by rerouting blood flow during physiological stress, a conserved response³⁹ that has never been imaged in real time. To assess this arterial constriction in more detail, we imaged a transgenic line that selectively labels blood vessels, Tg(flk1:dsRED-CAAX) (Fig. 4d) and confirmed the constriction of the main mesenteric artery during hypoxia (Methods, Fig. 4e,f). As the diameter of a blood vessel is an indirect indicator of blood flow, we imaged a transgenic line labelling red blood cells, Tg(gata1:dsRED) (Extended Data Fig. 11b), and found a near-complete cessation of blood flow to the gut during hypoxic conditions, which recovers upon returning to normoxic conditions (Extended Data Fig. 11c and Supplementary Video 9). By contrast, blood flow to the brain and muscle remained largely unaltered, potentially explaining the difference in baseline modulation between these tissues. Thus, consistent with mammals and certain aquatic animals^39,40, WHOLISTIC shows that hypoxia induces a rapid redistribution of blood away from the gut, and reveals the full temporal dynamics of the surprisingly fast redirection of oxygen as well as its spatial distribution, and concurrent activity changes in other cells across the body.

To test whether the reduction in gastric blood flow was neurally regulated, we inhibited the nervous system using the anaesthetic tricaine and observed that visceral blood flow was no longer diminished during hypoxia (Extended Data Fig. 11c), suggesting that the reduction in visceral blood flow is mediated by neurons, and showing that this brain–body circuit is already developed at 7 days post-fertilization. To more directly test this, we optogenetically inhibited the hindbrain during hypoxia while measuring blood flow (Fig. 4g and Extended Data Fig. 11d). We found that inhibiting the hindbrain induced, within a few seconds, the return of blood flow to the gut (Fig. 4h,i). In spite of the hypoxic condition, near-normal blood flow was maintained for as long as the hindbrain was inhibited, with no corresponding effect observed when directing the stimulation light to a control location (Fig. 4j). Thus, hindbrain activity is necessary for the constriction of the mesenteric artery during hypoxia, which we speculate serves to redistribute oxygen, thereby preserving energy for the brain and other tissues whose activity must be prioritized, and may also be a control mechanism to decrease enteric cellular activity.

During hypoxia, we also observed a set of cells in the medial plane, below the notochord, that showed a systematic and sustained increase in activity during hypoxic periods (Fig. 4k,l). Given their location, and that sympathetic neurons induce vasoconstriction via noradrenergic activation of vascular smooth muscle cells⁴¹, we hypothesized that this cluster represented the sympathetic ganglion. Using cell-type-specific imaging (Tg(th:gal4; UAS:GCaMP6f))⁴², along with WB-ExM-IF and WB-ExM-FISH, we showed that this group of cells indeed corresponds to sympathetic neurons and that they send axons along the mesenteric artery (Fig. 4m–o and Extended Data Fig. 11e). Together, these results establish a brain-to-mesenteric artery pathway already functional at 7 days post-fertilization. Combined with optogenetic perturbation, WHOLISTIC can thus causally dissect cell-type-specific brain–body pathways engaged during physiological stress. This work, moreover, establishes the young zebrafish as a powerful model for studying the consequences of stress-induced blood shunting and, more generally, organism-wide feedback loops engaged by stress.

The organism is a cohesive entity, comprising networks of cells that communicate and coordinate across multiple scales⁴³. Feedback loops within these networks enable animals to respond to, predict and prepare for both internal and external challenges. Although many regulatory pathways—such as those governing hunger, glucose balance and immune responses—are well characterized⁴⁴, the complete picture of organism-wide control systems remains a formidable challenge. Given the interdependence of cellular interactions, understanding the organism as an integrated whole is essential^3,45,46,47. Yet, methods that probe the entire vertebrate organism while retaining cellular-level access have remained elusive. This study introduces a novel approach for imaging and analysing genetically encoded sensors expressed ubiquitously across cells, providing unprecedented insights into cellular interactions and organism-wide control systems.

Future implementations of WHOLISTIC in freely behaving animals will enable the study of body-wide cellular control systems in naturalistic settings, rather than an embedded preparation. This will facilitate investigation of the interplay between physiological states and behavioural strategies, such as goal-driven navigation and sleep–wake transitions. A parallel route will involve experimental systems incorporating multimodal virtual-reality environments⁴⁸ composed of visual, thermal, chemical, mechanical and other stimuli that respond to the behaviour of an animal. Further advances in protein engineering and microscopy that improve closed-loop monitoring and activation of cellular signals across the body⁴ will open new avenues of scientific inquiry.

Calcium sensors^6,23 provide an effective proxy for cellular activity due to the ubiquitous role of calcium in cellular processes¹. The finding that most cells exhibit calcium-activity dynamics within the range of standard sensors has enabled extraction of many cellular dynamics. However, this remains a limited view of signalling within and between cells, which occurs via a multitude of molecular, voltage and mechanical pathways. Measuring additional molecules and physical cellular properties will provide deeper insights. WHOLISTIC can be extended to co-image signals within the rapidly expanding repertoire of molecular, voltage and mechanical sensors^{49,50,51,52,53,54}, providing access to additional information channels, including metabolic states, hormones, ATP and neuromodulators that mediate intercellular communication.

Given the scale and complexity of multicellular interactions, uncovering higher-order interactions will require advanced analytical and modelling frameworks leveraging modern artificial intelligence, including graph neural networks⁵⁵, large-scale biophysical models and model-free predictive methods such as empirical dynamic modelling⁵⁶.

Platform automation⁵⁷ and artificial intelligence-in-the-loop systems, combining measurements of whole-body cellular dynamics with online cell-targeted perturbation, will support gaining fundamental biological insights as well as prototyping real-time medical interventions^57,58,59.

The fact that all findings presented here—pertaining to different cell types, in different parts of the body, under different physiological assays and screens—arose from the same WHOLISTIC methodology underscores its broad applicability to fundamental biology and medical research. The methodology can serve as a discovery tool, whereas validation and deeper insight can be gained using more specific molecular techniques, as performed here to understand phenomena such as the brainstem’s control of blood flow redistribution during physiological stress and the molecular identification of the quiescence-related signals as originating from ependymal cells.

Cellular function is determined by molecular and biophysical properties, tissue environment and connectivity to other cells^60,61. Combining organism-wide cellular activity imaging with whole-body expansion microscopy offers a means to bridge function and mechanism at scale, but will require the development of body-wide registration methods, capable of matching cells between in vivo and ex vivo expanded data cell by cell.

Much information throughout the body flows through the central and peripheral nervous systems, whose connectivity is becoming increasingly accessible through light-based connectomics⁶² and machine learning algorithms for connectome reconstruction⁶³. This work will seed an atlas of the mechanistic and molecular basis of functional coupling at a cellular level across the organism.

Beyond fundamental biology, WHOLISTIC enables disease modelling that tracks the effects of pathologies and treatments across spatial scales, from cells to the whole body, and temporal scales, from seconds to days. Most drug-screening workflows examine selected aspects of drug action, such as molecular binding⁶⁴, effects on individual cells or tissue subsets⁶⁵, or behaviour^66,67. WHOLISTIC offers a complementary lens through which to observe the effect of drugs and genetic interventions on the entire body, revealing off-target effects and system-wide interactions that would otherwise be missed when studying cells, tissues or organs in isolation.

In conclusion, WHOLISTIC bridges cellular and organismal physiology, systems neuroscience and behaviour, opening a new frontier in systems biology. By embracing organismal complexity while retaining the individual cell as a fundamental unit of analysis, this work provides a holistic yet mechanistic approach to unravelling the cellular interactions governing health and disease across organism-wide networks.

Experimental model and subject details

Zebrafish husbandry

Zebrafish were reared at 28.5 °C in 14–10-h light–dark cycles (conductivity of 1,000 μS, adjusted via Instant Ocean Sea Salt (approximately 30 g l⁻¹), pH 7.0, adjusted using sodium bicarbonate)⁶⁹. Zebrafish from 5 to 14 days post-fertilization were fed rotifers and used for experiments. All experiments complied with protocols approved by the Institutional Animal Care and Use Committee of Janelia Research Campus. Zebrafish sex cannot be determined until approximately 4 weeks post-fertilization⁷⁰, so the sex of the experimental animals was unknown. Where relevant, fish were randomized across conditions.

No blinding was used in either data collection or analysis. Blinding during data collection was not possible because the experimental condition determined the acquisition protocol and was therefore necessarily known to the experimenter at the microscope. Blinding during analysis was not applied because all reported quantities were extracted by automated pipelines using identical parameters across conditions.

Danionella cerebrum husbandry

Danionella cerebrum were reared at 26.5 °C in 14–10-h light–dark cycles (conductivity of 450 μS, pH 7.5). Feeding protocols were adjusted according to age: (1) at 5–28 days post-fertilization, rotifers were administered once daily, (2) at 16–28 days post-fertilization, in addition to rotifers, GEMMA 75 was provided twice daily, and (3) at 29 days post-fertilization and beyond, the diet was composed of GEMMA 75 twice daily and Artemia once daily. Adult Danionella cerebrum were maintained in group housing with stock density of approximately 45 fish in 3.5-l tanks (Tecniplast). Fish younger than 6 weeks of age are sexually immature and could not be sexed; the sex of fish older than 6 weeks of age is mentioned in the main text. For egg collection, 10-cm-long custom-made acrylic tubes were used. All experiments complied with protocols approved by the Institutional Animal Care and Use Committee of Janelia Research Campus.

Zebrafish transgenics and transgenesis

Transgenic zebrafish were maintained in the Casper or Nacre background⁷¹. All lines were generated using the Tol2 system⁷² and genes were codon optimized using CodonZ⁷³. Codon-optimized GCaMP7f and jRGECO1b were synthesized (Twist) and used for subsequent cloning. For cloning, restriction digest cloning was used throughout and all genes (GCaMP7f, jRGECO1b and tTA) were cloned with a preceding Kozak sequence and followed by an SV40 poly(A) signal sequence. Plasmids (150 ng μl⁻¹), along with Tol2 transposase mRNA (50 ng μl⁻¹), were co-injected (0.5 nl total injected volume) into one-cell stage embryos. Embryos were screened at 7 days post-fertilization for expression, and positive embryos were reared to maturity. At maturity, these adults were individually screened for dense expression in progeny, and the best founders were retained. We note that due to the non-deterministic landing site of the transgene, founders have variation in expression and need to be carefully screened for dense expression (Extended Data Fig. 12).

The ubi:tTA and TRE elements were obtained from multiple plasmids, a gift from D. Feliciano and I. Espinosa-Medina. These included a ubiquitous promoter containing vector or p5E-ubi⁷ (Addgene 27320), a vector containing the tTA advanced Tet-off transcriptional activator from pTet-Off Advanced Vector (631070, Takara) inserted into the multiple-cloning site of the pME entry vector⁷⁴, and a vector containing the tetracycline-responsive element promoter p5E-TRE⁷⁵. To generate the Tg(ubi:tTA;TRE:jRGECO1b) animals, the ubi:tTA;TRE elements were cloned and a codon-optimized jRGECO1b sequence placed downstream. To generate the Tg(ubi:tTA);Tg(TRE:GCaMP7f) animals, plasmids containing ubi:tTA and TRE:GCaMP7f were independently cloned and co-injected at equimolarity. To generate the Tg(foxj1a:GCaMP7f) animals, the foxj1a promoter was cloned (Addgene plasmid 163829)³⁶ and a codon-optimized GCaMP7f sequence placed downstream using restriction digest cloning. To generate the Tg(elavl3:gtACR2-eYFP) transgenics, the promoter was cloned using a known elavl3 promoter sequence⁷⁶ and gtACR2-eYFP sequence placed downstream⁷⁷. To generate the Tg(β-actin2:mCherry-CAAX; myl7:GFP) transgenic line the bactin2 promoter was cloned (Addgene plasmid 82583)⁷⁸ and the mCherry-CAAX sequence placed downstream. For optogenetic activation of motor vagal neurons, the transgenic lines Tg(VAChTa:Gal4)⁷⁹ and Tg(UAS:CoChR-eGFP)^jf44 (ref. ¹³) were utilized. To record activity of the sympathetic ganglia, the transgenic lines Tg(th:Gal4)⁴² and Tg(UAS:GCaMP6f)^jf46 (ref. ¹³) were crossed and imaged. To quantify blood vessel diameter and track blood flow, the transgenic lines Tg(flk1:dsRED-CAAX)⁸⁰ and Tg(gata1:dsRED)⁸¹ were imaged. To check the colocalization of neurons and astrocytes with the identified ‘brain-border’ population, and its location relative to the brain’s basement membrane, the transgenic lines Tg(elavl3:H2B-jRGECO1b)⁸², Tg(gfap:jRGECO1b)¹³ and TgBAC(lamC1:lamC1-sfGFP)⁸³ were respectively used. Additional lines used for expansion microscopy were: Tg(isl1CREST-hsp70l:mRFP)⁸⁴, Tg(phox2bb:eGFP)⁸⁵ and Tg(foxj1a:eGFP)⁸⁶.

Danionella cerebrum transgenics and transgenesis

To generate pigmentless D. cerebrum mutants, we used the CRISPR–Cas9 genome-editing technique to disrupt the function of the mitfa gene, following established protocols⁵ and utilizing a mitfa-targeting guide RNA (sequence: CAGCATTATACACTAAGAGT). Mutations were confirmed via PCR amplification and subsequent sequencing, resulting in the establishment of D. cerebrum mitfa^−/− colonies.

To generate the Danionella cerebrum transgenic line, Tg(ubb^R:jGCaMP8m), a Tol2 vector was constructed containing the ubb^R promoter^87,88, provided by Balciunas and Lazutka. Following the promoter, a zebrafish codon-optimized jGCaMP8m sequence⁸⁸ and the SV40 polyadenylation signal were arranged sequentially. This plasmid (25 ng μl⁻¹), along with Tol2 transposase mRNA (20 ng μl⁻¹), was co-injected (0.5 nl total injected volume) into one-cell stage mitfa^−/− embryos. Embryos at 3 days post-fertilization were screened, and those exhibiting strong jGCaMP8m expression were reared to maturity. Upon reaching adulthood, these fish were screened collectively, and the most robust founders were selected for further study.

Experimental procedures

Sample preparation for functional imaging

Before imaging, zebrafish and D. cerebrum samples were embedded on their right side in a drop of 2% low-melting point agarose (A9414, Sigma) in a glass-bottom Petri dish (35 mm, P35G-1.5-14-C, Mattek). The right pectoral fin was moved away from the flank of the fish to point either tangentially or anteriorly such as to prevent it from covering visceral organs. Following agarose solidification, the dish was filled with E3 fish water. Samples were left to rest and settle for 30 min before the imaging session to minimize drift during the experiment. Cutaneous respiration is considered to be able to meet the oxygen needs of young zebrafish⁸⁹ and thus agarose was not removed from the gills, nor mouth, and no oxygen perfusion was required.

Microscope and data acquisition

Functional data were acquired using a Nikon spinning-disk inverted confocal microscope (CSU-W1) with a ×20 0.75 NA air objective (field of view of 800  μm × 600 μm, working distance of 1 mm) or a ×40 1.15 NA water immersion objective. For dual-colour acquisition of GCaMP7f and mCherry signals, excitation lasers at 488 nm and 594 nm were used (3–4% and 3–5% power). Emitted light was split using a 560 long-pass dichroic and passed through emission filters (525/36, 610 LP) before reaching the cameras (Hamamatsu ORCA-Fusion BT). For jRGECO1b imaging, a 561-nm laser line was used with a 610/75-nm emission filter. Camera exposure times were set between 120 ms and 200 ms. A piezo motor was used to acquire fast z-stacks with a typical inter-plane interval of 7–9 μm, resulting in a volumetric scan rate of approximately 0.3 Hz. Smaller z-steps were used for more targeted investigations (for example, during vascular imaging). There is a trade-off between the speed of volumetric imaging and sampling density (number of planes acquired), with the upper limit set by the sensor (for example, GCaMP), which acts as a low-pass filter. We have provided a quantitative means of estimating the present and resolvable temporal frequency content across tissues (Extended Data Fig. 13).

Ketamine treatment

Animals were imaged for a 25-min baseline period, after which ketamine was manually added to the imaging dish to achieve a final concentration of 100  μg ml⁻¹ (400 μM).

Tricaine treatment

Animals were treated with 750 μM tricaine (MS-222, E10521-10G, Sigma) diluted in E3 water. The samples were incubated with the drug for 15 min before and throughout the duration of the experiments.

Cold stimulus delivery

Fish water was cooled to 10 °C, and 350 μl of cold water was manually added to the imaging dish after 5 min of baseline imaging. To rule out any neural influence or motion, fish were anaesthetized before and during the experiment with tricaine (750 μM; MS-222, E10521-10G, Sigma).

Hypoxia treatment

Oxygen levels were programmatically varied using an Okolab O₂ controller module, controlled via the Nikon spinning-disk microscope software (Elements). The module varies O₂ levels by altering the ratio of N₂ and O₂ mixed and the oxygen levels at the chamber inlet were recorded. A baseline period (21% O₂) of at least 5 min was recorded, after which oxygen levels were lowered to 10% for 10–20 min, after which levels were returned to 21%. Oxygen levels were measured in the bath using an oxygen micro-optode (O₂ MicroOptode, Unisense) inserted in the agarose (Fig. 4a). We fit a mono-exponential model to the data using the scipy.optimize.curve_fit function. The model accounts for 86% of the variance with a time constant of approximately 1.2 min. We note that hypoxia is expected to cause changes in acid–base balance within cells⁹⁰. As the aim was to study the physiological consequences of hypoxia, we did not attempt to control tissue pH. We further note that even if the agarose is 98% water, and commonly used in the field, given that water is not flowing, this might still result in oxygen levels that are slightly lower than if there was no agarose, and if the animals were freely swimming. We therefore consider our assay to be one in which the change in oxygen levels are what is most important.

Optogenetics

A commercial integrated digital micromirror device (DMD) module within the Nikon spinning-disk confocal microscope was used for all optogenetic experiments in conjunction with a 488-nm LED. The duty cycle of the DMD was set to 10%, and LED power set to approximately 8–10%. A dual-colour reference image stack was acquired at the beginning of the experiment to get the anatomical location of the cells expressing CoChR–eGFP. On the basis of this 3D volume, regions of interest (ROIs) to be illuminated were defined in 3D via the graphical user interface and programmatically stored. For optogenetic activation of motor vagal cells (Fig. 2), stimulus intensity was initially calibrated by collecting a dose–response curve with progressively increasing stimulus intensity until a reliable response was observed in the neurons being stimulated (approximately 8% LED power). Once calibrated, the experimental protocol was programmed using Nikon software (Elements v5) with inter-stimulus intervals of 10 s, with ROIs selected randomly (choose without repeat) or sequentially. ROIs were illuminated for varying durations ranging from 10 to 500 ms during dose–response experiments, and for all other experiments for one frame (approximately 180 ms). The data and metadata were extracted using the Python package ND2. For each trial type (that is, specific ROI stimulated), the difference in activity pre-stimulus and post-stimulus was extracted and averaged over an approximately 5-s window (two volumes). To display the data within a single graph, the averaged stimulation-induced maps were combined, with each pixel coloured according to the ROI inducing the most change in fluorescence, alpha-weighted by the magnitude of change with alpha set to 0 when the change was below noise level (estimated by the average change induced by the control ROI). For optogenetic inhibition of the hindbrain (Fig. 4), the same stimulus intensity was used as in the optogenetic activation experiment and was confirmed to result in loss of all motor output. After 5 min of onset of lower oxygen levels (10%), an ROI over the hindbrain or the heart (control) was illuminated in alternation for 1 min followed by 1 min of no illumination. After three alternating trials of each location, oxygen levels were returned to 21%.

Data processing workflow

Registration

The registration pipeline was based on local iterative motion estimation.

Iterative patch-wise optical flow. We referred to the static image as the ‘template image’, I_temp, and referred to the data that we aimed to transform as the ‘moving image’, I_mov. The key part of the registration algorithm is the iterative patch-wise optical flow module. This module estimates motion using a modified version of the optical flow algorithm¹². It aims to minimize the intensity difference between the template and the moving image while encouraging the motion field to vary smoothly over space. We formulated an optimization problem defined by the following loss function (for clarity, the one-dimensional version is described; in practice, this is extended to three dimensions):

where ΔX ∈ R^N is the motion field to be estimated, Δxⁿ is the nth element of ΔX, N is the total number of patches the image is split into, Z is the number of pixels per patch, Δxⁿ is the estimate of the spatial displacement for patch n and \(\overline{\Delta {x}^{n}}\) is the average Δxⁿ for the neighbouring patches of n.

This original optimization problem is non-convex and cannot be efficiently solved. Therefore, we re-formulated the problem into one that can be solved iteratively through a set of convex problems using a modified version of the optical flow algorithm¹²:

Computing spatial intensity gradients can be prone to noise if done pixel-wise; to overcome this, a locally averaged estimate of the spatial intensity gradient was used. The image was divided into patches of ([2r + 1] × [2r + 1]) pixels and the average spatial gradient was computed (r = 5 pixels). The accumulated motion \(\Delta {x}_{a}^{n}\) was estimated and optimized per patch. Once convergence was reached, the motion between patch centres was linearly interpolated. The estimated displacement is also desired to be somewhat smooth over space. Therefore, the loss function was augmented with a penalty term weighted by β, which penalizes the difference between each motion estimate of a patch’s motion and the average of the immediately adjacent patches. This individual optimization step is a quadratic function and thus has a closed-form solution. Once a motion step ΔX_s was computed, the accumulated motion was updated and linearly interpolated between patch centres, then the accumulated motion was applied to the moving frame, and this was iterated until the loss ceased to decrease beyond a user-defined threshold or reached the maximum iteration number of each pyramid layer N_iter. To take into account that it is the accumulated motion that ought to be spatially smooth, the terms weighted by β are the accumulated patch and patch neighbourhood motion. Furthermore, to enable the parallelization of calculations, the gradients of the previous time steps were used for the neighbourhood motion estimation.

Image processing scheme and registration reliability mask. First, one needed to define the template and the moving image. For short datasets, the first frame served as the template, and each subsequent frame was treated as the moving image. This setup allowed for the alignment of all frames one by one to the initial frame. For longer-duration datasets, where bleaching occurs, a floating template was introduced along with an initialized motion field. The floating template was defined as the median of the last N_template selected motion-corrected frames, whereas the initial motion field was chosen as the one closest to the centre of these N_template motion fields.

To enhance the algorithm’s robustness to noise, the ‘foreground’ in both the template and the current frame was defined. This step excluded regions containing moving immune cells, which appear as small, bright, connected components that otherwise introduce distortions into the estimated motion field. After pre-processing the data, pyramid downsampling was performed, a technique that reduces image resolution by iteratively smoothing and subsampling the image. This reduced the size of the template, moving image and the initialized motion field to 1/2^L of their original dimensions in the x and y directions, where L is the number of levels. The original size along the z direction was maintained. On the basis of specified patch size parameters, iterative patch-wise optical flow registration was performed on the downsampled data, which gives a rough estimate of the motion field. After this, the data were upsampled and the motion field further refined, thereby enhancing its accuracy. Once fit, using the refined motion field, the calcium channel was corrected by relocating the moved cells to their original positions as seen in the first frame. This systematic approach provides robust frame alignment and motion correction, accommodating dynamic cell movements and variations in image intensity throughout the functional recording. Finally, to filter out any data that were either not sufficiently well or unreliably registered, a registration reliability mask was computed. To identify unreliable textures, characterized by low spatial intensity gradients, the amplitude of the spatial gradient of the template image was computed and standardized using z-scores. Regions exhibiting a score below 0 had little texture and were flagged as potentially unreliable; in addition, the mean squared error was computed for each patch across time, and all patches whose error was above a user-set threshold were discarded from downstream analysis.

Parameters and implementation. The algorithm was executed with β = 0.01 (smoothness parameters), r = 5 (motion correction patch size, measured in pixels), N_L = 4 (number of pyramid layers used for motion correction), N_template = 5 (number of frames incorporated in the floating template), N_moving = 5 (number of frames to consider when initializing motion) and N_iter = 10 (maximum iteration number of each pyramid layer). The code was developed and run in Python and MATLAB (R2024b) and is available on GitHub (https://github.com/vruetten/wholistic_registration, https://github.com/Weizheng96/WHOLISTIC-registration).

Validation. To validate the registration at single-cell resolution, the transgenic line used for the registration, pancellular membrane marker (Tg(β-actin2:mCherry-CAAX)) was crossed to a sparse transgenic line, (Tg(phox2bb:eGFP)), to use the latter as held-out ground truth. Dual-colour data were acquired for 1 h. Using the reference channel alone, the motion field was estimated and applied to the held-out green (sparse) channel. In both unregistered and registered data, individual cells within the sparse held-out channel were manually located in the ‘template’ (first time point) and 20 and 60 min into the recording within Fiji software. Displacement distance was subsequently extracted for individual cells.

Registration: method comparison

Synthetic dataset generation

To benchmark registration performance under controlled but challenging conditions, we created a library of three-dimensional simulated recordings in which the signal-to-noise ratio, motion smoothness and displacement amplitude were varied independently.

Seed volume. We selected an anatomically rich volume of 256 × 256 × 13 voxels and termed it the reference stack (I_ref).

Motion field synthesis. A zero-mean, unit-variance three-dimensional Gaussian random field G(x, y, z) was convolved with an isotropic Gaussian kernel of standard deviation σ_k. The smoothed field was normalized to unit peak magnitude and scaled by an amplitude factor A to yield the displacement field. Finally, we added a bias term R to simulate rigid deformation to the final motion field.

Larger σ_k produces more spatially coherent motion, whereas increasing A and R increases the maximum voxel displacement.

Noise injection and simulation parameters. Additive white Gaussian noise \(\eta \sim {\mathcal{N}}(0,{\sigma }_{\eta })\) was applied independently to the reference and the warped (moving) image: 9 noise levels σ_η = (1.4⁰…1.4⁸), 16 σ_k = (5–20), and 10 amplitudes A = (1–10 pixels) were explored. When one parameter was swept, the remaining two were fixed at σ_η = 5, σ_k = 20 and A = 5 pixels. R was a 2D vector, of length 15 pixels added to each plane, with direction randomized between realizations. Ten stochastic realizations were generated for every parameter combination and averaged to obtain a robust estimate.

Performance evaluation

For every simulated dataset, we measured two scores:

1. (1)

   Image fidelity: the mean squared error (MSE) between I_ref and the registered image, \(\widehat{I}\).

2. (2)

   Motion fidelity: the MSE between the ground-truth displacement field Δ and the field recovered by the algorithm \(\widehat{\Delta }\).

Data along the borders of the image (20 pixels) were not included in the MSE as none of the registration methods can extrapolate meaningfully and so borders were discarded in downstream analysis to avoid artefacts.

Benchmark algorithms

WHOLISTIC registration was compared against four widely used non-rigid registration frameworks, each tuned for volumetric calcium-imaging data. The results are presented in Extended Data Fig. 2.

WHOLISTIC registration. Method parameters were set to pyramid_layer = 3, smooth_penalty = 0.01 and r  = 5. In addition, we have provided the results without the pyramid registration.

Suite2p (MATLAB GPU port). Multi-plane registration was enabled; the field of view was partitioned into 40 × 40 xy blocks with 2-pixel overlap.

NoRMCorre. Parameters were set to grid_size = [16, 16, 1], mot_uf = [4, 4, 1], max_shift = [15, 15, 5], max_dev = [3, 3, 1] and overlap_pre = [4, 4, 1].

Demons (SimpleITK). Intensity histograms were matched (1,024 bins, 16 match points and background threshold = mean). The algorithm ran for 200 iterations with Gaussian-smoothed updates (σ = 0.5).

B-spline FFD (SimpleITK). A 24 × 24 × 13 control-point grid was optimized with L-BFGS-B (tol = 10⁻⁵, 100 iterations, 5 corrections, 1,000 evaluations; cost-function convergence factor = 10⁷). Correlation similarity and linear interpolation were used; the final transform was converted to a dense displacement field.

Cellular segmentation

In dense volumetric fluorescence imaging, camera pixels inevitably receive photons from multiple overlapping sources due to the optical point spread function and tissue scattering. To take into account the mixed sources, we used non-negative matrix factorization, a linear source separation method that explicitly models the data as the superposition of multiple spatially overlapping sources with distinct temporal dynamics.

The registered data were segmented on a cellular scale using Voluseg, a pipeline that implements locally constrained non-negative matrix factorization to transform neighbouring correlated pixels into functional segments¹³ (https://github.com/mikarubi/voluseg/). In brief, an intensity-based brain mask was created and divided the volume into spatially contiguous three-dimensional blocks, which overlap slightly to capture cells on the borders, and cell detection was run in parallel on these through the use of Spark, a distributed cluster-computing framework. The algorithm fits the following model:

where V is the full spatiotemporal fluorescence matrix for each block, W and H are, respectively, the spatial footprint and time series of segmented cells, and X and I are rank 1 spatiotemporal model of the background signal. The time series were then normalized (mean-subtracted, and divided by the standard deviation of the time series). This formulation explicitly represents each pixel’s fluorescence as a weighted linear combination of cellular sources contributing to that pixel, enabling computational separation of overlapping sources. This resultant normalized time series H is denoted F_norm.

Denoising

Synchronous Ca²⁺ bursts arising from muscle activation result in significant fluorescence. Camera pixels that ought to only receive light from non-muscle cells can, if situated in close proximity to muscle tissue, capture some scattered photons from muscle cells, which confounds downstream analysis. Although the contribution of a constant baseline fluorescence from muscle is effectively mitigated by mean subtraction of the data, fluctuations in such fluorescence remain problematic.

To address this challenge, the timepoints of muscle activity are located and nearby highly correlated cells identified. These potentially contaminated data points were masked, and linear interpolation was applied to the time series. This method may omit activity genuinely linked to muscle activations but is a conservative approach that could be refined in future work, or two-photon microscopy can be used instead of one-photon confocal to avoid out-of-focus excitation. More specifically, the functional tissue ensembles corresponding to muscle were labelled, separating (owing to per-plane variations in contamination) and averaging them by imaging plane to form ‘muscle group’ time series. Muscle activation timepoints are determined using the probabilistic oasis model from Suite2p⁹¹ with parameters (window = maxmin, win baseline = 120 s and sig baseline = 2). Activation timepoints were almost always isolated, that is, neighbouring timepoints contained no muscle activation by virtue of swimming being sparse in time in these experiments. Activation events with a probability over 0.6 are deemed real, and correlated cells above 0.35 close to the muscle group (within one plane above or below the muscle group) were masked at these timepoints, and the signal was then linearly interpolated between neighbouring time points.

Resolvable and present spatial-frequency estimation via FRC

To quantify the resolvable and present spatial frequencies across imaged volumes, we computed the Fourier ring correlation (FRC)⁹² of pairs of independent images of the same object (x, y).

The FRC is defined as the real-valued, normalized cross-power:

Here k denotes a radial shell in Fourier space, and F₁(f_x, f_y) and F₂(f_x, f_y) are the complex Fourier coefficients of the two independent images at frequencies f_x, f_y. Writing each coefficient in polar form \({F}_{j}=| {F}_{j}| \,{e}^{i{\phi }_{j}}\) and using \({\rm{Re}}[{F}_{1}{F}_{2}^{* }]=| {F}_{1}| | {F}_{2}| \cos ({\phi }_{1}-{\phi }_{2})\), we obtain:

This expression highlights that the FRC is the weighted mean cosine of the phase differences in shell k, with weights given by the product of Fourier magnitudes, normalized so that −1 ≤ FRC ≤ 1.

FRC maps were computed for volumes acquired dorsally and sagittally at 0.1625 × 0.1625 resolution. To ensure that the results were not limited by the true spatial frequency content of the specimen, we used a transgenic line labelling all cellular membranes, guaranteeing high-frequency content throughout most of the sample. We note that organs such as the ear, swim bladder and gallbladder containing acellular spaces have no high-frequency content and thus, regardless of the achievable resolution, will have a lower FRC score. The analysis was run with a patch size of 20 μm with 75% overlap. The spatial frequency at which the FRC fell below 1/7 was taken as the resolution cut-off, following established practice⁹², and was mapped across the sample (user adjustable).

Resolvable and present temporal-frequency estimation via FTC

To quantify the resolvable and present temporal frequencies in imaged volumes, we generalized the FRC to the time domain, which we term Fourier temporal correlation (FTC). Instead of taking two independent images of the same sample, we considered independent samples of a time series by splitting the time series into odd and even bins, x₁ and x₂.

We define FTC as the real-valued, normalized cross-spectrum:

where \({S}_{{x}_{1},{x}_{2}}\) is the cross-spectrum between time series x₁ and x₂ obtained by splitting the original time series x into odd and even frames. The spectrum is estimated using Welch’s method, with K being the number of time windows the time series is split into. Like FRC, −1 ≤ FTC ≤ 1.

Sensor bandwidth. We note that the published single-spike half-decay times (t_1/2 = 0.27 s)⁸⁸ indicate \(\tau ={t}_{1/2}/\text{ln}2=0.39\,{\rm{s}}\). The fluorescence impulse response of jGCaMP7f can be approximated by a mono-exponential kernel \(h(t)=\frac{1}{\tau }\,{e}^{-t/\tau }\), which has Fourier magnitude \(| H(f)| =\frac{1}{\sqrt{1+{(2\pi f\tau )}^{2}}}\). Substituting f = 3.5 Hz gives ∣H(3.5 Hz)∣_{τ = 0.39} = 0.12, that is, a nearly 10× loss in amplitude or ≈100× loss in power. Consequently, frequencies above approximately 3.5 Hz are already attenuated by nearly one order of magnitude.

Coherence-based spectral clustering for identification of functional tissue ensembles

To cluster the data, spectral clustering¹⁴ was performed using a new coherence-based distance measure to define the adjacency graph. The similarity function was defined to be: \({W}_{i,j}=\exp \,\left(-\frac{1-{\mathcal{C}}({{\boldsymbol{x}}}_{i},{{\boldsymbol{x}}}_{j})}{\tau }\right)\) where \({\mathcal{C}}({{\boldsymbol{x}}}_{j},{{\boldsymbol{x}}}_{j})(\omega )\) is the coherence between unit x_i and x_j. This can be demonstrated to be a valid positive semi-definite kernel. Coherence is estimated using Welch’s method, that is, interpreting it as the windowed magnitude-squared coherence estimator for stationary signals:

where K is the number of time windows.

This definition of coherence weights each \({\mathcal{C}}(\omega )\) at each frequency irrespective of the power in that frequency band, this can be detrimental as bands of low power can be dominated by noise. To avoid this, the measure was normalized by the total power across all frequency bands:

which remains bounded between 0 and 1. Having defined this adjacency graph, standard spectral clustering was performed. The Laplacian of the adjacency matrix was computed: L_norm = I − D⁻¹W where D is a diagonal matrix and d_i,i = ∑_jW_i,j. The lowest N eigenvectors of the matrix were computed and used as an encoding basis, and finally the k-means algorithm was run to identify clusters. The k-means clusters were initialized with k-means++⁹³, an algorithm for choosing good initializations of the cluster centroids. The model was fit with: N_clusters = 400, τ = 0.3, for coherence estimation, a time window of approximately 8 min was used, with 80% overlap and a Hanning tapering window.

We chose to use coherence rather than correlation as the similarity metric to accommodate phase lagged or delayed activity patterns that are prevalent in our data (Extended Data Fig. 5). 

Hierarchical clustering and lag-regression model

Spectral clusters mapping to muscle were identified manually based on their anatomy and time series, and the mean activity of each cluster was computed. Hierarchical clustering was performed on the cluster means using the linkage function from the Python SciPy package, utilizing correlation as a metric and Ward linkage. A dendrogram was plotted, with muscle groups with greater than 0.5 correlation colour coded in similar shades, which we have denoted hyper-clusters, and shown in anatomical space with the same colour scheme. To confirm the greater synchrony of the cervical epaxial muscle with the ventral abdominal muscle, rather than with the neighbouring hypaxial muscle, the cluster corresponding to the cervical epaxial muscle was identified in different samples and cross-correlation between the activity of hypaxial and abdominal muscle computed, and a one-sided Wilcoxon signed-rank test used.

The mean activity of the muscle hyper-clusters identified in hierarchical cluster analysis was used as the basis for the regressors in the lag-regression model. A lag-regression model is a predictive model for time-series data in which a regression equation is used to predict the current value of the dependent variable, y(t), based on both the current and the past (that is, lagged), values of an explanatory variable, x(t), x(t − 1),⋯, x(t − L) where L is the maximum number of lags. This is the equivalent to assuming that the observed dependent variable arises from the convolution of the independent variable with a learnt kernel, which is parametrized by a set of weights, one for each lag:

where M is the number of regressors or independent variables, and L is the maximum number of lags.

The kernel parameters (k_m) are fit to minimize the mean squared reconstruction error (that is, \(\parallel y-\bar{y}{\parallel }_{2}^{2}\), where y is a vector containing the observations, and \(\bar{y}\) is the model prediction), by finding the least squares solution using the numpy.linalg.solve function. This approach allows each tissue to have different temporal response properties (for example, neurons show fast responses, whereas other tissues tend to show slower, more prolonged responses to motor activity).

The mean activity of the muscle hyper-clusters (grouped clusters or cellular ensembles) were converted into a weighted binary trace by identifying the onset of motor contraction and weighting them by the power of the contraction. Contraction onset and offset were identified by finding positive deflections above the noise floor and the time at which the curve returned to the noise floor; power was defined by the area under the curve of these two time points. For each cell y_i ∈ Y, a lag-regression model was fit, \({y}_{i}(t)={\sum }_{m=1}^{M}{\sum }_{\tau =0}^{L}{x}_{m}(t-\tau )\cdot {k}_{m}(\tau )+{\epsilon }\) where M is the number of muscle regressors included (4–6), L is the number of time lags included (60 s), y_i(t) and x_m(t) are the activity of a cell i (y_i), and muscle regressor m (x_m), at time t. A time lag of 60 s was chosen, as we wanted to capture the relatively short-lived responses. The identified kernel values were threefold cross-validated and the average R² value on held-out data is reported.

Kymograph analysis

Peak oscillation frequency was defined as the frequency with largest power excluding the zero-frequency power. To compute kymographs of ependymal cell activity from imaging data acquired from Tg(foxj1a:GCaMP7f) imaging, the ventral and dorsal boundaries of the brain were first anatomically identified. Line integrals perpendicular to the ventral boundary were computed, resulting in a 1 × N vector with N being the number of spatial bins used. When repeated over frames, this results in a N × T matrix. The data were denoised through spatial smoothing using a Gaussian kernel with a standard deviation of 4 pixels and temporally filtered with a bandpass Butterworth filter of order two, using cut-off periodicities of 0.5–9 min.

Coherence k-means algorithm

Below is the derivation of the algorithm used in Fig. 3k. The k-means algorithm alternates between computing cluster centroids and identifying the clusters to which data points belong. We used coherence as a measure of similarity. We define the following: M as the number of latent clusters, K as the number of sub-samples used to average over to calculate coherence, T as the total number of time points, L as the number of time points per sub-sample, x_i ∈ R^T as sample i, x_ik ∈ R^L as sub-sample k of sample i of length L, \({{\boldsymbol{\mu }}}_{m}\in {{\mathbb{R}}}^{L}\) as the centroid of cluster m, and \({\mathcal{C}}(x,y)(\omega )\) as the coherence between x and y at frequency ω.

We wished to find a centroid μ_m such that the sum of the coherences of points within that cluster \({\sum }_{i}{\mathcal{C}}({{\boldsymbol{x}}}_{i},{{\boldsymbol{\mu }}}_{m})\) is maximal. This quantity is invariant to an arbitrary real scaling of \({\widehat{\mu }}_{mk}\), so the following constraint was added: S_μμ(ω) = 1, that is, \({\mathbb{E}}[| {\widehat{\mu }}_{m}(\omega ){| }^{2}]=1\) so \(\frac{1}{K}{\sum }_{k=i}^{K}| {\widehat{\mu }}_{mk}(\omega ){| }^{2}=1\).

Thus, we aimed to maximize:

where \({\lambda }_{i}=\frac{1}{{K}^{2}{S}_{{x}_{i}{x}_{i}}(\omega ){S}_{\mu \mu }(\omega )}=\frac{1}{{K}^{2}{S}_{{x}_{i}{x}_{i}}(\omega )}\)

We note that \({A}_{k,{k}^{{\prime} }}(\omega )={A}_{{k}^{{\prime} },k}^{* }(\omega )\) so the elements form a Hermitian matrix A(ω).

Vectorizing μ_m, the objective can be written to maximize as:

from which it follows that \({\widehat{{\boldsymbol{\mu }}}}_{m}(\omega )\) is the leading eigenvector of A(ω).

To initialize the group labels, the k-means++ algorithm⁹³ was first run to choose initial centroid values. Cluster centroids were calculated as above, and cluster allocation reassigned the cluster ID based on coherence. The process was iterated until the convergence of the labels.

Characterizing hypoxia-induced physiological changes

Estimation of blood vessel diameter. To measure blood vessel diameter, a transgenic line labelling vascular endothelium (Tg(flk1:dsRED-CAAX)) was imaged. To ensure accurate width estimation, volumetric stacks were acquired with 2-μm z-spacing spanning the entire width of the mesenteric blood vessel, and maximum intensity projection of each stack was computed. A line crossing the middle of the vessel orthogonally was defined, and the kymograph was computed using Fiji’s kymograph function, resulting in a matrix with T rows, where T is the number of time points. To denoise the data, the matrix was row-wise median filtered with a window size of 5 pixels. The matrix was thresholded at the midrange value. The largest connected component was identified using the scipy-ndimage label function, which corresponds to the main blood vessel. The width of the mask at each row was computed, and the resulting time series was smoothed over time using a rolling mean with a window length of 10 time points. Finally, the width was scaled by the resolution of the data to have units of microns. Baseline and hypoxic-state vessel widths were defined as the mean width over the 1-min interval preceding the onset of the gas switch or at steady-state hypoxia (10 min following the gas switch from 21% to 10% O₂), respectively.

Estimation of blood flow. To measure blood flow, a transgenic line labelling red blood cells (Tg(gata1:dsRED)) was imaged. For faster imaging, a single plane was acquired at 10 Hz. The plane was selected to cover the mesenteric artery, which runs parallel to the body. The frame rate was too low to track individual red blood cells; instead, local bulk blood flow was estimated by quantifying changes in fluorescence induced by red blood cells moving in and out of any specific ROI. ROIs over arteries feeding the brain, muscle and mesentery were manually defined. The absolute value temporal difference of the mean was computed for each of these ROIs. As passing red blood cells result in differing changes in magnitude depending on whether the entire or a fraction of the cell was in the ROI, the time series was clipped between 0 and 3× the median value of the beginning of the time series (normoxic period, 2 min). The time series was finally normalized by dividing by the standard deviation of this baseline period. The baseline and hypoxic-state blood flow were defined as the mean blood flow over the 1-min interval preceding the onset of the gas switch or at steady-state hypoxia (10 min following the gas switch from 21% to 10% O₂), respectively. Change in blood flow was defined as the difference between these values and baseline blood flow (mean blood flow before gas switch, approximately 2 min). To compare the effects of optogenetic inhibition of the hindbrain and control illumination of the heart, mean blood flow over each 1-min stimulation period was computed for each condition (hindbrain or control).

Baseline oxygen modulation score. To quantify the effect of hypoxia on baseline calcium levels, activity traces were passed through a minimum–maximum filter using centred windows of length 3 min, and were subsequently smoothed using a Gaussian kernel with a standard deviation of 15 s. The mean baseline activity during normoxia and hypoxia over a 3-min window (taken at the end of the normoxic–hypoxic phases to ensure steady-state oxygen levels had been reached) was computed on a per-cell basis and subtracted.

WB-ExM

Comprehensive, step-by-step protocols are available on protocols.io: WHOLISTIC ExM: whole-body expansion microscopy with immunofluorescence and histological stains (https://doi.org/10.17504/protocols.io.dm6gp9wxjvzp/v5) and WHOLISTIC ExM: whole-body expansion microscopy with fluorescence in situ hybridization (WB-ExM FISH; https://doi.org/10.17504/protocols.io.5qpvo9y89v4o/v1).

WB-ExM-IF

Fixation, immunofluorescence and agarose embedding. Whole larval zebrafish were fixed with 4% paraformaldehyde overnight at 4 °C on a shaker, then washed with 1× PBS four times for 15 min. Fixed fish were permeabilized for 5 h with 0.5% Triton X-100 in 1× PBS (PBST-0.5) at room temperature with gentle agitation. Permeabilization must be included for all specimens, whether or not they were stained with antibodies. For immunofluorescence, specimens were blocked for 3 h with blocking buffer (5% goat serum, 0.5% Triton X-100 and 0.1% Na-azide in 1× PBS) at room temperature with gentle agitation. Blocked specimens were incubated with primary antibodies diluted 1:100 in blocking buffer for 3 days at room temperature, followed by washing with PBST-0.5 three times for 2 h. Specimens were then incubated with secondary antibodies diluted 1:100 in blocking buffer for 2 days at room temperature, followed by washing with PBST-0.5 three times for 2 h. Samples were embedded in a thin layer of 1% low-melting temperature agarose to ensure the desired sample orientation. The following antibodies were used: rabbit anti-eGFP (A11122, Invitrogen), chicken anti-RFP (409006, Synaptic Systems), goat anti-rabbit Atto647N (40839, Sigma) and donkey anti-chicken Alexa568 (A78950, Invitrogen).

Protein anchoring. Specimens were incubated in Acryloyl-X SE at 20 µg ml⁻¹ in 1× PBS (anchoring solution) for 1 h at room temperature followed by washing in 1× PBS. Anchoring solution was prepared just before use from a 10 mg ml⁻¹ stock solution dissolved in anhydrous DMSO. Samples were washed in 1× PBS three times for 5 min.

Gelation and digestion. Specimens were incubated in the first gelation solution (10% acrylamide, 0.5 M sodium acrylate, 0.1% bis-acrylamide, 0.01% 4HT, 0.2% TEMED, 0.2% APS and 1× PBS) on ice three times for 10 min with gentle agitation. Gelation chambers were constructed using an uncharged glass slide as the bottom piece and side walls consisting of 11 layers of Scotch tape (approximately 600 μm thick to approximate the thickness of the agarose block) serving as spacers. Coverslips were placed on top of the chamber. The chambers were filled with the first gelation solution by pipetting in solution from the open side. Fully assembled chambers were carefully placed in a humidified incubator for 2 h at 37 °C for gelation. After gelation, the chambers were carefully disassembled, and extra gel was trimmed around the samples using a scalpel, leaving an approximately 2-mm margin around the fish. Gelled specimens were treated with disruption buffer (5% SDS, 50 mM Tris pH 7.5 and 200 mM NaCl in H₂O) at 100 °C overnight in Eppendorf tubes, then washed with 1× PBS three times for 20 min.

Re-embedding and staining. Gelled and disrupted specimens were incubated in the second monomer solution (10% acrylamide, 0.5 M sodium acrylate, 0.02% bisacrylamide, 0.01% 4HT, 0.2% TEMED, 0.2% APS and 1× PBS; this is the same as the first gelation solution except with bisacrylamide reduced from 0.1% to 0.02%) on ice for 3 × 10 min with gentle agitation. Gels imbued with the second gelation solution were placed on uncharged glass slides with side spacers as before, composed this time of 20 layers of Scotch tape (approximately 1.2 mm). A coverslip was positioned on top as the chamber top piece. The space surrounding the first gel was filled with additional second gelation solution. Fully assembled chambers were carefully placed in a humidified incubator for 2 h at 37 °C. After gelation, the chambers were gently disassembled, and excess gel was trimmed around the embedded fish with a scalpel, leaving an approximately 1-mm margin. Gelled specimens were stained with Alexa 488-NHS ester (A20000, Thermo Fisher Scientific) and/or Atto647N-maleimide (2857, AAT Bioquest) dyes. Samples were incubated in dye solution (1:1,000 in PBS from 10 mg ml⁻¹ stocks dissolved in anhydrous DMSO (D12345, Thermo Fisher)) for 1 h at room temperature with shaking, then washed three times for 1 h with 1× PBS. The re-embedding process can be repeated up to four times resulting in approximately 5× expansion.

Imaging, tile stitching and data registration. Samples were mounted onto a glass slide using polylysine and attached to a custom-built sample holder. Specimens were imaged on a Zeiss Z1 light-sheet microscope (×20 NA 1.0 water immersion objective), immersed in 1× PBS. Stitching and registration were done using the Imaris stitching and registration software using the default parameters.

Online data. An example dataset can be viewed online (https://neuroglancer-demo.appspot.com/#!gs://flyem-user-links/short/2025-01-03.165221.071748.json). This is the sample shown in Fig. 3g: double transgenic zebrafish animal labelling the ventricular and vascular systems (Tg(foxj1a:eGFP) × Tg(flk1:dsRED-CAAX)), stained against eGFP (magenta) and dsRED (green; 10 days post-fertilization, expanded approximately 2×).

WB-ExM-FISH

All solutions were made using molecular grade (RNase-free) water and reagents. A comprehensive, step-by-step protocol is provided in Supplementary Notes.

Anchoring stock solutions. Melphalan and Acryloyl-X SE (AcX) were prepared as 2.5 mg ml⁻¹ and 10 mg ml⁻¹ stocks, respectively, in anhydrous DMSO and stored desiccated. Melphalan-X was prepared by combining melphalan and AcX stocks at a ratio of 4:1 to yield a final concentration of 2 mg ml⁻¹ each and incubating at room temperature overnight with shaking. Melphalan-X was aliquoted and stored desiccated.

Fixation and permeabilization. Whole larval zebrafish were fixed with 4% paraformaldehyde overnight at 4 °C on a shaker, then washed 4 × 15 min in 1× PBS. Fixed fish were permeabilized for 1 h with PBST-0.5 at room temperature with gentle agitation.

RNA and protein anchoring. Fixed and permeabilized fish were washed with MOPS buffer (20 mM, pH 7.7) for 30 min at room temperature. Specimens were next treated with 1 mg ml⁻¹ melphalan-X supplemented with 0.1 mg ml⁻¹ extra AcX diluted freshly into MOPS buffer overnight at 37 °C, followed by washing in MOPS buffer (2 × 5 min) and then PBS (2 × 5 min). Anchored specimens were then mounted on poly-L-lysine-coated coverslips with the fish positioned on its side.

Gelation and digestion. Mounted specimens were incubated in the first gelation solution and gelled as in the protein method variant above. In brief, this included incubating in complete gelation solution, assembling the gelation chamber, gelling at 37 °C, disassembling the chamber and trimming the gel to leave an approximately 2-mm margin around the specimen. The trimmed gel was then incubated in the first digestion solution (500 mM NaCl, 0.3% SDS, 50 mM Tris-HCl pH 8.0 and 1 mM EDTA with proteinase K diluted 1:50 from 800 U ml⁻¹ stock) for 4 h at 50 °C, followed by incubation in the second digestion solution (50 mM NaCl, 1% SDS, 50 mM Tris-HCl pH 8.0 and 1 mM EDTA with proteinase K diluted 1:50 from 800 U ml⁻¹ stock) overnight at 50 °C. Digested specimens were washed 4 × 15 min in 1× PBS.

Re-embedding. Digested specimens were re-gelled as in the protein method above. Gelled specimens were recovered into 1× PBS and trimmed, leaving an approximately 1-mm margin around the specimen.

Probe hybridization and hybridization chain reaction. Gels were incubated in hybridization buffer (Molecular Instruments; https://www.molecularinstruments.com/hcr-rnafish-products) for 30 min at 37 °C, followed by incubation in primary probes designed by Molecular Instruments (6 μl in 600 μl hybridization buffer, 10 nM final concentration) overnight at 37 °C. Following hybridization, gels were washed with pre-warmed (37 °C) probe wash buffer (Molecular Instruments) 3 × 30 min, then washed with pre-warmed (37 °C) 1× PBS 3 × 1 h, and once overnight at room temperature. Gels were next incubated in amplification buffer (Molecular Instruments; https://www.molecularinstruments.com/hcr-rnafish-products) for at least 30 min at room temperature. Hairpins were diluted 1:50 in amplification buffer and snap cooled by heating to 95 °C for 90 s followed by cooling at room temperature for 30 min. Gels were incubated for 4 h at room temperature in the dark in hairpin–amplification buffer mix. Amplified gels were washed in 5× SSCT (5× SSC and 0.1% Tween) 2 × 20 min at room temperature, then in 0.5× SSCT (0.5× SSC and 0.1% Tween) 2 × 40 min at room temperature and finally equilibrated in 1× PBS for 2× expansion. The following probes were used: phox2bb_B3 (lot #RTG113), Th_B1 (lot #RTB474), hsd3b1_B5 (lot #RTG103), calca_B2 (lot #RTG105) and pomca_B5 (lot #RTG108).

Stripping and re-probing. For multi-round imaging, hybridization chain reaction amplification products and probes were stripped by digestion with DNase followed by re-probing and imaging as described for the first round, above.

Imaging, tile stitching and data registration. Samples were processed the same way as WB-ExM-IF samples.

Whole-body gel embedding to preserve endogenous fluorescence

Fixed fish were embedded in 2% low-melting-point agarose and permeabilized with 0.1% saponin in PBS overnight at 4 °C. Permeabilized fish were treated with AcX gel anchor solution (20 μg ml⁻¹ Acryloyl-X-NHS in 1× PBS for 1 h at room temperature, diluted freshly from 10 mg ml⁻¹ stock in anhydrous DMSO). Fish were incubated in 1× gelation solution (10% acrylamide, 5% N,N′-diallyltartardiamide, 1× PBS, 0.05% APS and 0.05% TEMED) on ice with shaking for 30 min, followed by gelation at 37 °C for 2 h under humidified nitrogen, with a coverslip placed on top of the agarose block. Embedded fish were stained in 1× PBS with DAPI and 0.3% saponin for 2 days.

PhotoMap

To quantify deformation introduced by expansion in an unbiased manner, we developed PhotoMap, a method that measures the gel’s deformation field from a pre-imposed reference pattern, inspired by GelMap³². In brief, a fluorescent gel (Extended Data Fig. 10a) was photobleached with a regular grid pattern using a two-photon microscope before expansion (Extended Data Fig. 10b), and re-imaged post-expansion (Extended Data Fig. 10c). The resulting grid deformation can be readily visualized and quantified (Extended Data Fig. 10d,e). A comprehensive step-by-step protocol is provided in Supplementary Notes.

Fluorescent gel formation. A fluorescent monomer was created by conjugating AF488-NHS (20 mg ml⁻¹, 31 mM) to 3-aminopropyl methacrylamide (10 mg ml⁻¹, 56 mM) in the presence of triethylamine (10%, 720 mM). The reagents were incubated overnight at room temperature, protected from light. Gel preparation followed the protocol described in WB-ExM-IF, except that 1:100 of the fluorescent monomer was added to the first monomer solution to render the gel intrinsically fluorescent. In brief, samples were fixed, agarose mounted, chemically anchored and polymerized, and then kept unexpanded in the gelation chamber for photobleaching.

PhotoMap grid formation through photobleaching. Using a large field-of-view two-photon microscope, a 50 × 50 × 50 μm grid was photobleached into the gel before expansion. Custom code was written to drive the microscope via ScanImage⁹⁴. Regions corresponding to the eyes, which contain light-absorbing pigments, were excluded from the photobleaching pattern to prevent excessive local heating. After photobleaching, the sample was disrupted overnight, expanded and re-imaged using the same microscope at 1 × 1 × 5 μm resolution.

For PhotoMap deformation field analysis, to extract the resulting grid intersections from post-expansion images, we used a two-step procedure combining zero-mean normalized cross-correlation (ZNCC) with peak detection.

For template matching via ZNCC, we first manually selected a representative intersection from the raw image and used this as a template \(T\in {{\mathbb{R}}}^{P\times Q}\). The full image \(I\in {{\mathbb{R}}}^{N\times M}\) was then scanned using ZNCC, defined as:

where μ_I and μ_T are the local mean intensities of the image and the template, respectively, over the ROI. This was efficiently computed via fast Fourier transform (FFT)-based cross-correlation and using sliding windows.

For local peak detection, the resulting ZNCC map was passed through a local top-percentile filter, keeping the brightest 10% of pixels per patch (80 × 80 pixels). A Laplacian-of-Gaussian filter was then applied to enhance bright-on-dark blobs corresponding to grid intersections. Peak centres were finally extracted using the peak_local_max function from scikit-image, with minimum distance (50 pixels) and relative threshold parameters empirically tuned to match the grid geometry.

For quantifying local grid distortion, to assess spatial distortions in the photobleached grid post-expansion, two geometric features were computed for each grid point:

• Length deviation: the mean relative deviation of edge lengths (horizontal and vertical) from an ideal spacing L = 90 pixels, computed as \(1-| | {\rm{vec}}(v)| /L| \).

• Angle deviation: the cosine of the angle between adjacent horizontal and vertical vectors at each point, ideally zero for orthogonal axes (that is, deviation from 90°).

These values were computed across the entire sample. To visualize and compare the spatial distortion statistics, kernel density estimates of each component (length and angle) were computed separately for on-sample and off-sample points. The analysis pipeline of PhotoMap is available (https://github.com/vruetten/PhotoMap.git).

A comprehensive, step-by-step protocol is available on protocols.io: PhotoMap: unbiased mapping of expansion microscopy deformation fields (https://doi.org/10.17504/protocols.io.261ge169wv47/v1).

ExM data modelling and quantification

Building 3D ExM model. To build the main whole-fish model, WB-ExM-Histo zebrafish datasets were acquired (total protein stain Alexa488-NHS) and imaged on a Zeiss Z1 light-sheet microscope (×20 NA 1.0 water immersion objective). The organs were manually segmented within the Amira software and the meshes imported in the software Blender. To incorporate cellular populations that have too complex morphology or spatial distribution for manual annotation, and that can be molecularly defined, such as ependymal cells and the vasculature, WB-ExM-IF data were utilized. To build meshes, WB-ExM-IF data were up-sampled in z by a factor of 2, converted to 8-bit format, denoised (Fiji’s Remove Outliers function), thresholded (Fiji’s Threshold function, Otsu algorithm), opened in Fiji’s 3D viewer as a surface and exported as an stl binary file, ready to be imported into Blender. Meshes from different samples were manually aligned in Blender to a common reference space using the total protein stain present in all datasets.

Quantification of signal-to-noise ratio in ExM data. Scattering across the sample was quantified by comparing signal-to-noise ratio within muscle fibres between the near and far-side of the sample, defined as the average difference between minimum and maximum intensity signals across muscle sarcomeres.

Quantification of enteric motor vagus innervation. To estimate motor vagus innervation density along the gastrointestinal tract, an animal expressing RFP under the isl1 promoter (Tg(isl1CREST-hsp70l:mRFP)) was stained and expanded using WB-ExM-IF. The gastrointestinal tract was manually segmented using the Amira software. The immunofluorescence signal was averaged radially resulting in a one-dimensional histogram of average innervation density as a function of position along the gastrointestinal tract.

Tracing of motor nerves. A high-resolution confocal stack at 4 days post-fertilization (Tg(VAChTa:eGFP)) was acquired using a spinning-disk confocal microscope (0.1625 μm × 0.1625 μm × 1 μm voxel size). Nerve bundles and muscle outlines were traced manually using Amira, and the resulting tracts were imported into Blender.

Statistics and reproducibility

Images in Figs. 1b,d,g,h, 2h,i, 3g,i and 4b,h,k,n–o and Extended Data Figs. 1a,c,d,f,g, 4a–e, 6b, 7a,h, 8e–j, 9a,b,d–k, 11d–e and 12a–c are representative examples of data collected in at least four animals.

Ethics statement

All animal procedures were approved by the Institutional Animal Care and Use Committee of the Howard Hughes Medical Institute, Janelia Research Campus and were conducted in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Raw WHOLISTIC data (https://s3.janelia.org/ahrens-lab/ruettenv/WHOLISTIC/) and raw WB-ExM-immunohistochemistry data (https://s3.janelia.org/ahrens-lab/ruettenv/ExM/) are available. The data can also be browsed on the accompanying website (https://wholistic.janelia.org/); an example WB-ExM-immunohistochemistry dataset can be viewed (https://wholistic.janelia.org/data/exm-demo-dataset-vasculo-ventricular-system/). The plasmid map for Tg(ubb^R:jGCaMP8m) is available from Addgene (plasmid 232485).

1. Clapham, D. E. Calcium signaling. Cell 131, 1047–1058 (2007).

   Article  CAS  PubMed  Google Scholar 

2. Chen, F., Tillberg, P. W. & Boyden, E. S. Expansion microscopy. Science 347, 543–548 (2015).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

3. Buchman, T. G. The community of the self. Int. Congress Ser. 1255, 3–5 (2003).

   Article  Google Scholar 

4. Sofroniew, N. J., Flickinger, D., King, J. & Svoboda, K. A large field of view two-photon mesoscope with subcellular resolution for in vivo imaging. eLife 5, e14472 (2016).

   Article  PubMed  PubMed Central  Google Scholar 

5. Schulze, L. et al. Transparent Danionella translucida as a genetically tractable vertebrate brain model. Nat. Methods 15, 977–983 (2018).

   Article  CAS  PubMed  Google Scholar 

6. Dana, H. et al. High-performance calcium sensors for imaging activity in neuronal populations and microcompartments. Nat. Methods 16, 649–657 (2019).

   Article  CAS  PubMed  Google Scholar 

7. Mosimann, C. et al. Ubiquitous transgene expression and Cre-based recombination driven by the ubiquitin promoter in zebrafish. Development 138, 169–177 (2011).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

8. Li, F. et al. Generation of gcamp6s expressing zebrafish to monitor spatiotemporal dynamics of calcium signaling elicited by heat stress. Int. J. Mol. Sci. 22, 5551 (2021).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

9. Gossen, M. & Bujard, H. Tight control of gene expression in mammalian cells by tetracycline-responsive promoters. Proc. Natl Acad. Sci. USA 89, 5547–5551 (1992).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

10. Ahrens, M. B., Orger, M. B., Robson, D. N., Li, J. M. & Keller, P. J. Whole-brain functional imaging at cellular resolution using light-sheet microscopy. Nat. Methods 10, 413–420 (2013).

    Article  CAS  PubMed  Google Scholar 

11. Tomer, R., Khairy, K., Amat, F. & Keller, P. J. Quantitative high-speed imaging of entire developing embryos with simultaneous multiview light-sheet microscopy. Nat. Methods 9, 755–763 (2012).

    Article  CAS  PubMed  Google Scholar 

12. Horn, B. K. & Schunck, B. G. Determining optical flow. Artif. Intell. https://doi.org/10.1016/0004-3702(81)90024-2 (1981).

13. Mu, Y. et al. Glia accumulate evidence that actions are futile and suppress unsuccessful behavior. Cell 178, 27–43.e19 (2019).

    Article  CAS  PubMed  Google Scholar 

14. Von Luxburg, U. A tutorial on spectral clustering. Stat. Comput. 17, 395–416 (2007).

    Article  MathSciNet  Google Scholar 

15. Teixido-Trujillo, S. et al. Measured GFR in murine animal models: review on methods, techniques, and procedures. Pflügers Arch. Eur. J. Physiol. 475, 1241–1250 (2023).

    Article  CAS  Google Scholar 

16. Szebényi, K. et al. Visualization of calcium dynamics in kidney proximal tubules. J. Am. Soc. Nephrol. 26, 2731–2740 (2015).

    Article  PubMed  PubMed Central  Google Scholar 

17. Jimenez, Y. E., Parsons, J. W. & Brainerd, E. L. Epaxial and hypaxial co-contraction: a mechanism for modulating strike pressure and accuracy during suction feeding in channel catfish. J. Exp. Biol. 226, jeb244714 (2023).

    Article  PubMed  Google Scholar 

18. Tulenko, F. J. & Currie, P. Zebrafish myology. In The Zebrafish in Biomedical Research: Biology, Husbandry, Diseases, and Research Applications (eds Cartner, S. C. et al.) 115–121 (Academic Press, 2020).

19. Schilling, T. F. & Kimmel, C. B. Musculoskeletal patterning in the pharyngeal segments of the zebrafish embryo. Development 124, 2945–2960 (1997).

    Article  Google Scholar 

20. Jimenez, Y. E. & Brainerd, E. L. Dual function of epaxial musculature for swimming and suction feeding in largemouth bass. Proc. R. Soc. B 287, 20192631 (2020).

    Article  PubMed  PubMed Central  Google Scholar 

21. Wang, N. et al. Mechanotransduction pathways in articular chondrocytes and the emerging role of estrogen receptor-α. Bone Res. 11, 13 (2023).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

22. Jurisch-Yaksi, N., Yaksi, E. & Kizil, C. Radial glia in the zebrafish brain: functional, structural, and physiological comparison with the mammalian glia. Glia 68, 2451–2470 (2020).

    Article  PubMed  Google Scholar 

23. Dana, H. et al. Sensitive red protein calcium indicators for imaging neural activity. eLife 5, e12727 (2016).

    Article  PubMed  PubMed Central  Google Scholar 

24. Ran, C., Boettcher, J. C., Kaye, J. A., Gallori, C. E. & Liberles, S. D. A brainstem map for visceral sensations. Nature 609, 320–326 (2022).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

25. Kaneko, T., Boulanger-Weill, J., Isabella, A. J., & Moens, C. B. Position-independent functional refinement within the vagus motor topographic map. Cell Rep. 43, 114740 (2024).

26. Li, X. et al. The periaxonal space as a conduit for cerebrospinal fluid flow to peripheral organs. Proc. Natl Acad. Sci. USA 121, e2400024121 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

27. Leung, L. C. et al. Neural signatures of sleep in zebrafish. Nature 571, 198–204 (2019).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

28. Steib, E. et al. TissUExM enables quantitative ultrastructural analysis in whole vertebrate embryos by expansion microscopy. Cell Rep. Methods 2, 100311 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

29. Sim, J. et al. Nanoscale resolution imaging of whole mouse embryos using expansion microscopy. ACS Nano 19, 7910–7927 (2025).

30. Wang, K. et al. TSA-PACT: a method for tissue clearing and immunofluorescence staining on zebrafish brain with improved sensitivity, specificity and stability. Cell Biosci. 13, 97 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

31. Mao, C. et al. Feature-rich covalent stains for super-resolution and cleared tissue fluorescence microscopy. Sci. Adv. 6, eaba4542 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

32. Damstra, H. G. J. et al. GelMap: intrinsic calibration and deformation mapping for expansion microscopy. Nat. Methods 20, 1573–1580 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

33. Deng, S. et al. Roles of ependymal cells in the physiology and pathology of the central nervous system. Aging Dis. https://doi.org/10.14336/AD.2022.0826-1 (2022).

34. Furube, E. et al. Neural stem cell phenotype of tanycyte-like ependymal cells in the circumventricular organs and central canal of adult mouse brain. Sci. Rep. 10, 2826 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

35. Miranda-Negrón, Y. & García-Arrarás, J. E. Radial glia and radial glia-like cells: their role in neurogenesis and regeneration. Front. Neurosci. 16, 1006037 (2022).

    Article  PubMed  PubMed Central  Google Scholar 

36. Zhang, J. et al. Wnt-PLC-IP3-Connexin-Ca²⁺ axis maintains ependymal motile cilia in zebrafish spinal cord. Nat. Commun. 11, 1860 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

37. Luo, Z. et al. Hypoxia signaling in human health and diseases: implications and prospects for therapeutics. Signal Transduct. Target. Ther. 7, 218 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

38. Berna, N., Arnould, T., Remacle, J. & Michiels, C. Hypoxia-induced increase in intracellular calcium concentration in endothelial cells: role of the Na⁺-glucose cotransporter. J. Cell. Biochem. 84, 115–131 (2002).

    Article  Google Scholar 

39. Szabo, J. S., Stonestreet, B. S. & Oh, W. Effects of hypoxemia on gastrointestinal blood flow and gastric emptying in the newborn piglet. Pediatr. Res. 19, 466–471 (1985).

    Article  CAS  PubMed  Google Scholar 

40. Eliason, E. J. & Farrell, A. P. Effect of hypoxia on specific dynamic action and postprandial cardiovascular physiology in rainbow trout (Oncorhynchus mykiss). Comp. Biochem. Physiol. A Mol. Integr. Physiol. 171, 44–50 (2014).

    Article  CAS  PubMed  Google Scholar 

41. Reid, J. L. α-Adrenergic receptors and blood pressure control. Am. J. Cardiol. 57, 8E–12E (1988).

    Google Scholar 

42. Li, J. et al. Intron targeting-mediated and endogenous gene integrity-maintaining knockin in zebrafish using the CRISPR/Cas9 system. Cell Res. 25, 634–637 (2015).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

43. Cannon, W. B. The Wisdom of the Body (W.W. Norton, 1939).

44. Alon, U. Systems Medicine: Physiological Circuits and the Dynamics of Disease (CRC Press, 2023).

45. Dantzer, R. Neuroimmune interactions: from the brain to the immune system and vice versa. Physiol. Rev. 98, 477–504 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

46. Billman, G. E. Homeostasis: the underappreciated and far too often ignored central organizing principle of physiology. Front. Physiol. 11, 200 (2020).

    Article  PubMed  PubMed Central  Google Scholar 

47. Palmquist, K. H., Ko, C. S., Shyer, A. E. & Rodrigues, A. R. Biological theories of morphogenesis based on holistic biophysical thinking. Biol. Theory https://doi.org/10.1007/s13752-024-00477-1 (2024).

48. Vladimirov, N. et al. Light-sheet functional imaging in fictively behaving zebrafish. Nat. Methods 11, 883–884 (2014).

    Article  CAS  PubMed  Google Scholar 

49. Koberstein, J. N. et al. Monitoring glycolytic dynamics in single cells using a fluorescent biosensor for fructose 1,6-bisphosphate. Proc. Natl Acad. Sci. USA 119, e2204407119 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

50. Wang, H. et al. A tool kit of highly selective and sensitive genetically encoded neuropeptide sensors. Science 382, eabq8173 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

51. Liu, W., Liu, C., Ren, P. G., Chu, J. & Wang, L. An improved genetically encoded fluorescent cAMP indicator for sensitive cAMP imaging and fast drug screening. Front. Pharmacol. 13, 902290 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

52. Lobas, M. A. et al. A genetically encoded single-wavelength sensor for imaging cytosolic and cell surface ATP. Nat. Commun. 10, 711 (2019).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

53. Abdelfattah, A. S. et al. Bright and photostable chemigenetic indicators for extended in vivo voltage imaging. Science 365, 699–704 (2019).

    Article  ADS  CAS  PubMed  Google Scholar 

54. Farrants, H. et al. A modular chemigenetic calcium indicator for multiplexed in vivo functional imaging. Nat. Methods 21, 1916–1925 (2024).

55. Allier, C. et al. Decomposing heterogeneous dynamical systems with graph neural networks. Preprint at https://doi.org/10.48550/arXiv.2407.19160 (2025).

56. Park, J., Pao, G. M., Sugihara, G., Stabenau, E. & Lorimer, T. Empirical mode modeling: a data-driven approach to recover and forecast nonlinear dynamics from noisy data. Nonlinear Dyn. 108, 2147–2160 (2022).

    Article  Google Scholar 

57. Chang, T. Y., Pardo-Martin, C., Allalou, A., Wählby, C. & Yanik, M. F. Fully automated cellular-resolution vertebrate screening platform with parallel animal processing. Lab Chip 12, 711–716 (2012).

    Article  CAS  PubMed  Google Scholar 

58. Hovorka, R. Closed-loop insulin delivery: from bench to clinical practice. Nat. Rev. Endocrinol. 7, 385–395 (2011).

    Article  CAS  PubMed  Google Scholar 

59. Pardo-Martin, C. et al. High-throughput hyperdimensional vertebrate phenotyping. Nat. Commun. 4, 1467 (2013).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

60. Wang, Y. et al. EASI-FISH for thick tissue defines lateral hypothalamus spatio-molecular organization. Cell 184, 6361–6377.e24 (2021).

    Article  CAS  PubMed  Google Scholar 

61. Dorkenwald, S. et al. Neuronal wiring diagram of an adult brain. Nature 634, 124–138 (2024).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

62. Tavakoli, M. R. et al. Light-microscopy-based connectomic reconstruction of mammalian brain tissue. Nature 642, 398–410 (2025).

63. Januszewski, M. et al. High-precision automated reconstruction of neurons with flood-filling networks. Nat. Methods 15, 605–610 (2018).

    Article  CAS  PubMed  Google Scholar 

64. Zhou, G. et al. An artificial intelligence accelerated virtual screening platform for drug discovery. Nat. Commun. 15, 7761 (2024).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

65. Do, A., Zahrawi, F. & Mehal, W. Z. Therapeutic landscape of metabolic dysfunction-associated steatohepatitis (MASH). Nat. Rev. Drug Discov. https://doi.org/10.1038/s41573-024-01084-2 (2024).

66. Gendelev, L. et al. Deep phenotypic profiling of neuroactive drugs in larval zebrafish. Nat. Commun. 15, 9955 (2024).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

67. Patton, E. E., Zon, L. I. & Langenau, D. M. Zebrafish disease models in drug discovery: from preclinical modelling to clinical trials. Nat. Rev. Drug Discov. 20, 611–628 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

68. Stringer, C. et al. Rastermap: a discovery method for neural population recordings. Nat. Neurosci. 28, 201–212 (2025).

69. Westerfield, M. The Zebrafish Book. A Guide for the Laboratory Use of Zebrafish (Danio rerio) 5th edn (Univ. Oregon Press, 2007).

70. Anderson, J. L. et al. Multiple sex-associated regions and a putative sex chromosome in zebrafish revealed by RAD mapping and population genomics. PLoS ONE 7, e40701 (2012).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

71. White, R. M. et al. Transparent adult zebrafish as a tool for in vivo transplantation analysis. Cell Stem Cell 2, 183–189 (2008).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

72. Kawakami, K. Transposon tools and methods in zebrafish. Dev. Dyn. 234, 244–254 (2005).

    Article  CAS  PubMed  Google Scholar 

73. Horstick, E. J. et al. Increased functional protein expression using nucleotide sequence features enriched in highly expressed genes in zebrafish. Nucleic Acids Res. 43, e48 (2015).

    Article  PubMed  PubMed Central  Google Scholar 

74. Kwan, K. M. et al. The Tol2kit: a multisite gateway based construction kit for Tol2 transposon transgenesis constructs. Dev. Dyn. 236, 3088–3099 (2007).

    Article  CAS  PubMed  Google Scholar 

75. Espinosa-Medina, I. et al. TEMPO enables sequential genetic labeling and manipulation of vertebrate cell lineages. Neuron 111, 345–361.e10 (2023).

    Article  CAS  PubMed  Google Scholar 

76. Park, H.-C. et al. Analysis of upstream elements in the HuC promoter leads to the establishment of transgenic zebrafish with fluorescent neurons. Dev. Biol. 227, 279–293 (2000).

    Article  CAS  PubMed  Google Scholar 

77. Mohammad, F. et al. Optogenetic inhibition of behavior with anion channelrhodopsins. Nat. Methods 14, 271–274 (2017).

    Article  CAS  PubMed  Google Scholar 

78. Fowler, D. K. et al. A multisite gateway toolkit for rapid cloning of vertebrate expression constructs with diverse research applications. PLoS ONE 11, e0159277 (2016).

    Article  PubMed  PubMed Central  Google Scholar 

79. Taniguchi, A., Kimura, Y., Mori, I., Nonaka, S. & Higashijima, S. I. Axially confined in vivo single cell labeling by primed conversion using blue and red lasers with conventional confocal microscopes. Dev. Growth Differ. 59, 741–748 (2017).

    Article  PubMed  PubMed Central  Google Scholar 

80. Fujita, M. et al. Assembly and patterning of the vascular network of the vertebrate hindbrain. Development 138, 1705–1715 (2011).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

81. Traver, D. et al. Transplantation and in vivo imaging of multilineage engraftment in zebrafish bloodless mutants. Nat. Immunol. 4, 1238–1246 (2003).

    Article  CAS  PubMed  Google Scholar 

82. Yang, E. et al. A brainstem integrator for self-location memory and positional homeostasis in zebrafish. Cell 185, 5011–5027.e20 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

83. Yamaguchi, N. et al. Rear traction forces drive adherent tissue migration in vivo. Nat. Cell Biol. 24, 194–204 (2022).

    Article  MathSciNet  CAS  PubMed  PubMed Central  Google Scholar 

84. Grant, P. K. & Moens, C. B. The neuroepithelial basement membrane serves as a boundary and a substrate for neuron migration in the zebrafish hindbrain. Neural Dev. 5, 9 (2010).

    Article  PubMed  PubMed Central  Google Scholar 

85. Nechiporuk, A., Linbo, T., Poss, K. D. & Raible, D. W. Specification of epibranchial placodes in zebrafish. Development 134, 611–623 (2007).

    Article  CAS  PubMed  Google Scholar 

86. Grimes, D. T. et al. Zebrafish models of idiopathic scoliosis link cerebrospinal fluid flow defects to spine curvature. Science 352, 1341–1344 (2016).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

87. Bakūnaitė, E. et al. Recombinant ubb^R promoter enables highly efficient tamoxifen-inducible Cre recombination in embryonic and adult zebrafish. Genetics https://doi.org/10.1093/genetics/iyag155 (2026).

88. Zhang, Y. et al. Fast and sensitive GCaMP calcium indicators for imaging neural populations. Nature 615, 884–891 (2023).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

89. Jonz, M. G. Respiratory system. In The Zebrafish in Biomedical Research: Biology, Husbandry, Diseases, and Research Applications (eds Cartner, S. C. et al.) 103–107 https://doi.org/10.1016/B978-0-12-812431-4.00010-5 (Academic Press, 2020).

90. Swenson, E. R. The many acid-base manifestations and consequences of hypoxia. Curr. Opin. Physiol. 7, 72–81 (2019).

    Article  Google Scholar 

91. Pachitariu, M. et al. Suite2p: beyond 10,000 neurons with standard two-photon microscopy. Preprint at bioRxiv https://doi.org/10.1101/061507 (2017).

92. Koho, S. et al. Fourier ring correlation simplifies image restoration in fluorescence microscopy. Nat. Commun. 10, 3103 (2019).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

93. Arthur, D. & Vassilvitskii, S. k-means++: the advantages of careful seeding. In Proc. 18th Annual ACM-SIAM Symposium on Discrete Algorithms 1027–1035 (SIAM, 2007).

94. Pologruto, T. A., Sabatini, B. L. & Svoboda, K. ScanImage: flexible software for operating laser scanning microscopes. Biomed. Eng. Online 2, 13 (2003).

    Article  PubMed  PubMed Central  Google Scholar 

95. Ruetten, V. Whole body model. GitHub https://github.com/vruetten/wholebodymodel (2025).

96. Ruetten, V. & Chi, Y. WHOLISTIC registration. GitHub https://github.com/vruetten/wholistic_registration (2025).

97. Zheng, W. WHOLISTIC registration. GitHub https://github.com/Weizheng96/WHOLISTIC-registration (2025).

98. Ruetten, V. & Chi, Y. WHOLISTIC registration: WHOLISTIC Nature paper (version v1.0.0) [computer software]. Zenodo https://doi.org/10.5281/zenodo.21460339 (2026).

99. Rubinov, M., Tauffer, L., Xu, C. & Dichter, B. Voluseg. GitHub https://github.com/mikarubi/voluseg (2018).

100. Rubinov, M., Tauffer, L., Xu, C. & Dichter, B. Voluseg: WHOLISTIC Nature paper (version v1.0.0) [computer software]. Zenodo https://doi.org/10.5281/zenodo.21460399 (2026).

101. Ruetten, V. PhotoMap. GitHub https://github.com/vruetten/PhotoMap (2025).

102. Ruetten, V. PhotoMap: WHOLISTIC Nature paper (version v1.0.0) [computer software]. Zenodo https://doi.org/10.5281/zenodo.21460352 (2026).

Download references

Author notes

1. These authors jointly supervised this work: Maneesh Sahani, Misha B. Ahrens

Authors and Affiliations

1. Janelia Research Campus, HHMI, Ashburn, VA, USA

   Virginie M. S. Ruetten, Igor Siwanowicz, Brett D. Mensh, Mark Eddison, Amy Hu, Andrew L. Lemire, Caiying Guo, Yisheng He, Kari Close, Gudrun Ihrke, Aparna Dev, Alyson Petruncio, Yinan Wan, Rongwei Zhang, Mikail Rubinov, Philipp J. Keller, Chie Satou, Paul W. Tillberg & Misha B. Ahrens

2. Gatsby Computational Neuroscience Unit, UCL, London, UK

   Virginie M. S. Ruetten & Maneesh Sahani

3. Virginia Tech, Blacksburg, VA, USA

   Wei Zheng & Guoqiang Yu

4. Tsinghua University, Beijing, China

   Yunfeng Chi & Guoqiang Yu

5. Charité Universitätsmedizin Berlin, Berlin, Germany

   Mykola Kadobianskyi, Marc Renz & Benjamin Judkewitz

6. Harvard University, Cambridge, MA, USA

   Sara Lelek-Greskovic, Mark C. Fishman & Florian Engert

7. Biozentrum, University of Basel, Basel, Switzerland

   Yinan Wan

8. Vanderbilt University, Nashville, TN, USA

   Mikail Rubinov

Contributions

V.M.S.R., M.S. and M.B.A. conceptualized the study. V.M.S.R. developed the methodology, undertook the data collection and performed the data analysis. V.M.S.R., W.Z., Y.C., M. Rubinov, and G.Y. conducted registration and segmentation. I.S. built the 3D model and performed the rendering. V.M.S.R., M.E., A.H., Y.H., K.C., G.I., A.P., A.D. and P.W.T. conducted ExM. V.M.S.R., C.G., A.L.L., C.S., M.K., M. Renz, B.J., Y.W. and P.J.K. performed the genetics studies and created transgenics. S.L.-G., R.Z., F.E. and M.C.F. contributed to conceptualization and testing the hypoxia findings. V.M.S.R., B.D.M., M.S. and M.B.A. wrote, reviewed and edited the manuscript. M.S. and M.B.A. supervised the study and acquired funding.

Corresponding authors

Correspondence to Virginie M. S. Ruetten or Misha B. Ahrens.

Competing interests

V.M.S.R. and M.B.A. are co-inventors on a provisional patent application no. 64/011,419, filed by the Howard Hughes Medical Institute relating to the WHOLISTIC imaging method described in this work. The remaining authors declare no competing interests.

Peer review information

Nature thanks the anonymous reviewers for their contribution to the peer review of this work. Peer reviewer reports are available.

Reporting Summary (download PDF )

Peer Review File (download PDF )

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions