Whole-brain serial-section electron microscopy in larval zebrafish

Animal care

Adult zebrafish (Danio rerio) for breeding were maintained at 28 °C on a 14 hr:10 hr light:dark cycle following standard methods³¹. The Tg(elavl3:GCaMP5G)a4598 transgenic line³² used in this study was of genotype elavl3:GCaMP5G^+/+; nacre (mitfa^−/−), conveying nearly pan-neuronal expression of the calcium indicator GCaMP5G³³ and increased transparency due to the nacre mutation³⁴. The larval zebrafish samples described in this study were raised in filtered fish facility water³¹ until 5–7 days post-fertilization (dpf).

Mice from which support tissue was collected had been previously euthanized for other experiments. Only unused, to-be-discarded tissue was harvested to serve as support tissue.

The Standing Committee on the Use of Animals in Research and Training of Harvard University approved all animal experiments.

Two-photon laser-scanning microscopy

Larval zebrafish were immobilized by immersion in 1 mg mL⁻¹ α-bungarotoxin (Invitrogen) and mounted dorsum-up in 2% low-melting-temperature agarose in a small dish containing a silicone base (Sylgard® 184, Dow Corning). Upon agarose hardening, E3 solution (5 mM NaCl, 0.17 mM KCl, 0.33 mM CaCl₂, and 0.33 mM MgSO₄) was added to the dish. In vivo structural imaging of elavl3-driven GCaMP5G signal was conducted with a custom-built two-photon microscope equipped with a Ti:Sapphire laser (Mai Tai®, Spectra-Physics) excitation source tuned to 800 nm. Frames with a 764.4×509.6 µm² field of view size (1200×800 px²) were acquired at 1 µm intervals (0.637×0.637×1 µm³vx⁻¹) at ~1 Hz with a scan pattern of four evenly spaced, interlaced passes³⁵. A low-noise anatomical snapshot of brain fluorescence was captured in 300 planes, each the sum of 50 single frames. All light-based imaging was performed without any intentional stimulus presentation.

Dissection and tissue preparation

Initial attempts at high-quality larval zebrafish brain preservation were impeded by skin and membranes, which prevented sufficient fixation with whole-fish immersion alone (Extended Data Fig. 1a). To overcome this, the skin and membranes covering the brain³⁶ were dissected away.

Each larval zebrafish, previously immobilized and embedded for two-photon laser-scanning microscopy, was introduced to a dissection solution (64 mM NaCl, 2.9 mM KCl, 10 mM HEPES, 10 mM glucose, 164 mM sucrose, 1.2 mM MgCl₂, 2.1 mM CaCl₂, and pH 7.5; ref. 37) containing 0.02% (w/v) tricaine mesylate (MS-222, Sigma-Aldrich). Flow of red blood cells through the vasculature was confirmed before proceeding as an indicator of good health. A portion of agarose was removed to expose the dorsum from the posterior hindbrain to the anterior optic tectum. The dissection was initiated by puncturing the thin epithelial layer over the rhombencephalic ventricle above the hindbrain³⁸ with a sharpened tungsten needle. Small incremental anterior-directed incisions were made along the midline as close to the surface as possible until the brain was exposed from the hindbrain entry to the middle of the optic tectum (Extended Data Fig. 1b). The majority of damage associated with this dissection was restricted to medial tectal proliferation zone progenitor cells³⁹ that are unlikely to have integrated into functional neuronal circuits.

