Subvoxel light-sheet microscopy for high-resolution high-throughput volumetric imaging of large biomedical specimens

Abstract. A key challenge when imaging whole biomedical specimens is how to quickly obtain massive cellular information over a large field of view (FOV). We report a subvoxel light-sheet microscopy (SLSM) method enabling high-throughput volumetric imaging of mesoscale specimens at cellular resolution. A nonaxial, continuous scanning strategy is developed to rapidly acquire a stack of large-FOV images with three-dimensional (3-D) nanoscale shifts encoded. Then, by adopting a subvoxel-resolving procedure, the SLSM method models these low-resolution, cross-correlated images in the spatial domain and can iteratively recover a 3-D image with improved resolution throughout the sample. This technique can surpass the optical limit of a conventional light-sheet microscope by more than three times, with high acquisition speeds of gigavoxels per minute. By fast reconstruction of 3-D cultured cells, intact organs, and live embryos, SLSM method presents a convenient way to circumvent the trade-off between mapping large-scale tissue (>100  mm3) and observing single cell (∼1-μm resolution). It also eliminates the need of complicated mechanical stitching or modulated illumination, using a simple light-sheet setup and fast graphics processing unit-based computation to achieve high-throughput, high-resolution 3-D microscopy, which could be tailored for a wide range of biomedical applications in pathology, histology, neuroscience, etc.


Introduction
In optical microscopy, high-resolution (HR) volumetric imaging of thick biological specimens is highly desirable for many biomedical applications, such as development biology, tissue pathology, digital histology, and neuroscience. To obtain information on cellular events from the larger organism, e.g., a live embryo, intact tissue, or an organ, spatiotemporal patterns from the micro-to mesoscale must be in toto determined and analyzed. [1][2][3][4][5] Thus, there is a growing need to develop HR, highthroughput imaging methods that can map entire large-volume specimens at high-spatiotemporal resolution. 6,7 Recently, lightsheet microscopy (LSM) has emerged as a technique of choice that can image samples with low phototoxicity and at high speed. [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] However, similar to conventional epifluorescence methods, LSM remains subject to the fundamental trade-off between high illumination/detection numerical apertures (NAs) and wide imaging fields of view (FOVs). In addition, an accurate digital sampling by the camera is also compromised by the need for large pixel size with high-fluorescence sensitivity. Therefore, the achievable resolution of current LSM systems is often pixel-limited under large FOVs, yielding inadequate optical throughput for digital imaging of mesoscale organisms at the cellular resolution. Tile imaging-based LSM systems have been developed to artificially increase the space-bandwidth product (SBP), 24 hence, realizing HR imaging of large specimens. 18,[25][26][27][28][29] Despite the compromised speed induced by repetitive mechanical stitching, the high illumination/detection NA configuration in tile imaging induces increased phototoxicity for increasing sample size and limits fluorescence extraction from deep tissue. In addition, several techniques, such as Fourier ptychographic microscopy, 30,31 synthetic aperture microscopy, [32][33][34][35] contact-imaging microscopy, 36,37 wavelength scanning microscopy, 38 and lens-free digital holography, [39][40][41] have recently provided a computational means of reconstructing a wide-FOV, HR image based on a number of low-resolution (LR) frames having certain correlations in the space, frequency, or spectrum domain. [42][43][44] However, the majority of these methods target twodimensional (2-D) bright-field microscopy and are not compatible with volumetric fluorescence imaging of thick samples.
Here, we present an imaging method capable of providing three-dimensional (3-D) super-resolution and increased SBP for conventional LSM, without the involvement of mechanical stitching or complicated illumination modulation. This method, termed "subvoxel light-sheet microscopy (SLSM)," shares its roots with pixel super-resolution (PSR) techniques 36,37,[41][42][43]45 and works by efficiently computing a number of LR, undersampled, and shift-modulated LSM volumes in the spatial domain to reconstruct an output with significantly higher resolution throughout the entire sample. Unlike either PSR or conventional LSM, our SLSM method super-resolves large-scale fluorescence images in terms of voxels for the first time, expanding the system SBP, which was originally limited by the low-magnification/NA optics. Furthermore, in SLSM, these special image sequences carrying high-frequency spatial information can be obtained through a simple retrofit of the z-scan apparatus employed in conventional LSM systems. Moreover, segmentation of the raw image sequence, modeling of the spatial shifts beyond the system resolution, and reliable estimation of the HR output are all based on a graphics processing unit (GPU)-accelerated parallel computation flow. Therefore, rapid and convenient generation of raw image data, in conjunction with efficient image modeling and reconstruction, realizes highthroughput, HR imaging of large biological specimens via the SLSM technique.
In this study, this SLSM capability is broadly verified by imaging various samples, such as 3-D cultured cells, a live zebrafish embryo, and intact mouse organs. In addition, we demonstrate that SLSM can be combined with multiview data fusion 46 to achieve complete imaging of scattering samples, realizing an isotropic resolution of ∼1.6 μm (compared with the original ∼6.5 and 26 μm) throughout a volume of >100 mm 3 . In the following, we elucidate the SLSM experimental setup, discuss implementation of the subvoxel-resolving (SVR) computation, and demonstrate applications to wholeorganism imaging.

