Overview of the experimental and computational pipeline for in toto imaging of organoids at cellular resolution and multiscale analyses.

Dual-view multi-color two-photon imaging.

(a) Multiple gastruloids are mounted in glycerol clearing medium in between two glass slides using spacers. (b) Mounting gastruloids in 80% glycerol significantly improves imaging at depth compared to PBS mounting. Insets show images at different depth in a central square region of 60×60 µm2. (c) Quantification of intensity, FRC-QE image quality metric and density of segmented nuclei as a function of depth in the central region depicted in (b) for both glycerol- and PBS-mounted samples. Shaded areas show standard deviation for six (PBS) and nine (glycerol) gastruloids. (d-f) Top and bottom positions are assigned using a pattern matching algorithm that is robust to residual movements of individual gastruloids, as illustrated for position 6 (d). Respective image stacks are registered using rigid 3D transformations and fused applying sigmoid fusion weights (e). This leads to in toto coverage of whole gastruloids, in which individual cells are still visible in the mid-plane (e), as shown for a gastruloid of typical size of 250 µm diameter. The average intensity per plane and FRC-QE are shown in (f) for one- and dual-view images. (g) Polar histogram plots of rotation angles around X-,Y-, and Z-axis from the 3D registration of dual-view top and bottom acquisitions, pooled for 37 aggregates. The sample slide was flipped around the X-axis. (h-j) Multi-color detection using two-photon excitation can lead to significant spectral cross-talk, as shown in the apparent emission spectra detected on gastruloids stained with a single fluorophore species (h). Using spectral unmixing, false-positive signal due to spectral cross-talk is effectively removed, as shown in the raw and filtered images (here at 50 µm depth) of gastruloids stained with T-Bra-AF568 and Sox2-AF647 (i) and the quantification of the intensity along a line (10 pixels width) across few nuclei for raw and spectrally filtered images (j). (k) Four-color images of the mid-plane of 77 h gastruloids stained with Hoechst, FoxA2-AF488, Oct4-AF568 and T-Bra-AF647, obtained with dual-view four-channel two-photon imaging and spectral unmixing. The gastruloids were acquired in the same experiment. Scale bars are 50 µm in (b,i) and 100 µm in (k).

Description of the datasets used for training our custom StarDist3D model to detect nuclei.

All datasets were acquired with a two-photon microscope.

3D nuclei segmentation with StarDist3D.

(a) Exemplary segmentation (contours) of Hoechst stained nuclei (gray) using custom trained StarDist3D model. Insets show segmentation results in inner and outer cell layers, at 150 µm depth. (b) Three different datasets were annotated and used for training. Images show exemplary annotated patches (top row: intensity images, bottom row: label images). Before training, images were resampled to a common voxel size of 0.62 × 0.62 × 0.62 µm3. Image contrast when enhanced using local histogram clipping and normalization. During training, data augmentation was used to maximize training input. (c) To evaluate segmentation performance, entire z-planes of a gastruloid were annotated at different depths in 2D, on a 250 µm size gastruloid imaged from two views and fused, and the segmentation result was compared with the ground truth. (d) Top : Precision and Recall at IoU threshold 0.5 as a function of depth. Bottom: F1 score as a function of IoU thresholds, color represents depth. (e) Cell density obtained from 3D segmentation of gastruloids of different diameters. The corresponding images at the mid-plane are shown at the bottom. (f) 3D rendering of the detection of mitotic cells using StarDist3D. Top image shows ph3 stained sample, bottom image shows detected cell divisions overlayed with the ph3 signal.

Cell-to-tissue scale morphological analysis.