Dissections lasted 1–2 min, upon which time the complete dish was immersed in a 2.0% formaldehyde and 2.5% glutaraldehyde fixative solution (Electron Microscopy Sciences) overnight at room temperature (Extended Data Fig. 1d). Following washes, larval zebrafish were cut out from the dish in a block of agarose with a scalpel and moved to a round-bottom microcentrifuge tube. Specimens were then incubated in post-fixation solution containing 1% osmium tetroxide and 1.5% potassium ferricyanide for 2 hr (Extended Data Fig. 1e), washed with water, washed with 0.05M maleate buffer (pH 5.15), and stained with 1% uranyl acetate in maleate buffer overnight (Extended Data Fig. 1f). During the subsequent wash step with maleate buffer, larval zebrafish were freed from the surrounding agarose block and moved to a new microcentrifuge tube. Next, specimens were washed with water, dehydrated with serial dilutions of acetonitrile in water (25%, 50%, 70%, 70%, 80%, 90%, 95%, 100%, 100%, 100%) for 10 min each, and infiltrated with serial dilutions of a diepoxyoctane-based low viscosity resin⁴⁰ in acetonitrile (25%, 50%, 75%, 100%) for 1 hr each. The samples were then embedded in the diepoxyoctane-based resin with surrounding support tissue and hardened for 2–3 days at 60 °C (Extended Data Fig. 1g–h). Aqueous solutions were prepared with water passed through a purification system (typically Arium 611VF, Sartorius Stedim Biotech). This process resulted in high-quality ultrastructure preservation (Extended Data Fig. 1i).

Additional solution, washing, and timing details were described previously in a step-by-step protocol⁴¹.

Serial sectioning

Consistent ultrathin sectioning was difficult to achieve in samples containing heterogeneous tissues but imperative for reconstructing three-dimensional (3-D) structure from a series of two-dimensional (2-D) sections. Tests revealed that errors occur primarily when the sample composition changed dramatically (e.g., borders between tissue and empty resin). We overcame this by embedding samples within a surrounding support tissue of mouse cerebral cortex (Extended Data Figs. 1g–h,2f).

We preferred sectioning perpendicular to most axon and dendrite paths for ease and reliability in reconstructing neuronal morphology. For this, our cutting plane was oriented perpendicular to the long (anterior-posterior) axis, which required ~2.5× more sections than alternative orientations. This was made possible by customising an automated tape-collecting ultramicrotome²⁶ by extending the device’s main mounting plate and enlarging its reels (compare Extended Data Fig. 2a with Fig. 1E from ref. 26) to accommodate one long tape stretch capable of collecting all sections.

Sections were continuously cut with a diamond knife (Extended Data Fig. 2b–c) affixed to an ultramicrotome (EM UC6, Leica) and collected onto 8mm-wide and 50–75 µm-thick tape (Kapton® polyimide film, DuPont). Restarts were occasionally required for three reasons: fine-tuning of tape positioning or settings is necessary at the beginning of a run; the ultramicrotome design is constrained by a cutting depth range of ~200 µm; and diamond knives must be shifted after cutting several thousand sections to expose the sample to a fresh edge before dulling impairs sectioning quality. When necessary, restarts were completed as quickly as possible (typically 1–2 min) to minimize possible thermal, electrostatic, or other fluctuations. For the same reason, tape reels were fed continuously without ever being reloaded or exchanged. This combination of fast restarts and continuous tape feeding successfully maintained a steady state across restarts.

We sectioned two larval zebrafish specimens and these represent the only two samples we have attempted to cut with the surrounding support tissue approach. The primary focus of this study was a 5.5 dpf larval zebrafish sectioned with a 45° ultra diamond knife (Diatome) and a nominal sectioning thickness that averaged 60 nm with a variable setting ranging from 50–70 nm depending on sectioning consistency. Restarts occurred after sections 276, 3,669, 6,967, 10,346, 12,523, 12,916, and 15,956. Knife shifts occurred after sections 6,967 and 12,916. After sectioning, the tape was cut into segments with a razor blade between collected sections and adhered with double-sided conductive carbon adhesive tape (Ted Pella) to 4 in-diameter silicon wafers (University Wafer), which served as an imaging substrate. A total of 17,963 × ~60 nm-thick sections were spread across 80 wafers (Extended Data Figs. 2d–e, 3).