Off-Detection-Axis Scanning Setup
In this study, SLSM imaging was implemented based on a home-built selective plane illumination microscope (SPIM), which is a well-known LSM modality (see Methods and Supplementary Fig. S1 for the full optical layout). This approach functions under a low-magnification setup with a wide-FOV covering the entire specimen. Unlike a regular SPIM with z-scan parts, the SLSM additionally contains a customized tilting plate that can adjust the scanning axis, allowing sample scanning in a direction with a certain deviation angle θ relative to the z axis (detection axis) [ Fig. 1(a)]. As the sample is continuously moved across the laser sheet in this nondetection-axis direction, a camera records images of the sequentially illuminated planes at a high acquisition speed [ Fig. 1(b)]. By matching the frame rate r and scanning velocity v, a fine step size s, typically hundreds of nanometres (s ¼ v∕r), is generated between adjacent frames [ Fig. 1(b)]. Thus, every nonaxial incremental step provides unit subvoxel shift components s x , s y , and s z simultaneously in both the lateral and axial directions, through a one-directional scan [ Fig. 1(b)]. In addition, this continuous scanning mode matched by fast camera acquisition provides SLSM with a high imaging throughput of up to hundreds of megavoxel per second. Depending on the sample size and fluorescence intensity, SLSM acquisition usually requires only a few seconds to a few minutes, finally generating a raw image stack containing thousands of spatially modulated frames.

