Figures and data
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).
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 renormalized 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 pattern. (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.
Purpose, requirements, applicability, and limitations of each step in the processing and analysis pipeline.
Description of the datasets used for training our custom StarDist3D model to detect nuclei.
All datasets were acquired with a two-photon microscope.
Benchmarking of dual-view two-photon imaging and spectral unmixing.
(a) Two-photon and confocal images of Hoechst-stained 77 h gastruloids at 100 µm depth. Insets show nuclei in the central 100×100 pixels. Scale bars are 50 µm. (b) Normalized intensity and FRC-QE in the center of gastruloids described in (a), as a function of depth. Shaded regions show standard deviation across six gastruloids of similar size. (c) Normalized intensity and FRC-QE obtained with two-photon imaging in the center (60 × 60 µm2 region) of 72 h Hoechst-stained gastruloids as a function of depth, mounted in either PBS, 20% optiprep, ProLong gold antifade, or 80% glycerol medium. (d) Normalized relative intensities detected in four detector channels as a function of imaging depth on aggregates stained with either Hoechst, AF488, AF568, or AF647. (e,f) Spectral unmixing of images acquired on gastruloids expressing transgenic H2B-GFP, additionally stained with Ecad-AF488. Raw data show significant cross-talk, particularly in the second channel. Unmixing using the spectra shown in (f) strongly reduces nuclei cross-talk in the Ecad-AF488 channel.
Quantification of StarDist3D performance.
(a) Training dataset of StarDist3D. Five samples with different characteristics were annotated in 3D, all the nuclei in a cropped region, with dimensions indicated in the table. (b) Comparison for each sample of the F1 scores obtained using either 3D or 2D ground-truths. 2D annotations were obtained on XY slices sampled every 20 µm. For each slice, we quantified the F1 score at IoU threshold = 0.5. To check for eventual biases due to the original anisotropy of the image, we also plot the F1 scores using 2D annotations in YZ and XZ (yellow and pink bars). Light blue bars represent the average of scores along all 2D planes. The F1 score using 3D annotations (dark blue bars) is well approximated by the F1 score using 2D annotations in XY, with a 5% difference at most. (c) As a consequence, we could quantify the segmentation score on new samples only by annotating in 2D the XY planes of two samples. 7 to 8 slices were annotated regularly throughout the volumes and the plots (d-g) show the 2D quantification of nuclei segmentation performance in depth. (d) F1 score as a function of depth in sample 1, comparing our contrast enhancement method in a local window (solid line) versus using a classical contrast enhancement method based on histogram statistics from the full image (dotted line), see Methods. The arrow shows the midplane, where the SNR is the lowest after fusion of both sides. (e) XZ view of the images and corresponding segmentation using either global or local contrast enhancement methods. The dotted line indicates the midplane. (f) F1 score as a function of depth in sample 2, comparing two-views imaging (solid line) versus one view imaging on the same sample (dotted line). We enhanced contrast locally in both samples to increase the signal in deeper z planes, but two-views imaging is still necessary to segment the nuclei deeper than 200 µm in the sample. (g) Cropped regions at 210 µm depth of a sample imaged from two sides, registered, reconstructed in 3D, and locally enhanced, versus the same image from only one side, far from the objective.
The wavelength-dependent normalization method.
(a) We evaluated how the intensity decreases in an ubiquitous signal as a function of depth and wavelength. Samples are immunostained with histones markers in the four colors. The intensity is measured in a central column of the sample, to simplify the geometry and remove border effects. In this region, we compute the median intensity in each plane. (b) Plot of the intensity in each channel, normalized to 1. Shaded regions indicate standard deviation. An exponential function is fitted to the normalized intensity to extract the characteristic length of decrease d. (c) This length d is plotted for each sample as a function of their wavelength. (d) We extract the ratio r = dblue/d, dblue being the median of the characteristic length d in the blue channel. Mean and standard deviation are shown. (e) Plot of the intensity in the far-red channel as a function of the intensity in the blue channel, both axis are in logarithmic scale. Images are binned so that each black dot represents a 10 · 10 µm2 square. N=5 samples are stacked onto the plot.The slope of 0.46 validates the value 
Illustration of the intensity model and of the normalization procedure based on a ubiquitous signal.
(a) In equation (9), we model any arbitrary intensity signal Ib by decomposing it into three terms: (i) a continuous field A representing local intensity artifacts, (ii) a field Sb representing a continuous depiction of the biological signal, and (iii) an instance mask Minst, which localizes the biological objects (e.g. nuclei or membranes) from where the biological signal originates. The first two fields are assumed to have variations at length scales equal or above the typical cell size. Finally, an additive term ε which varies at sub-cellular length-scales and encompasses all the remaining voxel-scale signal (e.g. measurement noise or intra-cellular gene expression heterogeneities) completes the model. (b) An ubiquitous signal (here Hoechst) can be used to approximate the field of intensity artifacts, which can be removed from the original intensity signal.
Illustration of the sigmoid function that is used for weighted fusion of dual-view images,
see equation (3). The coordinate along the overlapping region between the two views is mapped to the range [0,1] to define the parameters of the sigmoid function independently from the gastruloid size. l is the decay length of the sigmoid, and we define the effective fusion width δ as the fixed fraction of the overlapping region where both signals are substantially present to create the fused part.
(a) Illustration of the reproducibility of the spatial colocalization of a sharp low-to-high density boundary with a region of high nuclei alignment in XY and orthogonal views from several gastruloids. The leftmost column shows locally enhanced Hoechst signal in the XY plane, and a plane perpendicular to it. The red lines give the positions of the complementary plane. 120 h gastruloid, stained with ph3 and Hoechst, processed and segmented (nuclei and divisions) using the pipeline to generate the division density map and the cell density map. With these two fields, we compute a 3D proliferation map (shown here is a 2D section in the mid-plane).
Direct comparison of depth-dependent correlation of intensity between immunostained T-Bra signal and the intensity signal of a T-Bra reporter (H2B-GFP expressed under a T-Bra promoter) in gastruloids grown from a transgenic cell line fixed at 72 h.
(a) Image of a sample expressing the endogenous T-Bra reporter (green) and immunostained T-Bra (magenta). Signals are correlated almost everywhere (1) except for cells in mitosis (2) (b) Comparison of XZ views of Hoechst, endogenous T-Bra and immunostained T-Bra, before and after wavelength-dependent intensity normalization using the Hoechst signal as reference. (c) Intensity correlation heatmaps between the immunostained and endogenous T-Bra signals, in the full sample. After normalization, the two signals are directly equal up to a small additive bias. Vertical and horizontal black lines represent the median of each signal. (d) Intensity correlation heatmaps computed in three 3D slices of width 50 µm at varying depths. Normalization leads to a totally homogeneous correlation relationship between the two signals, with the same slope and additive bias at all depths.
Effect of spectral filtering and local intensity normalization on 3D gene expression analysis.
(a) As illustrated in Fig.5d for FoxA2, we analyze the radial distribution of Sox2 in 60 h gastruloids. Below, we present 3D visualizations of the normalized Sox2 signal alongside the segmentation of Sox2+ cells. The resulting 3D maps of the local fraction of Sox2+ cells reveal no preferential radial pattern. (b) Comparison of lateral intensity profiles at different depths. We compute the intensity across rectangular sections of each image in 2D, as shown on the left. The graphs show the line profiles at different depths for the Hoechst channel (top) and the T-Bra channel (bottom), comparing profiles before (left plot) and after (right plot) normalization. Color code stands for depth (darker colors correspond to deeper z-sections). We observe for non-normalized data that z-planes at mid-depth have a low intensity that is corrected in the normalized case. Hoechst profiles are flattened throughout the section after normalization, as expected. (c) Illustration of plots (d) and (e) : for each z-plane, we measure the median intensity in the Hoechst channel, within the sample, and plot it as a function of its depth. We use this method to quantitatively check that the intensity in the nuclei channel is constant in depth and does not present gradients after normalization. (d) Using the method described in (c), we show the median intensity as a function of depth for 4 samples from different acquisitions, corresponding to different timing and/or diameters. We compare the intensities before (dotted lines) and after normalization (solid lines), and observe that the normalization scheme robustly flattens the intensity profile. (e) For one sample, we compare different values of σ which is the standard deviation of the Gaussian kernel used to normalize the data. The images are respective maps of nuclei intensities for σ = 10, 20 and 40 µm. Using the method described in (d), we show the median nuclei intensity as a function of depth for these three σ values. f Correlation maps between FoxA2 and Oct4 for spectrally non-filtered data (top) versus filtered data (bottom), and the corresponding images shown on the left. Besides filtering, images were still registered and normalized the same way. (g) Correlation maps between FoxA2 and T-Bra for non-normalized data (top) versus normalized data (bottom), and the corresponding images are shown on the left. Images were still filtered and registered the same way.
Qualitative benchmark of our segmentation model.
All images were segmented in 3D, and we show here one XY plane for simplicity. (a) Qualitative assessment of the performance of different segmentation models against stardist-tapenade on gastruloid images (Hoechst staining). The pre-trained Cellpose-SAM [70] model gives high-quality results, which are further investigated in Fig. S10. (b) Qualitative comparison of our custom Stardist3D segmentation strategy (stardist-tapenade) on diverse published 3D nuclei datasets proves the versatility of our model. (1) Gastruloid stained with the nuclear marker DRAQ5 imaged with an open-top dual-view and dual-illumination LSM [62]. (2) Breast cancer spheroid [29]. (3) Primary pancreatic ductal adenocarcinoma organoids imaged with confocal microscopy [29]. (4) Human colon organoid imaged with LSM laser scanning confocal microscope [29]. (5) Monolayer HCT116 cells imaged with LSM laser scanning confocal microscope [29]. (6) Fixed zebrafish embryo stained for nuclei and imaged with a Zeiss LSM 880 confocal microscopy [63].
Validation of the Cellpose–SAM pretrained model for nuclear segmentation and morphology in gastruloids.
(a) The Tapenade preprocessing pipeline improves Cellpose–SAM [70] outputs, especially in deep planes (z > 150 µm), by using local contrast enhancement. (b) Compared to our StarDist3D model (stardist-tapenade), Cellpose–SAM better preserves elongated nuclei (green arrows), while StarDist3D tends to under-segment by rounding objects. Cellpose–SAM can, however, yield partial segments in low–signal-to-noise regions (red arrows). Overall this model seems highly competent for nuclei segmentation in gastruloids. Here and in the rest of this figure, we used the same sample as in Fig. 4 to allow for direct comparison. (c) Tissue-scale, packing-related metrics from Cellpose–SAM labels qualitatively match those from stardist-tapenade, though the more accurate boundaries increase estimated nuclear volume and volume fraction. The true-strain major-axis field remain similar because StarDist3D mostly preserves orientations. (d) Quantification of the region where the cell-density gradient and the true-strain major axes are anti-aligned is consistent when using Cellpose–SAM. Vector magnitudes of both fields are displayed to aid interpretation in areas with small gradients or deformations.