One potential limitation of the 5.5 dpf larval zebrafish series is the section thickness. Minimizing section thickness is an important factor in the success of axon and dendrite reconstructions¹. Small neuronal processes (on the order of the section thickness) are difficult to reconstruct in thicker sections, especially when they are running roughly parallel to the plane of the section. To be sure that our approach was not fundamentally limited to thicker sections, we sectioned the second sample—a 7 dpf larval zebrafish—with a nominal sectioning thickness that remained constant at 50 nm throughout the entire cutting session using a 45° histo diamond knife (Diatome). Restarts occurred after sections 296, 312, 4,114, 8,233, and 12,333. Knife shifts occurred after sections 4,114 and 12,333. A total of 15,046 × ~50 nm-thick sections were obtained from 15,052 attempted (Extended Data Fig. 4) and spread across 70 wafers. The thinner sections did not result in more lost material: this series contained 6 losses (0.04%; Extended Data Fig. 4d upper), 25 partial sections (0.17%; Extended Data Fig. 4d middle), no adjacent losses, and 6 adjacent lost-partial or partial-partial events (0.04%; Extended Data Fig. 4d lower).

The nominal section thickness of ~60 nm made it possible to span the entire 5.5 dpf larval zebrafish brain in ~18,000 sections, as determined by finding the location of the spinal cord-hindbrain boundary⁴². Though the 7 dpf sample was sectioned at 50nm, it was not made the focus of subsequent imaging because of an over-trimming error that caused less of the brain to be captured. However, improved reliability for this sample despite a ~17% reduction in nominal sectioning thickness suggests that yet higher axial resolution is attainable. A section thickness of ≤30 nm would increase confidence in the ability to reconstruct complete neuronal circuit connectivity, and thicknesses of ≤30 nm are known to be possible for mammalian brain sections of comparable sizes^{27, 43} using a similar approach.

Once wafers contained tape segments, they were made hydrophilic by glow discharging very briefly, post-section stained for 1–2 min inside a chamber containing sodium hydroxide pellets using a stabilised lead citrate solution (UltroStain II, Leica) filtered through a 0.2 µm syringe filter, and then washed thoroughly with boiled water. A thin layer of carbon was then deposited onto each wafer to prevent charging during scanning electron microscopy.

Electron microscopy

WaferMapper software was then used with light-based wafer overview images to semi-automatically map the positions of all sections and relate them to fiducial markers. This enabled targeted section overview acquisition (758.8×758.8×60 nm³vx⁻¹ for 5.5 dpf; 741.5×741.5×50 nm³vx⁻¹ for 7 dpf). Semi-automated alignment of section overviews in WaferMapper then permitted targeting for imaging at higher resolutions²⁶.

Field emission scanning electron microscopy of back-scattered electrons was primarily conducted on a Zeiss Merlin equipped with a large-area imaging scan generator (Fibics) and stock back-scattered electron detector. An accelerating voltage of 5.0 kV and beam current of 7– 10 nA were used for most acquisition. Imaging of back-scattered electrons at the highest resolutions (4.0×4.0×60 nm³vx⁻¹) was performed on an FEI Magellan XHR 400L with an accelerating voltage of 5.0 kV and beam current of 1.6–3.2 A. Field of view sizes acquired from a given section varied depending on the cross-sectional area occupied by tissue. All acquisition was performed with a scan rate at or under 1 Mpx s⁻¹. For the 5 dpf larval zebrafish, this resulted in overhead-inclusive acquisition times of 5.4 days for section overviews (758.8×758.8×60 nm³vx⁻¹), 97 days for full transverse cross-sections (56.4×56.4×60 nm³vx⁻¹), and 100 days for high-resolution brain images (18.8×18.8×60 nm³vx⁻¹).

Continued development of faster ssEM technologies⁴⁴ will hasten the re-imaging process and permit whole-brain studies in a fraction of the time required here.

Image alignment and intensity normalization