(a) Segmenting individual nuclei instances in 3D allows a complete characterization of object packing by quantifying cell density, nuclear volume fraction, and nuclear volume. Each quantity is shown at three resolutions ranging from the cell scale to the tissue scale, which correspond to convolution with Gaussian averaging kernels with σ = 6, 12, and 25 µm, respectively. (b) Divisions stained with ph3 and nuclei stained with Hoechst are segmented using StardDist3D and used to compute respectively maps of division density and cell density. The proliferation is defined as the ratio of the division density with the cell density, showing the local fraction of cells undergoing mitosis. (c) Cell-scale fields of deformation, quantified by the major axis of the nuclei true strain tensors. A version sampled on a regular grid is provided for better readability. The deformation field reveals a region of spatially persistent alignment that overlaps with the density boundary and highlights nuclei forming the lining of a lumen. (d) We compute the tissue-scale field of squared cosine of the angle between the cell density gradient and the true strain tensor major axis to quantify local regions of anti-alignment. The vector magnitudes of the two fields are shown to add perspective to the analysis in regions where the gradient or the local deformation is small. All quantities shown in this figure are computed in 3D, but we show a single z-slice (vector quantities are projected onto the slice) for convenience.

Coarse-grained and cell-level gene expression analysis.

(a,b) Images of Hoechst and T-Bra in the mid-plane of 77 h gastruloids. Using the nuclei segmentation, maps of averaged nuclei intensity were generated and used to locally re-normalize the T-Bra signal. The graphs on the right side of each panel show the intensity across rectangular sections of the mid-plane, shown in the images as blue/ brown boxes, before and after intensity normalization. This intensity profile is smoothed using a Gaussian kernel on the images to extract the large-scale intensity variations. (a) example of a radially symmetric gradient of nuclei intensity, from the exterior to the interior. (b) example of an asymmetric gradient, the left side of the sample has a high nuclei intensity that is corrected after normalization, changing the T-Bra profile (bottom right graph). (c) Images of 77 h gastruloids (mid-plane) of different size, immunostained for T-Bra and FoxA2 gene expression markers. From the re-normalized 3D images, radial intensity plots were generated for big (291 ± 15 µm average diameter, n=5) and small (193 ± 20 µm average diameter, n=8) gastruloids as a function of distance to the border. Shaded regions show standard deviations. (d) For one of the previously described samples fixed at 77 h, FoxA2 positive cells are visualized from the re-normalized data using the interactive napari tools and the image thresholded accordingly. Using the binary image of FoxA2 positive cells, a map of local fraction of FoxA2 positive cells, rescaled by the local cell density, was computed, showing a radial patern. (e) Analysis of the radial distribution, i.e. fraction of FoxA2 positive cells as a function of distance to the border,for big (278 ± 16 µm average diameter, n=3) and small (193 ± 20 µm average diameter, n=8) aggregates, computed from the coarse-grained maps as shown in (d). Shaded regions show standard deviations. (f) Three-color image of a 77 h gastruloid immunostained for T-Bra, Oct4, and FoxA2. The three plots show a cell-level correlation analysis of the three markers for the entire gastruloid. Each dot in the correlation plot represents an individual cell. Shown are all three marker combinations. Quadrants show regions of background signal (−) and actual expression (+) for each marker in the respective color, defined by Otsu thresholding.

Interactive napari plugins for bioimage handling, manual registration, preprocessing, data exploration and analysis.

(a) A napari plugin for exploring and processing large nD bioimages. Users can load datasets lazily, inspect slices along a chosen dimension, and export them as separate TIFF files. Compatible with common bioimage formats, it simplifies handling complex datasets. (b,c) To fuse data acquired with our dual-view setup, we provide a napari plugin to assist users in manually defining a 3D rigid transformation either when it is trivial (e.g a single axial translation) or to initialize an automatic registration tool close to the optimum. The first mode (left) allows the user to select salient landmarks and match them between the two views. An optimal rigid transformation is found automatically. The plugin’s second mode (right) lets the user define each element of the transformation (translations and rotations) until a satisfactory match is observed in 3D or on 2D planes. (d) A complete preprocessing suite is provided as a napari plugin to transform raw 3D and 3D+time datasets into datasets ready for analysis. The napari implementation allows for quick and visual exploration of parameters best suited to a given dataset. A recording feature can be toggled to save a complete user-defined pipeline into a JSON file, that can be used to process large datasets in batch or for sharing. (e) Interactive analysis of multiscale correlation heatmaps in napari is provided by allowing users to change the analysis length scale and to draw regions of interest directly on the correlation plot to see which cells contributed to the region’s statistics on a 3D view of the data.