Principles of Subvoxel-Resolving Procedure
An SVR algorithm was designed to specifically segment the acquired shift-encoded, wide-view, LR image stack (denoted by P) into multiple substacks, model them as probability distributions with subtle spatial correlations, and reconstruct a final output image that encompasses significantly improved resolution over the entire large-volume sample [ Fig. 1(c)]. In this approach, the raw image stack P is first divided into a number of LR 3-D images P k (k ¼ 0,1; 2; …; n). In each LR P k , the voxel width wðx; yÞ is simply given by the ratio of the camera pitch size and magnification factor. The voxel depth d, denoting the axial spacing implemented when extracting the slices from the raw stack, is set to approximately one-third of the laser-sheet longitudinal extent lð1∕e 2 Þ to satisfy the Nyquist sampling principle. The segmentation is, therefore, performed by reslicing the raw image every l∕3s z frames as shown in Figs. 1(b) and 1(c) (step 1). The number of total LR images n, which implies the amount of subvoxel information extractable from the raw image, is determined by dividing one LR w by s x ; this indicates that these n images are correlated with each other in terms of the SVR displacements [Figs. 1(b) and 1(c)]. Each segmented LR P k can be considered as a blurred and undersampled "regular SPIM" image presenting standard resolutions determined by the system optics. Meanwhile, all the P k can be spatially registered to the HR image I to be solved (also the reference image P 1 ), according to a nonaxial, SVR shift s k ¼ ðk − 1Þ Ã s [ Fig. 1(c): red-, green-, and purple-bordered images]. The SVR procedure then reconstructs the HR I based on computation of these LR images modeled with their known subvoxel shifts, the system optical blurring, and the camera discretization.
The design of the SVR procedure follows a straightforward principle similar to PSR techniques functioning in the 2-D Fig. 1 Principles of SLSM. (a) SLSM geometry with wide laser-sheet illumination and large-FOV detection. The sample is scanned in one direction (red arrow) with the deviation angle θ being 10 deg to 20 deg relative to the z axis. (b) Sequence of spatially modulated, LR images recorded under high-frame rate. The step size s is much smaller than the laser-sheet longitudinal extent lð1∕e 2 Þ. The raw image stack is split into multiple substacks (indicated by different colors), with a voxel depth d that is one-third the laser-sheet thickness l. Each LR substack P k is correlated to the first reference stack P 1 (red), with a nonaxial, subvoxel shift s k ¼ ðk − 1Þ Ã s. For instance, P 2 (yellow) can be registered to P 1 with unit displacement s, which consists of 3-D components s x , s y , and s z simultaneously. (c) Block diagram of iterative SVR procedure. The spatially correlated, wide-view, LR stacks (step 1) are input to an MLE (step 2) to iteratively reconstruct a voxel-superresolved image. Block G k represents a gradient backprojection operator that compares the k 'th LR image to the estimate of the HR image in the m'th steepest descent iteration. At the final step (3), voxel realignment is applied to recover the sample shape accurately from the slight deformation caused by the tilted scan. spatial domain. 36,42 It seeks the super-resolved estimate I that is most consistent with multiple measurements P k after a series of degradation operators that reasonably model the digital imaging process being successively applied. 36,43 Here, I can be specifically solved by minimizing the following cost function: (1) where ρ is the difference between the model and measurements, S k is the geometric motion operator between the HR estimate I and the k'th LR P k , the point spread function (PSF) of the SLSM system is modeled by the blur operator O, and D k is the decimation operator that models the camera digital sampling. Theoretically, the computation estimates an HR image having maximum likelihood with the LR inputs after given degradations S k , O k , and D k are applied. Figure 1(c) (step 2) shows an iterative maximum-likelihood estimation (MLE) solution procedure for I. The program first generates an initial guess of HR IðI 1 Þ, which is simply the interpolation of P 1 . Each P k is then compared with the warped, blurred, and decimated I 1 using the gradient back-projection operator (block G k ). Their differences are summed and thereafter weighted by a factor β, for the calculation of the second estimate. This process is iterated until a converged solution of the m'th I is approached after feeding the (m − 1)'th estimate. A voxel realignment is applied in the final step to recover an accurate reconstruction from the slight deformation caused by the nonaxial scan