Producing anatomically consistent image registration over ~18,000 sections required control of region of interest drift, over-fitting, magnification changes, and intensities. In order to quickly assess the quality of the dataset and begin reconstructions, we initially performed affine intra- and inter-section image registrations with Fiji⁴⁵ TrakEM2 alignment plug-ins⁴⁶. These results revealed that additional nonlinear registration was required in order to compensate for distortions likely caused by section compression during cutting and sample charging during imaging. While the state-of-the-art elastic registration method⁴⁷ also provided in Fiji⁴⁵ as a TrakEM2 alignment plug-in achieved excellent local registration, we experienced difficulty—at least without modification to the existing implementation—in achieving an anatomically consistent result that preserved the overall larval zebrafish structure, largely due to struggles with constraining region of interest drift across magnification changes and correcting for shearing caused by sectioning. We also determined that the similar AlignTK⁹ method, which uses Pearson correlation as the matching criterion coupled with spring mesh relaxation to stabilize the global volume, was likely to suffer from similar problems and would require substantial additional data handling to operate on our multi-resolution dataset.

Therefore, in order to preserve the overall larval zebrafish structure and simultaneously achieve high-quality local registration, we turned to a new Signal Whitening Fourier Transform Image Registration (SWiFT-IR) method^{43, 48}. Compared to conventional Pearson or phase correlation-based registration approaches, SWiFT-IR produces more robust image matching by using modulated Fourier transform amplitudes, adjusting its spatial frequency response during matching to maximize a signal-to-noise measure as its indicator of alignment quality. This alignment signal better handles variations in biological content and typical data distortions. Additionally, SWiFT-IR achieves higher precision in block matching as a result of the signal whitening, improves computational speeds with the computational complexity advantages of fast-Fourier transforms, and reduces iterative convergence from thousands to dozens of steps. Together, these capabilities enable a model-driven alignment in place of the usual approach of comparing and aligning a given section to a pre-selected number of adjacent sections.

The SWiFT-IR model we used consisted of an estimate of local aligned volume content formed by a windowed average, typically spanning ~6 µm along the axis orthogonal to the sectioning plane (z, anterior-posterior). Damaged regions, in particular partial sections, were removed from the model to avoid adversely influencing alignment results. This model then served as a registration template, where raw images were matched to the current model rather than nearby sections. Alignment proceeded in an iterative fashion starting at 758.8×758.8×60 nm³vx⁻¹ (section overviews) and progressing incrementally to 56.4×56.4×60 nm³vx⁻¹ for most regions outside the brain and 18.8×18.8×60 nm³vx⁻¹ for most regions inside.

At each resolution, source images were iteratively aligned to the current model until no further significant alignment gain could be achieved, as indicated by the SWiFT-IR signal-to-noise figure of merit. The model was then transferred to higher resolution data by applying the current warpings to source data for that scale. Iterative model refinement then continued at this subsequent level. Although most computations were locally affine, residual nonlinear deformations, particularly at the highest resolutions, were represented by a triangulation mesh that deformably mapped raw data onto the model volume.

Importantly, access to the lowest resolution section overview data for each section permitted us to build an initial model that constrained subsequent registration steps to the overall larval zebrafish structure. Although their resolution and signal quality were intentionally low in favour of rapid acquisition, the fact that overviews were quickly captured with the same microscope settings and included support tissue provided key constraints for model refinement that resulted in a more accurate global result.

More specifically, the 17,963-section overview image volume was processed using SWiFT-IR to produce an initial model at 564×564×600 nm³vx⁻¹. Although the lowest resolution section overview images were each captured at 758.8×758.8×60 nm³vx⁻¹, the relative oversampling orthogonal to the sectioning plane enabled a geometrically accurate model at 564×564×600 nm³vx⁻¹. This initial model was then cropped and warped using SWiFT-IR–driven matching across the midline axis to remove cutting compression, rotations, and other systematic variations in the specimen pose. The 16,000-section 56.4×56.4×60 nm³vx⁻¹ volume was next downsampled to 564×564×600 nm³vx⁻¹ and aligned to the initial overview model, resulting in an improved model. The matching and remodeling process was iterated at this scale until there was no further improvement in SWiFT-IR match quality. The final model at this scale was then expanded to 282×282×300 nm³vx⁻¹ and similarly aligned in an iterative fashion. This model volume (~6 Gvx, 1600×1400×2667 vx) was convenient for rapid viewing to identify and manually correct defects and refine the pose. Further scales at 169.2×169.2×180 nm³vx⁻¹ and 56.4×56.4×60 nm³vx⁻¹ were similarly processed by successively expanding the model and aligning until no significant improvement in the figure of merit was reached. The 12,546-section highest resolution 18.8×18.8×60 nm³vx⁻¹ image set was then registered using the final 56.4×56.4×60 nm³vx⁻¹ volume as its model.

