Imaging cellular activity across all organs reveals body-wide circuits

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.