Imaging Characterization
Simulation results of the proposed SVR procedure are provided in Supplementary Fig. 2(a). A sequence of images from consecutively illuminated planes was rapidly recorded for 1 min under a 100-frames/s acquisition rate, targeting a specific volume-of-interest containing dense cell spheroids [ Fig. 2(a)]. Note that each LR image subdivided from the raw stack simply adopts the limited resolution from the system optics; hence, only the rough shape and distribution of the selected spheroid were identified, whereas the single cells remained unresolvable [ Fig. 2(b1)]. The SVR procedure then began with an initial guess, which was simply a 4× interpolation of the first reference image [ Fig. 2 SLSM provides a wide-FOV of ∼23 mm 2 from its low-magnification DO (4×/0.13) and low NA PI (0.022), whereas its achieved resolutions are similar to those of 20×-SPIM (20×/ 0.45-DO/0.07-NA PI). Thus, the SLSM can be regarded as a light-sheet microscope combining the FOV advantage of 4×-SPIM with the resolution advantage of 20×-SPIM. From another perspective, its stitching-free, continuous scanning mode exhibits a much higher acquisition throughput as well as lower photobleaching than the stitching results of the 10×-SPIM, 20×-SPIM, and confocal microscope (Supplementary Notes).
We calculated the imaging resolution, speed, and photobleaching of 4×-SPIM, 10×-SPIM tile imaging, 20×-SPIM tile imaging, 20×-confocal microscopy, and 4×-SLSM. As shown in Fig. 2(f), the SLSM exhibited the highest effective throughput at ∼25 megavoxel SBP per second, which is >10 times higher than those of the other modalities. In addition to the increased SBP (FOV divided by resolution), the SVR computation, to some weak fluorescing extent, improves the signal-to-noise ratio after multiple image fusion ( Supplementary Fig. S7). Furthermore, circumventing the use of high-NA, high-maintenance optics renders SLSM considerably less vulnerable to chromatic aberration occurring in multicolor illumination 47 (Supplementary Notes and Supplementary Fig. S2), and resistant to the spherical aberration 48,49 that causes severe image deterioration for deep tissue (Supplementary Fig. S3). In the applications described below, this underlying robustness allows the SLSM prototype to image thick specimens at high spatial-temporal performance while retaining a relatively simple setup.

Multicolor Three-Dimensional Whole-Organism
Imaging at High Throughput SLSM does not require special engineering of the fluorescence emission to generate image modulation, rendering it compatible with various labeling techniques. In Fig. 3, we demonstrate SLSM imaging of two types of specimen: an optically cleared intact mouse heart (neonate, D2) exhibiting endogenous autofluorescence in cardiomyocytes [ Fig. 3(a)] and a two-color transgenic zebrafish embryo [3-day postfertilization (dpf)] tagged with green and red fluorescence proteins at the motor neurons (Islet1-GFP) and somite fast muscles (mlcr-DsRed), respectively [ Fig. 3(b)]. Furthermore, a time-course study on an    anesthetized live fish embryo from 48-to 72-h postfertilization (hpf) was implemented in this work to investigate the neuron/ muscle development in the hindbrain region as shown in Supplementary Fig. S9. The raw data acquisition was very fast, at a rate above 200 megavoxels per second. Meanwhile, to match the high-speed acquisition, we developed a GPU-based parallel computation flow to greatly accelerate the SVR procedure simultaneously, super-resolving the large-scale data at a high-processing throughput above 100 megavoxels per second. Hence, SLSM can rapidly image (acquisition + computation) these millimetre-size organisms on a scale of dozens of gigavoxel, providing super-resolved cellular structure-function information, such as myocardium architectures and interactions of developing motor neurons with somite muscles, within a few minutes to an hour. In Supplementary  Figs. S10 and S11, broader demonstrations of SLSM applications are provided, for whole-organism imaging of a Tg (cmlc2:GFP) adult zebrafish heart (60 dpf) and an αMHC Cre ; R26 VT2∕GK murine heart (three colors, P1), which are spatially or temporally more challenging to regular light microscopes.

Multiview SLSM for Isotropic Super-Resolution Imaging
Despite the use of chemical clearing, light scattering from deep tissue continues to pose a challenge for optical microscopy. Both laser excitation and fluorescence emission experience deflection and attenuation, which together deteriorate signals severely. In addition, even light-sheet setup remains anisotropic for the imaging of mesoscale samples, showing suboptimal axial resolution for several applications, such as cell phenotyping and neuronal tracing. 18 The multiview fusion method 46,47 was previously developed to address these problems. It functions by registering, weighting, and fusing a number of image stacks recorded under different views and finally recovers a stack that shows complete signals with improved axial resolution. [46][47][48][49][50] Here, we demonstrate that the SVR procedure can be combined with a multiview approach to in toto image thick and scattering samples at near isotropically improved resolution. We developed a four-view SLSM prototype by imaging a human umbilical vein endothelial cell (HUVEC) and human dermal fibroblast (HDF) 3-D sprouting network (Sec. 4). The cells were cocultured in a fibrinogen-based hydrogel mixed with dextrancoated Cytodex 3 microbeads (SiO 2 ), finally forming a complex and light-scattering cellular network containing fibrinogen, HUVEC beads, and HDF cells.
First, the sample was rotated by 90 deg for each set of acquired 2×-SLSM data (2×/0.06 detection, 0.015-NA sheet illumination, 580-nm step size, 7000 frames in 140 s), with four views of the raw stacks being obtained in total. Each set was subvoxel-resolved separately based on 32 groups of segmented LR inputs, generating four views of the anisotropic SVR image (4 × 4 × 2 enhancement). Multiview registration followed by weighted fusion was included in the final step to produce an output image with isotropic SVR. This multiview SVR (mv-SVR) workflow is shown in Fig. 4(a).
In Fig. 4(b), we compare the reconstructed volume renderings of a 0-deg SPIM image (1.1 gigavoxels, voxel size: 3.25 μm × 3.25 μm × 9 μm) and four-view SLSM (190 gigavoxels, reconstructed voxel size: 0.81 μm × 0.81 μm × 0.81 ). Because of the strong scattering of the fibrinogen-microbead substrate, the single-view SPIM visualizes an incomplete structure with poor resolution [ Fig. 4(b), left]. In contrast, the four-view SLSM reconstructs the entire sample, showing details of HUVEC sprouts from two adjacent dextran beads [ Fig. 4(b), right]. HR local views containing dense HDF fibers are shown in Figs. 4(c) and 4(d), also compared with single-view SPIM results. In addition to the isotropically improved resolution [resolvable distance: <1.6 μm; Fig. 4(d), inset] that enables clear identification of single sprouting HDF cells, the multiview SLSM (mv-SLSM) also restores the highly light-scattered area in the reference view [ Fig. 4(d), left]. Note that multiview subvoxel computation relies on the image registration and weighted fusion in a Fourier space and, therefore, the processing speed of four-view-fused SVR is ∼1 order lower than the computation for a single-view, reconstructing 18 megavoxels per second on average. Finally, by creating an isotropic, HR, panoramic visualization that encompasses 190 gigavoxels across a volume exceeding 100 mm 3 (total processing time: ∼3 h), the mv-SVR procedure can accurately analyze vast numbers of cells, such as sprouting branches and tip filopodia, over a piece of mesoscale tissue. Hence, mv-SVR provides a solid foundation for the exploration of mature vessel formation and vessel anastomosis processes, both of which are crucial for studying tissue regeneration.

Fast Quantitative Mapping of Mouse-Brain Neuronal Networks
We reconstructed a P30 Thy1-GFP-M mouse brain using an eight-view SLSM [ Fig. 5(a)]. HR volumetric renderings of the (i) hippocampus, (ii) thalamus, and (iii) cortex regions are shown in Figs. 5(b)-5(d), respectively. Compared with conventional SPIM with suboptimal quality, the eight-view SLSM imaging resolved fine neuronal substructures such as dendrites and axons at an isotropic resolution of ∼1 μm (350-gigavoxel SBP). Benefitting from remarkably improved visualization, neuronal subsets such as clustered astrocyte neurons in the thalamus region could be segmented with clearly separated nerve fibers [ Fig. 5(e), Imaris software]. Three long-projection neurons were also successfully identified and registered in the P30 reference mouse brain. The pathways of these projection neurons were subsequently annotated according to the standard mice brain atlas 51 as shown in Figs. 5e(ii)-5e(iv). The HR, highthroughput features of SLSM enable fast and accurate tracing of both long-distance projections and pathways within the neuronal circuits in the brain of a fluorescent protein transgenic mouse.

Discussion
SLSM is based on an unconventional off-axis scanning together with subvoxel reconstruction, so as to computationally surpass the resolution limit of a regular light-sheet microscope. This technique can be applied to most existing light-sheet microscopes by retrofitting with a readily available tilting stage, thereby expanding the optical throughput for fast HR mapping of large biomedical specimens. In principle, the SVR procedure estimates a 3-D image that best suits a certain conditional probability in the spatial domain. Provided the aperture function and subvoxel motion are accurately characterized, the maximumlikelihood link between the actual sample profile and recorded data allows the SVR to iteratively render an HR image using a limited-NA, undersampled configuration, which was originally incapable of providing such a small PSF (Fig. 1). Unlike most super-resolution fluorescence microscopy methods, which image a single cell or a few cells beyond the diffraction limit through processing of multiple frames acquired by patterned illumination or stochastic activation of molecules under highly specialized optics, SLSM is free from either illumination modulation or particular fluorescence labeling, being designed to rapidly enhance the inadequate accuracy obtained by unraveling large organisms under an ordinary small-NA/ large-view configuration. Its stitching-free, high-speed image acquisition followed by a parallelized GPU processing flow provides super-resolved 3-D visualization at a quasireal-time throughput. In this work, by imaging a variety of biological specimens from 3-D cells to developing embryos and intact organs, our SLSM robustly exhibited improved performance at low hardware and time cost, using efficient computation to circumvent the trade-off between larger imaging volume and more resolvable detail. The ability to rapidly accomplish cellular imaging of mesoscale organisms at hundreds-of-gigavoxel SBP renders SLSM a valuable tool for wide biomedical applications, such as phenotype screening, embryogenesis and brain connectome in histology, development, and neuroscience research, for which both large-scale statistics and HR details are highly desired. In addition to a single-view mode, SLSM can also be expanded by combination with multiview fusion. Through a rational balance between higher throughput and increased views, mv-SLSM can image thick and scattering samples with isotropic super-resolution and at moderately high-throughput. More broadly speaking, this method is currently optimized for LSM, which has a relatively high-image contrast and fast acquisition rate. However, the nonaxial scan as well as the SVR computation may be equally suited to other 3-D microscopy methods, such as confocal microscopy and 3-D deconvolution microscopy. Furthermore, we believe this subvoxel imaging strategy could prove to be transformative as it provides a general iterative resolution recovery method that could potentially be employed to improve other 3-D imaging modalities, which are possibly limited by inadequate sampling and poor focusing capability.

Experimental Setup
A four-wavelength, fiber-coupled semiconductor laser (CNI Laser, RGB 637/532/488/405, China) was used as an excitation source. The laser was first transformed into a collimated Gaussian beam with diameter of ∼8 mm (1∕e 2 value). Then, an adjustable mechanical slit (0-to 8-mm aperture) was used to truncate the beam in the horizontal direction and thereby tune the aperture of PI. The illuminating cylindrical lens (focal length f ¼ 40 mm) finally formed a wide laser sheet that optically sectioned large specimens with an NA adjustable from 0 to 0.1. Then, a 2× or 4× (Nikon Plan Apo Fluor 2×/0.06 or 4×/ 0.13 objective) infinity-corrected, wide-field detection path was constructed orthogonally to the PI path to collect the fluorescent signals. A four-degree-of-freedom motorized stage (x; y; z translation and rotation around the y axis, Thorlabs) integrated with a pair of customized tilting plates was constructed for sample mounting and scanning across the laser sheet in an off-detection-axis direction (Supplementary Fig. S1). For a given imaging configuration, the choice of θ was balanced between generating a small lateral shift component for a larger enhancement factor and scanning a short axial distance when completing an entire voxel shift. Usually, the diagonal of the LR voxel is chosen as the nonaxial scanning axis, with θ being 10 deg to 15 deg. A scientific complementary metal-oxide (e) Analysis of neural structures and connectivity. Three clustered astrocyte neurons in the thalamus region are segmented with different colors encoded in (i). Three registered projection neurons across the brain block are annotated according to the mice brain atlas (ii, iii, and iv). Abbreviations: cg, cingulum; ec, external capsule; dhc, dorsal hippocampal commissure; CPu, caudate putamen; ic, internal capsule; Rt, reticular thalamic nucleus; Prs, presubiculum; RSG, retrosplenial granular cortex; S1HL, primary somatosensory cortex hindlimb; RSA, retrosplenial agranular cortex; and scale bar: 500 μm.
semiconductor camera (Hamamatsu Orca Flash 4.0 v2 or Andor Zyla 5.5, pitch size: 6.5 μm) continuously recorded the planar images from the consecutively illuminated planes at a high speed of up to 200 frames/s.

SPIM, SLSM, and Multiview SLSM Acquisition
For the conventional SPIM imaging, the sample was scanned along the z axis with 2-to 9-μm step size depending on the different illumination NAs. To avoid motion blur generated by these incremental moves of several microns, it was necessary to implement the z-scan in a step-by-step manner by synchronizing the laser, motor, and camera using LabVIEW software (version 2014, National Instrument). This stepwise z-scan limited the acquisition speed to a maximum rate of 13 frames/s in our system (Thorlabs stepper ZST225B). Further, 10×-and 20×-tile SPIM imaging was automatically implemented in a z − x − y raster scanning method, and the acquired mosaic volumes were subsequently stitched together using the grid/ collection stitching plugin (ImageJ Software). Each time, only a small region of the entire illuminated plane was imaged; thus, the tile imaging was subjected to a photobleaching rate higher than those for the SPIM and SLSM.
For SLSM imaging, the camera was synchronized with a continuous sample scan under LabVIEW control, recording images in a sequential mode. The step size could, thus, be determined by matching the scanning velocity (10 to 30 μm∕s) and camera frame rate (20 to 200 frames/s). Depending on the sample dimension, tilt angle, and enhancement factor, this value varied from 140 to 600 nm. Compared with the voxel size of several microns in the segmented LR volumes, the motion blur caused by the continuous scan was negligible under such a small step size and did not affect the SVR computation accuracy. At the same time, the SLSM image acquisition could be fast under continuous mode, taking 30 to 300 s to obtain a sequence containing thousands to tens of thousands of frames. For multiview imaging, the stepper rotated the thick and scattering samples four to eight times, acquiring a number of raw SLSM image stacks under different views. The data stream was transferred from the camera to a redundant arrays of independent drives 0 volume of solid-state drives (4× crucial M550 1TB) in real time. The raw images were finally saved in a 16-bit TIFF format on the gigavoxel to tens of gigavoxel scale.

Efficient GPU-Accelerated SVR Reconstruction
In practice, a steepest descent method was used in the SVR computation to iteratively approach a converged super-resolved solution at high efficiency. Herê where S T k ; O T k ; D T k represent the inverse operations of S k ; O k ; D k , respectively. A MATLAB script was first developed as a proof of principle for the SVR (Supplementary Notes and Supplementary Fig. S4). Using a desktop workstation with dual Intel E5-2630 v4 central processing units (no GPU involved), the processing throughput did not exceed 0.3 megavoxels per second, almost 3 orders slower than the image acquisition.
Reconstructing data with volumes of tens to hundreds of gigavoxels for large specimens would require an impractically long time, i.e., multiple days. Therefore, we developed a GPU-based parallel workflow to accelerate the processing, achieving near 3-order acceleration (Supplementary Table 1). Depending on the degree of parallelization and the number of views provided, the processing throughput varied from tens to hundreds of megavoxels per second, depending on the power of the dual NVidia TESLA P100 graphical cards. For example, in a singleview configuration, the program explicitly finished an HR reconstruction of an intact heart in ∼12 min (90 gigavoxels) and a dual-color reconstruction of an entire zebrafish embryo in ∼4 min (34 gigavoxels). It should be noted that this speed could be further increased by employing additional GPUs.
For mv-SVR of thick and scattering HUVEC-HDF sprouting, 32 groups of LR images were extracted from the raw sequence of each view. The unit lateral and axial shifts were 100 and 558 nm, respectively. Considering the optical blurring and camera undersampling, the effective lateral and axial resolutions of each LR image were ∼7 and 26 μm, respectively, yielding ∼1.1-gigavoxel SBP (voxel size: 3.25 μm × 3.25 μm × 9 μm) over ∼100-mm 3 volume.
Applying the aforementioned SVR procedure, the superresolved image I single obtained under each view could be reconstructed as intermediate results. In each reconstructed I single , an increased SBP of 34 gigavoxels was presented anisotropically (0.81-μm × 0.81-μm × 4.5-μm voxel spacing). According to the reported multiview fusion method, we interpolated all I single with isotropic 0.81-μm pitch size, rotated the second view, and registered it to the 0-deg reference view by iteratively matching their histograms. An initial fusion of the images was generated by taking the weighted average of the two registered views in a Fourier space: where _ _ I is the Fourier transform of I and the w terms are the corresponding weights that reasonably average the two SVR views. Then, this initial fusion was used for the second registration and fusion iteration, in which the reference became the fusion of the first stage. By repeating this process, we obtained the fused image I merged that accumulated the effective information from all the SVR views. 46 This procedure can be expressed as follows: _ _ I ð1;…;jÞ;merged ¼ w ð1;…;j−1Þ · _ _ I ð1;…;j−1Þ;merged þ w j · _ _ I j;registered ; where j ¼ 3; 4 in a four-view SVR fusion. The weighting of the average was determined by the expected signal-to-noise ratio as follows: A multiview deconvolution was applied in the last step to obtain the deblurred output, which exhibited a final SBP of ∼190 gigavoxels with isotropically improved resolution. Implementation of the SVR and mv-SVR methods is also detailed in the Supplementary Notes. At present, as a proof-ofconcept, the mv-SVR is based on histogram registration and weighted fusion, and has not been performance optimized. Note that more efficient bead-based registration 50 and Bayesian deconvolution 52 have recently been reported. We believe that SVR can be combined with these techniques to yield even better image quality and higher speed. Further, the capacity was verified using four-view data in this work; however, it is quite certain that the result can be further optimized through provision of additional views. 46

Image Visualization and Analysis
The visualizations of SPIM, SLSM, and confocal data were performed using Amira (Visage Imaging). Planar images were presented in their original formats, unless otherwise mentioned, with no sharpening, interpolation, or registration applied. Maximum intensity projections and volume renderings were performed using the ProjectionView and Voltex functions in Amira with built-in colormaps. Long-projection neuronal tracing of the clarified mouse brain was performed using the Filament module in Imaris (Bitplane). See the supplementary figures for additional details and spatiotemporal visualizations of the developing zebrafish embryo, adult zebrafish fish heart, and neonate mouse heart.

Preparation of Large-Scale, 3-D Cultured Cells
NHBE cells (Lonza, Walkersville, Maryland) were prepared. The growth factor-reduced reconstituted basement membrane known as Matrigel (BD Biosciences, Bedford, Massachusetts) was used for the 3-D culture experiments. The cells were seeded into 300 μL of Matrigel in a liquid state and gently mixed with a pipette. The cells were then placed in a 24-well plate and incubated for 25 min at 37°C to gelatinize. Next, 1 mL of bronchial epithelial growth medium (Lonza, Walkersville, Maryland) was added to the 24-well plate and the medium was exchanged every other day.

Fluorescent Staining
Direct and indirect immunofluorescent stainings were applied to the cells to visualize the branching structure. For fixation, 4% paraformaldehyde (Electron Microscopy Science, Hatfield, Pennsylvania) was applied to the Matrigel at room temperature for 20 min. After washing with phosphate-buffered saline (PBS), PBS-containing 0.5% Triton X-100 was applied for cell permeabilization for 10 min at 4°C. This was followed by three 10-min washes with PBS. Then, the gels were blocked with 10% goat serum and 1% goat antimouse immunoglobulin G (Sigma-Aldrich, St. Louis, Missouri) in immunofluorescence buffer (0.2% Triton X-100; 0.1% bovine serum albumin and 0.05% Tween-20 in PBS). As a primary antibody, rabbit E-cadherin monoclonal antibody (Life Technologies, Grand Island, New York) was incubated with Matrigel overnight at 4°C and the gel was rinsed three times with 10% goat serum for 20 min each. Then, Alexa Fluor 555 goat antirabbit IgG was incubated for 2 h at room temperature followed by three 20-min rinses with PBS. For nuclear staining, DAPI was incubated for 20 min at room temperature followed by three 20-min rinses.
HUVECs were mixed with dextran-coated Cytodex 3 microcarriers at a concentration of 400 HUVEC per bead in 1 mL of endothelial cell growth medium (EGM)-2. The beads with cells were shaken gently every 20 min for 4 h at 37°C and 5% CO 2 , and then transferred to a 25-cm 2 tissue culture flask and left for 12 to 16 h in 5 mL of EGM-2 at 37°C and 5% CO 2 . In the following day, the beads with cells were washed three times with 1 mL of EGM-2 and resuspended at a concentration of 500 beads/mL in 2-mg/mL fibrinogen, 1-U/mL factor XIII, 0.04-U/mL aprotinin, and 80,000-cell/mL HDF at a pH of 7.4. Then, 250 μL of this fibrinogen/bead solution was added to 0.16 units of thrombin in one well of each of the glass-bottom 24well plates. The fibrinogen/HUVEC bead/HDF cell solution was allowed to clot for 5 min at room temperature and then at 37°C and 5% CO 2 for 20 min. EGM-2 was added to each well and equilibrated with the fibrin clot for 30 min at 37°C and 5% CO 2 . The medium was removed from the well and replaced with 1 mL of fresh EGM-2 and was later changed every other day. The coculture assays were monitored for 7 days and then fixed and stained with Alexa Fluor 488 for imaging of the cytoskeletons.

Zebrafish Embryo Culture
Transgenic zebrafishes (Islet1:GFP-mlcr:DsRed and cmlc2: GFP) were raised in the zebrafish core facility of the University of California, Los Angeles (UCLA). All experiments were performed in compliance with the approval of the GLA Institutional Animal Care and UCLA Institutional Animal Care and Use Committee protocols. To maintain transparency of the zebrafish embryos, they were incubated with egg water containing 0.2-mM 1-phenyl-2-thio-urea (Sigma) to suppress pigmentation at 24 hpf. The live fish embryos were anesthetized with lowconcentration tricaine (0.04 mg/mL, MS-222, Sigma) before being mounted in a fluorinated ethylene propylene tube for sustained imaging.

Optical Clearing of Thick Organs
The adult mouse brain (Thy1-GFP-M), neonate mouse heart (wild-type, αMHCCre; R26VT2/GK), and adult zebrafish heart (cmlc2-GFP) were originally turbid organs. As a result, it was necessary to perform tissue optical clearing before fluorescence imaging. An organic-solvent-based clearing method (uDISCO) 53 was used to clarify the adult mouse brain and zebrafish heart, and a hydrogel-based clearing (CLARITY) 2 method was used to clarify the neonate mouse hearts.