Image intensity was next adjusted across sections to achieve a consistent background level to match the average over a tissue-free region defined by a 256×256 px² area. Many images were acquired at 16-bit depth and were converted in this process to 8-bit depth. The target background level was mapped to intensity 250, which left headroom for bright pixels while keeping tissue of interest from saturating. Next, a linear intensity fit between the background and a second level, typically the average grey level of a continuous trajectory region on the right side of the brain, was made to adjust the intensity values for each section.

Correspondence across light and ssEM datasets

Correspondence of individual neurons or functional reference atlas regions across imaging modalities was achieved with landmark-based 3-D thin-plate spline warping of each fluorescence dataset to the ssEM dataset using BigWarp⁴⁹.

For matching in vivo two-photon laser-scanning microscopy data from the same specimen, we primarily chose landmarks consisting of distinctive arrangements of low-fluorescence regions where GCaMP5G was excluded and could be easily matched to similar patterns of nuclei in the ssEM dataset. This process was difficult in regions with low fluorescence signal (Extended Data Fig. 7e), where many cells were packed closely together (Extended Data Fig. 7f), and at locations where new neurons were likely added between light microscopy and preparation for ssEM (Extended Data Fig. 7g). In the future, improving the light-level data with specific labelling of all nuclei and faster light-based imaging approaches should improve the ease and accuracy of matching neuron identity.

Two functional reference atlases with many separate labels were also registered to the ssEM dataset. For matching the Z-Brain atlas⁵⁰, we chose landmarks based on identifiable structures in the Z-Brain averaged elavl3:H2B-RFP or anti-tERK fluorescence image stacks that were also observed in the ssEM dataset. These structures primarily consisted of region boundaries, known clusters of neurons, midline points, ganglia, and the brain outline. The same Z-Brain landmarks were used for transforming a version⁵¹ of the Zebrafish Brain Browser⁵² that was previously registered into the Z-Brain atlas.

Image annotation and neuron reconstruction

Reconstructions across multi-resolution ssEM image volumes profits from being able to simultaneously access and view separate but co-registered datasets. Without this, some of the time benefits of our imaging approach would be offset by the need to register and track each structure across volumes that span both low-resolution, large fields of view and high-resolution, specific regions of interest. With this in mind, we added a feature to the Collaborative Annotation Toolkit for Massive Amounts of Image Data (CATMAID) neuronal circuit mapping software^{53, 54} to overlay and combine image stacks acquired with varying resolutions in a single viewer (Extended Data Fig. 6). This is made possible by rendering using WebGL. Additionally, this new feature combines stacks via a configurable overlay order, introduces blending operations for each overlaid stack, and enables programmatic shaders for dynamic image processing. When overlaid stacks resolutions differ, the nearest available zoom level for each stack is interpolated. Missing data regions can be omitted or rendered with interpolation. To account for the increase in data storage and bandwidth when viewing multiple image stacks, the CATMAID image data hierarchy was also extended with a shared graphics card memory cache of image tiles using a least-recently-used replacement policy. All additions and modifications to the CATMAID software are now incorporated into the main open-source release.

Manual reconstruction was conducted using our modified CATMAID by placing nodes near the centre of each neuronal structure on every section in which it could be clearly identified. This led to a wire-frame model for each annotated structure. Starting points for reconstruction (“seeds”) of myelinated processes were manually identified by searching for profiles surrounded by the characteristic thick, densely stained outline associated with staining of the myelin sheath⁵⁵ (Fig. 1e,g,i). The search protocol required viewing all tissue on a given section from the upper-left corner to lower-right at the highest available resolution. To obtain seeds for projectome reconstruction, searching was repeated every 50 sections throughout all 16,000 sections acquired at or higher than 56.4×56.4×60 nm³vx⁻¹. Many annotations were produced in an affine-only alignment space before being mapped into the final SWiFT-IR alignment space. The reconstructions reported here represent ~450 days of human annotation.

For visualization and reported length measurements, each mapped wire-frame was smoothed using custom python-based implementation of a Kalman smoothing algorithm on a space defined by manually annotated points within unique segments. The initial state variables for smoothing were derived by an optimization of point-to-line distance to connected reconstruction segments. Other variables were tuned with the Estimation Maximization algorithm of the pykalman library to compensate for a lack of human input where data was unavailable due to lost or partial sections. Because the final image alignment was of good quality, smoothing in this manner should produce a slight underestimate in reported reconstruction path lengths.

Neurons with known projection patterns or identities were named in the CATMAID database. For example, the reconstruction of a neuron innervating the right anterior macula (utricle) might be named as “Ear_AnteriorMacula_R_01”, while an identified neuron such as the left Mauthner neuron was named “Mauthner_L”. Two identifiable left-right neuron pairs belonged to the “MeM” class, which emanates from the nucleus of the medial longitudinal fasciculus (nucMLF). On each side, these were differentiated into dorsal (MeMd) and ventral (MeMv) subclasses based on consistent soma positioning.

Visualization

Image volumes and reconstructions were primarily visualized using Vivaldi⁵⁶, a domain-specific language for rendering and processing on distributed systems, because it provides access to the parallel computing power of multi-GPU systems with language syntax similar to python.

For volume visualizations, we used a direct volume rendering ray-casting technique in which an orthogonal or perspective iterator was marched along a viewing ray while sampled voxel colours were accumulated using an alpha compositing algorithm. We screened out regions containing only support tissue during rendering with labelled volumes constructed by interpolating between manually produced masks indicating which image voxels belong to each separate tissue region. In cases where separate image volumes of the same region were rendered together (e.g., ssEM and fluorescence combined), direct volume rendering was performed by combining front-to-back colour and alpha compositions formed from the different transfer function belonging to each image volume.

For volume visualizations including reconstructions, direct volume rendering of image data was combined with streamline rendering of reconstructions using two different techniques. The first combined an OpenGL framebuffer with the Vivaldi volume rendering. In this case, each streamline was rendered using OpenGL as a tube into an off-screen buffer (i.e., Framebuffer Object). Vivaldi then compared the resulting render and depth buffers to perform direct volume rendering of only the image data above the streamline depth value. This made it possible to ignore image voxels obscured by streamlines, which were treated as opaque. The second technique involved generating a complete streamline volume by 3-D rasterization. This streamline volume was then combined with the image volume for direct volume rendering. The former technique is faster and can cope with dynamic streamline changes, but the latter was found to yield better overall rendering quality. Visualizations of reconstructions without the image volume context were rendered either in the CATMAID 3-D WebGL viewer or plotted in MATLAB. When reconstructions are shown without specific labelling, colours were assigned randomly from a custom palette.

Reference plane (e.g., horizontal, sagittal, and section) images were rendered with Vivaldi by detecting the zero-crossing of each viewing ray and the plane. Support for viewing opaque data views in some spatial regions alongside the semi-transparent volume visualization views in other regions was introduced as a new Vivaldi function, clipping_plane. Similarly, contour (nonplanar) reslice support was developed to illustrate a flattened view along a specific reconstruction path consisting of vertical line segments extracted from the image volume.

For many cases, the size of the volume being rendered was larger than available memory. In order to support out-of-core processing, we developed and integrated into Vivaldi a slice-based streaming computing framework using the Hadoop distributed file system (HDFS) that will be reported elsewhere.

Symmetry analyses

Initial observations of apparent myelinated axon symmetry were found during visual inspection (Fig. 3; Supplementary Videos 8–9). To quantitatively assess the extent of symmetry, we developed a 3-D symmetry plane fitting method and two symmetry analyses: one that produces a cost associated with the 3-D shape and position similarity between reconstructed structures and another that compares the relative 2-D (cross-sectional) positioning of two identified neuron axons on one side with that for the contralateral axons with the same identities Only the longest reconstructed path from the soma through the myelinated axon projection was considered in plane fitting and symmetry analyses. Dendrites or short axonal branches were ignored. Each resulting reconstruction path (skeleton) was represented as an ordered list of nodes (points) taken directly from manual reconstructions. Sidedness (left or right) was determined by soma position.

The new 3-D symmetry plane fitting and 3-D symmetry comparison analyses approaches were described elsewhere⁵⁷. The symmetry plane fitting, in brief, involves choosing an approximate symmetry plane, reflecting the complete set of points belonging to the reconstruction subset of interest with respect to this plane, registering the original and reflected point clouds with an iterative closest point algorithm, and inferring the optimal symmetry plane from the reflection and registration mappings. The subset of reconstructions from which this plane fitting was performed consisted entirely of identified neurons whose axon projections formed part of the ~30µm-diameter medial longitudinal fasciculus (MLF), recognized with the help of refs. 29, 58.

The 3-D symmetry comparison for each template reconstruction on one side, in brief, involved reflecting all contralateral skeletons and computing a matching cost via dynamic time warping (DTW) between the template and each reflected skeleton. The reconstruction subset analysed in this fashion was restricted to identified neuron classes with 1-2 members per side whose axons formed part of the MLF. For our purposes, the DTW cost was taken as the sum of the Euclidian distances between all matched points normalized to the number of matched point pairs (Extended Data Fig. 9a–c). The DTW gap cost parameter for matching a point in one sequence with a gap in another was set to zero because our data was sampled at a nearly constant rate and we sought the optimal subsequence match even in cases where one is shorter than or offset with respect to the other. To compensate for unmatched regions (i.e., overhangs), the DTW cost was then multiplied by a penalty factor proportional to the sequences lengths remaining unmatched (total length divided by matched length). Comparing each reconstruction on one side to all reconstructions from the opposite side formed a cost matrix (Fig. 4b) from which an optimal pairwise assignment could be determined without any bias introduced from the previously determined identities. The Munkres algorithm⁵⁹ was then used to compute a globally optimal pairwise assignment.

We next sought to compare the relative 2-D positioning for each set of two axons on one side with the contralateral set having the same identities. The reconstruction subset analysed in this fashion was restricted either to the Mauthner cell and nucMLF neurons (Fig. 4e–g) or the larger set of 44 identified reticulospinal neurons (Extended Data Fig. 9i–k).

To start, we compensated for a small angle offset in the sectioning plane relative to the true transverse plane by projecting the point coordinates of reconstructions such that the previously computed symmetry plane became the plane x = 0. Given the transverse planes Z₀, Z₁, and a projected skeleton S containing points s = (s_x, s_y, s_z), we let S_Z0,Z1 = {s ∈ S : Z₀ ≤ s_Z < Z₁}. That is, S_Z0,Z1 was taken as the subset of points from S whose coordinates s_Z are contained in the interval [Z₀, Z₁). We refer to the subset of ℝ³ bounded by Z ∈ [Z₀, Z₁) as the slice [Z₀, Z₁). For each slice [Z₀, Z₁) and skeleton S, we defined 〈S_{Z0, Z1}〉 as the mean of the elements in S_{Z0, Z1}. This mean was then taken as representative of the skeleton S in slice [Z₀, Z₁) for analysis and plotting. Note that all analysis and plotting presented in static form was based on a slice thickness corresponding to a single section (~60 nm), where each slice consisted simply of adjacent sections. Larger slice sizes were used for dynamic presentation (Supplementary Video 10) in order to reduce video duration and size.

For comparing a set of two axons with its contralateral counterpart, we then took s₁, …, s_n to be the set of representative points in a fixed slice for skeletons S₁, …, S_n and took t₁, …, t_n to be the representative points (for the same slice) of the respective skeletons T₁, …, T_n that were previously matched to S₁, …, S_n by the Munkres algorithm assignment after 3-D symmetry analysis.

To quantify the degree of similarity, we devised two measures (Extended Data Fig. 9e). The first, termed the angle difference, a_i,j, between a set of two axons and their contralateral counterparts, was defined as:

The second, termed the distance difference, d_i,j, between a set of two axons and their contralateral counterparts, was defined as:

where i, j were skeleton indices, $\{t_i^s\}$ was the reflections of {t_i} with respect to the computed plane of symmetry, and M is the maximum of $|\Vert s_j-s_i\Vert -\Vert t_j^s-t_i^s\Vert |$ across all axon sets and all slices. Note that a_i,j and d_i,j were normalized such that they could vary from 0 (no difference, 0° or 0 µm) to 1 (maximum difference, 180° or ~8 µm). Further, when the points s_i and s_j were perfectly symmetric with respect to points t_i and t_j, then a_i,j = 0 and d_i,j = 0.

To visualize this quantification, a difference matrix, D, was generated for each slice such that D(i,j) = a_i,j if j > i and D(i,j) = d_i,j if j < i (Fig. 4f; Extended Data Fig. 9f; Supplementary Video 10). Calculating the variance for each element in D across all slices showed which axon sets deviated most with respect to the reflection of their contralateral counterparts (Fig. 4k). Heatmaps of the vectorised upper (j > i) and lower (j < i) triangles of D across slices additionally revealed locations with differences between axon sets and their contralateral counterparts (Fig. 4h–i; Extended Data Fig. 9f–h; Supplementary Video 10). Plotting the sum of all a_i,j and d_i,j values for a given slice further illustrated where differences were present across the subset (Fig. 4j). Finally, the same analysis was performed after artificially swapping the identities (assignment) of the two axon reconstructions with the lowest 3-D symmetry analysis costs (MeLc and Mauthner) to provide a basis for comparison (Fig. 4j).

Code availability

Custom software tools generated for data handling, visualization, and analysis are publicly available (http://zebrafish.link/hildebrand16/code). Our modifications to CATMAID^{53, 54} software are included in the main open-source release (http://github.com/catmaid/catmaid). More information on SWiFT-IR alignment software is publicly available (http://www.mmbios.org/swift-ir-home).

Data availability

All aligned ssEM data, reconstructions, transformed functional reference atlases, and an introductory guide are publicly available (http://zebrafish.link/hildebrand16). Image data is served as a collection of 8-bit 1024×1024 px² PNG images with an optional tRNS value of 255 specified to enable transparency. The original resolution for each image stack was down-sampled multiple times to create a resolution hierarchy that provides a smooth visualization experience, where each level in the hierarchy corresponds to an image that is half the size as in the previous level. The entire aligned image dataset requires ~2.7 TB of disk space as compressed PNG images (607 GB for 56.4×56.4×60 nm³vx⁻¹ ssEM data, 1,824 GB for 18.8×18.8×60 nm³vx⁻¹ ssEM data, 355 GB for 4.0×4.0×60 nm³vx⁻¹ ssEM of dorsal neuromasts, 1 GB for 600×600×1200 nm³vx⁻¹ Z-Brain data, and 3 GB for 600×600×1200 nm³vx⁻¹ Zebrafish Brain Browser data). Data and reconstructions are served to end users via Amazon Web Services (AWS), with an instance of our modified CATMAID^{53, 54} software deployed on the Elastic Compute Cloud (EC2) that points to static images hosted by the Simple Storage Service (S3) built-in web server.