Introduction

Understanding organoid development at multiple scales

The biology of organoids and tumoroids represents a frontier in fundamental biology and biomedical research, driven by the need for accurate 3D models complementary to 2D cell cultures and animal studies [1, 2]. The context of the “3Rs” (Replacement, Reduction, and Refinement) aims at minimizing the use of animals in research and enhancing ethical experimental practices [3]. Developing these models in 3D is crucial for mimicking the complex architecture and functionality of tissues, enabling deeper insights into developmental or disease mechanisms and treatment responses [4, 5]. Organoids models exhibit less stereotypic development compared to embryos, where cell identification is typically straightforward based on position within the embryo at specific developmental stages [6, 7]. In contrast, organoids present more complex developmental trajectories [8–11], characterized by variability in morphology, cell type composition, and differentiation levels. This variability may stem from differences in the initial cell differentiation state or number [12], as well as from a less controlled biochemical and mechanical environment than that of embryos [13], which grow confined within maternal tissues for mammals or in a rigid shell for insects [14]. This inherent variability underscores the need for detailed characterization of a sufficient number of organoids in parallel to properly characterize the diversity of their developmental processes. To this end, coarse-grained methods have been developed to classify phenotypes of intestinal organoids [15] and gastruloids [16, 17] by analyzing maximum projections from 3D immunofluorescence staining using high-throughput confocal imaging. These approaches are effective for systematically screening the impact of chemical compounds or signaling pathway perturbations on organoid morphology, as well as on the localization and proportion of various cell types and tissues within the organoid.

However, understanding how organoids or tumoroids develop and self-organize also requires capturing how individual stem cells or cancerous cells differentiate and modulate their behavior and morphology in response to their local microenvironment [18, 19]. Simultaneously, it is essential to quantify tissue characteristics on a coarser scale (where each measurement integrates data from groups of several neighboring cells) to capture the mechanical and biochemical interactions among different tissues and cell populations within the overall geometry of the organoid [20]. To achieve this, a challenge is to image these 3D tissues in toto at a resolution sufficient to identify each cell and its neighbors both in live and fixed samples [21]. On fixed samples, multiple specific fluorescent markers can be used to capture cell fate by measuring the expression levels of particular genes or proteins by immunostaining or in situ hybridization. Complementary, live imaging is essential for understanding how local cellular events (division, migration, rearrangements) underlie the morphogenesis of these systems, as it has been done previously for embryonic development [22, 23].

Imaging organoids in toto at cell resolution

Recently, light-sheet and confocal imaging systems have been optimized for parallel imaging or long-term live imaging of organoids [24–27]. Using deep-learning-based segmentation, tracking, and classification tools, these protocols generate digital atlases of developing organoids at single-cell resolution. For example, [25] employed a dual illumination inverted light-sheet microscope to live-image gut organoids over several days, enabling cell lineage reconstruction from nuclei tracking. In [26], a single-objective light-sheet microscope was used to image tens of 3D tissues, such as gut, hepatocyte, and neuroectoderm organoids, in parallel for phenotypic classification.

However, light-sheet and confocal imaging approaches remain limited to relatively small organoids (typically under 100 micrometers in diameter) [24, 28, 29] or hollow organoids with cavities, lumens, and epithelial structures [25], making them less optically opaque than larger organoids like gastruloids [30, 31], neuromuscular organoids [32], or cancer spheroids [33], which can reach diameters of 300 microns and more.

For large organoids, multiphoton microscopy provides a powerful alternative due to its ability to penetrate deep into thick tissues with minimal photodamage [34]. This technique utilizes longer wavelengths of light to excite fluorescent molecules within the specimen, allowing the visualization of cellular structures and interactions in high resolution [35]. Furthermore, it avoids drawbacks of confocal or light-sheet microscopy on large, densely packed samples, such as strong intensity gradients, image blurring, and reduced axial information due to light scattering, aberrations, degradation or divergence of the light-sheet [36].

Gastruloids: a toy model for an integrated pipeline to analyze organoids development in toto at cell resolution

We have recently developed two-photon live imaging protocols to image developing gastruloids over multiple hours [31, 37]. Gastruloids are mouse embryonic stem cells that self-organize into 3D embryonic organoids, sharing similarities with their in vivo counterparts [30]. Recent studies have shown that, within a few days, gastruloids undergo significant morphological changes, developing structures that closely resemble organs both genetically and morphologically, such as neural tube-, [38] gut-, and cardiac-like structures [39]. Live two-photon imaging of gastruloids around their mid-plane has enabled us to propose biophysical models for their early symmetry breaking and elongation, using coarse-grained 2D analysis of collective tissue flows and gene patterning [37]. However, investigating in a full 3D context the relation between cellular fate and local tissue architecture during gastruloids morphogenesis [40] is still a key challenge. Indeed, it necessitates imaging in toto these dense 3D cell aggregates, which have a typical diameter above 200 microns and are highly light-diffusive objects.

In this work, we present an integrated pipeline to: (1) image multiple immunostained and cleared gastruloids in whole-mount configuration using two-photon microscopy, capturing both in toto and cellular scale details; (2) process and filter the resulting 3D images to correct optical artifacts from deep imaging and segment individual cell nuclei; and analyze gene expression patterns, cellular events, and morphologies in 3D across multiple spatial scales (Fig.1).

Fig. 1.

I. Results

An experimental and computational pipeline for in toto organoid imaging

We have developed a versatile pipeline to image, process and analyze tens of immunostained gastruloids (with diameters ranging from 100 to 500 µm). The experimental pipeline consists of sequential opposite-view multi-channel imaging of cleared samples, imaged with a commercial two-photon microscope. Afterwards, multiple processing steps are performed to extract meaningful information from the acquired images: i) spectral unmixing to remove signal cross-talk, ii) dual-view registration and fusion to reconstruct in toto images, iii) sample and single cell segmentation, and iv) signal normalization across depth and channels. Finally, reconstructed in toto images are analyzed at different scales, from cell-level correlations of different genes to coarse-grained maps of gene expression patterns, cell shapes, densities, and division events. All processing and analysis tools are open-source and Python-based. They are openly accessible as notebooks and interactive napari plugins.

In toto multi-color two-photon imaging

To facilitate dual-view imaging, we mounted immunostained gastruloids between two glass coverslips using spacers of defined thickness (typically 250−500 µm, adapted to match the size of the gastruloids without compressing them), allowing us to image the sample iteratively from two opposing sides (Fig.2a, see Methods III A 3). To optimize imaging performance, we compared several refractive index matching mounting mediums, glycerol [41], ProLong Gold Antifade mounting medium (ThermoFisher), and optiprep [42], on gastruloids labeled with Hoechst nuclei stain (see Methods III A 1). We identified 80% glycerol as the mounting medium with best clearing performance (Fig.2b, Fig.S1c), leading to a 3-fold/ 8-fold reduction in intensity decay at 100 µm/ 200 µm depth compared to mounting in phosphate-buffered saline (PBS) and superior performance to gold antifade and live-cell compatible optiprep medium (Fig.S1c). At these depths, information content was significantly improved by 1.5- and 3-fold, quantified using Fourier ring correlation quality estimate (FRC-QE) [43] (Fig.2c). Segmenting cell nuclei in images acquired on glycerol-cleared samples reliably detected cells at depth up to 200 µm, whereas a continuous decline in cell density was observed for PBS mounting, with four times fewer cells detected at 200 µm depth (Fig.2c, see details on the segmentation in the next section). These degradation effects were much more pronounced in confocal imaging even in glycerol clearing, which showed 2× lower intensity and 8× lower FRC-QE than two-photon imaging at 100 µm depth (Fig.S1a,b).

Fig. 2.

Nevertheless, two-photon microscopy in glycerol beyond 200 µm depth suffered from decreased signal-to-noise-ratio (SNR) (Fig.2b,c). To properly image larger gastruloids in toto, a second detection side was required. We, therefore, flipped the sample slide and re-imaged the same set of gastruloids from the opposing side. Whereas the flipping was done manually, all mounted gastruloids were imaged automatically after defining their positions in the imaging software. A typical acquisition took about 5-10 min per gastruloid per side for 1 µm z-spacing and full field-of-view (318 µm with 0.62 µm pixel size at 40× magnification). To reconstruct in toto image stacks, we automatically identified opposite views of the same gastruloid by applying a pattern matching algorithm to gastruloids positions in the two acquisitions (Fig.2d, see Methods III F). Next, we registered each pair of image stacks acquired on the same aggregate using a rigid 3D transformation determined by a content-based block-matching registration algorithm [44]previously applied for multi-view and time registration of light-sheet imaging data [6, 7]. Finally, the registered image stacks were fused using a sigmoid decay function to weight contributions of the two sides (Fig.2e,f and S5). The registration was performed on the ubiquitous nuclei channel, i.e. Hoechst staining. In some instances, flipping of the slide resulted in additional minor rotations (on average (−13 ± 11)^°, (3 ± 18)^° and (5 ± 10)^° around the X-, Y-, and Z-axis for n=37 aggregates), which occasionally prohibited convergence of the registration algorithm (Fig.2g). For these cases, we implemented a napari plugin to pre-register the images manually and visualize the result in 3D (Fig.6). Fused gastruloid images exhibit about 4-fold lower intensity in the mid-plane compared to the outer ends, but approximately constant FRC-QE across depth (Fig.2f).

A common limitation in fluorescence imaging is the low number of targets that can be excited simultaneously and spectrally discriminated. To image multiple targets in parallel, we excited four fluorophore species (Hoechst, AF488-, AF568-, AF647-tagged secondary antibodies on immunostained gastruloids) simultaneously. Thanks to the broad two-photon excitation spectra, the four dye species could be excited with just two laser lines, 920 nm and 1040 nm. However, detecting their emission simultaneously on four non-descanned detectors resulted in significant signal cross-talk, due to spectral overlap (Fig.2h-j). We circumvented this problem by applying spectral decomposition [45]. We calibrated spectral patterns, i.e. apparent emission spectra for each separate fluorophore species, in gastruloids stained with a single fluorophore species, and used these as reference patterns to decompose multi-color data (Fig.2h). Because of chromatic effects due to dispersion, spectral patterns were not constant at different depths in the sample. For example, shorter wavelength emission decays stronger with imaging depth than longer wavelength fluorescence (Fig.S3b). To account for this, we performed the unmixing with depth-dependent spectral patterns (see Methods III E). Thereby, we could significantly remove cross-talk, for example, false positive cells in the far-red detection channel that resulted from cross-talk of bright cells visible in the red channel (Fig.2i,j). Notably, spectral decomposition with four-channel detection was even possible for spectrally strongly overlapping fluorophores, which we verified on gastruloids with transgenic H2B-GFP expression in the nuclei, stained with E-cadherin (Ecad)-AF488 at cell membranes (Fig.S1e,f). We thus applied spectral unmixing to all data presented throughout this work. This way, we could generate cross-talk free four-color images of several biologically relevant markers imaged simultaneously in multiple gastruloids in the same experiment (Fig.2k, Suppl.Movies 1-3).

Single cell segmentation

Quantifying cell-scale gene expression and nuclei morphological properties requires accurate instance segmentation of individual nuclei. In early-stage gastruloid datasets acquired with two-photon microscopy, accurate segmentation of nuclei is hindered by high nuclei packing (Fig.2b), shape heterogeneity, and spatially varying SNR due to optical aberrations and scattering, especially deeper in the sample. We thus trained a custom Stardist3D [46] model to segment nuclei, as it excels with separating star-convex objects in close contact.

A mid-plane section of the StarDist prediction is shown Fig.3a, exhibiting a qualitatively good result on outer layers as well as inner part of the sample. To increase robustness, the network was trained on three datasets acquired with two-photon microscopy in different experimental conditions, namely on (1) live-imaging of Histone 2B (H2B)-GFP gastruloids, (2) live-imaging of sparsely labeled mosaic gastruloids (composed of a mix of 75% of non-fluorescently labeled cells and 25% of cells expressing H2B-GFP under the Brachyury gene promoter [31]), and (3) 96 h fixed samples stained with Hoechst. This combined dataset presented variations in labeled nuclei density, nuclei texture, and in voxel dimensions (Fig.3b, see table I and Fig.S2a). In total, we annotated 4414 nuclei in 3D and used data augmentation (see Methods III G) to enrich the training datasets (Fig.3b). Both during training and inference, we applied local contrast enhancement and normalization to all datasets by clipping intensity histograms computed in local boxes at the 1st and 99th percentile and by mapping intensities between 0 and 1. This step was crucial to homogenize the intensity distribution across and inside datasets during training, and improved segmentation quality by reducing the number of false detections at depths at which the SNR is reduced.

TABLE I.

Fig. 3.

To assess the performance of the training, we manually annotated 2D planes at different depths throughout a sample (see Fig.S2c, sample 1) and compared the result of StarDist3D to the ground truth, both in terms of quality of nuclei detection and volume reconstruction (Fig.3c). To evaluate the performance, we computed the recall and the precision (resp. quantifying how many relevant items are retrieved and how many retrieved items are relevant), and the F1 score as the harmonic average between the two. We show in Fig.3d that our trained model recovers a constant precision and recall across the full depth, and the F1 score is superior to 0.8 up to intersection over union (IoU) threshold 0.5. Comparing the recall and precision to the ones obtained with global contrast enhancement, i.e., histogram clipping based on the percentiles of the whole image, our local contrast enhancement scheme mitigates both false positives and false negatives in deeper z planes (Fig.S2d-e). It is still necessary to image from two views to recover all the nuclei, as shown Fig S2f-g, even using the local contrast enhancement method. Our score measurement is based on computing the IoU of 3D objects based on 2D images. As manually annotating many dense organoid datasets for validation is a slow and cumbersome process, we designed a validation method based on comparing the 3D prediction of StarDist3D with 2D annotations in several XY planes along the z-axis. We validate this approach on Fig.S2b using four different samples in which we chose XY planes spaced evenly along the z-axis, and on which nuclei were annotated in 3D. The F1 score obtained by comparing the 3D prediction and the 3D ground-truth is well approximated by the F1 score obtained by comparing the 2D predictions and the 2D sliced annotated segments, with at most a 5% difference between the two scores. Finally, we observed that the model preserves a constant cell density across depth in spherical gastruloids of up to 475 µm in diameter (Fig.3e).

Overall, our custom model achieved an F1 score of 85 ± 3% at 50% IoU. Using a variety of annotated datasets, dual-view imaging and a local contrast enhancement algorithm, we trained a powerful StarDist3D model that shows robust segmentation performances regardless of depth, sample size, and experimental conditions.

Maps of proliferation and cell density

Differential division and growth rates in tissues are key drivers of morphogenetic changes [47]. Therefore, quantifying tissue-scale gradients in cell proliferation is crucial for understanding the morphogenesis of gastruloids. To detect division events, we stained gastruloids with phosphohistone H3 (ph3) and trained a separate custom Stardist3D model using 3D annotations of nuclei expressing ph3. This allowed us to detect all mitotic nuclei in whole-mount samples for any stage and size (Fig.3f and Suppl.Movie 4). To probe quantities related to the tissue structure at multiple scales, we smooth their signal with a Gaussian kernel of width σ, with σ defined as the spatial scale of interest. From the segmented nuclei instances, we compute 3D fields of cell density (number of cells per unit volume), nuclear volume fraction (ratio of nuclear volume to local neighborhood volume), and nuclear volume at multiple scales. We plot each field on Fig.4a at three spatial scales ranging from the nuclei scale, i.e σ set to the average nuclear radius, to the tissue scale, i.e σ set to six times the average nuclei radius. Taking this smoothing approach, we computed fields of division density (number of division events per unit volume), that represent the local density of cells in mitosis, from the division labels detected on ph3 as described above. An increased volumetric division density can result from either elevated cell proliferation or higher cellular volumetric density. To distinguish between these effects, we computed a proliferation probability field that represents the local fraction of cells undergoing mitosis within the 3D sample (by dividing the division density field by the cell density field). A cross-section of each of these 3D maps is represented in Fig.4b, showing a typical heterogeneous pattern of division density and cell density. These maps can be visualized in parallel for several samples, which allows system-atic comparison of patterns for high-throughput acquisitions. Analyzing proliferation fields in 3D gives us insight on growth patterns during development and is key to understand and quantify gastruloid morphogenesis. Importantly, this quantification can be applied across various developmental stages, providing quantitative insights into growth phases under different stages or media conditions. A proliferation pattern for a strongly elongated 120 h gastruloid grown in a modified culture protocol is shown in Fig.S6b.

Fig. 4.

Morphometric analysis

Tissue-scale patterns of cell/nuclei deformations in developing embryos or organoids are good indicators of the existence of mechanical stresses which can contribute to tissue morphogenesis [48] and affect in return cell differentiation [49]. To explore changes in tissue structure and organization during gastruloid development, we studied spatial heterogeneities in cell densities and cell deformations emerging during the first morphogenetic phase of gastruloids development, which consists in a global tissue elongation between 72 and 96 h after their initial aggregation. We have shown in previous work [31, 37] that such elongation is associated with the polarization of the gastruloid, exhibiting a posterior pole rich in E-cadherin and constituted of a mixed population of both endoderm and nascent mesoderm tissues, and an anterior pole poor in E-cadherin constituted of cardiac and vascular mesoderm tissues which are mesenchymal. Here, to prevent any relaxation of cellular deformation caused by tissue fixation, we performed 3D acquisition of live gastruloids with an H2B-GFP marker. We illustrate the full approach on a single sample in Fig.4c-d, and show that the results are reproducible across different samples (see Fig.S6a). The cell density field revealed a high level of heterogeneity at the nuclei scale, but analysis at different scales always showed the presence of two tissue-scale populations of low and high cell densities located respectively in the posterior and anterior part of the sample, separated by a sharp boundary characterized by the length scale of a few nuclei. These two populations and the boundary were also visible, although less clearly, from the nuclear volume and nuclear volume fraction fields, with the high-density population showing a higher nuclear volume fraction but lower nuclear volumes. The nuclear volume fraction indicates regions of high packing located in the high-density population, with a maximum of 61%, close to the theoretical value for random close packing of spheres (64%). As extra-cellular space volume fraction is low in such sample [37], this indicates high nuclei-to-cell volume ratios, which suggests that nuclei shapes are a good proxy for cell shapes, and thus nuclei could be relevant in cell-scale mechanotransduction processes.

To quantify cell deformation, we first computed the inertia tensor from each segmented nuclei instance and obtained the true strain tensor [50] for each cell (see Methods III K), which is a size independent measure of object deformation. This justifies averaging these tensors among nuclei with large size variations by effectively removing the size bias. We plot the major axis (associated with the largest eigenvalue) of the true strain tensor for each nucleus on Fig.4 at the nuclei scale. The major axis is oriented along the direction of the largest deformation of the nuclei. Close to the low-to-high density boundary, we observed that nuclei are elongated in the direction of the boundary, i.e perpendicular to the gradient of cell density. We further quantified this observation by computing the angle θ between the gradient of cell density and the major axis of the true strain tensors. In Fig.4d, we plot the continuous field representing the quantity cos² (θ) averaged at the tissue scale. For random orientations of the major axis with respect to the gradient of cell density, the average value of cos² (θ) is 1/3, and the value approaches 0 or 1 for respectively perpendicular and parallel orientations. The field reveals a wide region with values below 0.1 around the boundary, confirming that nuclei are elongated along the tissue boundary. Other regions show comparatively smaller values of the local gradient amplitude or no tissue-scale alignment. These descriptions appear reproducibly in several samples, see Fig.S6a, and are consistent with an interpretation of the boundary as a region where mechanical stresses are generated between gastruloids anterior and posterior poles. This analysis also shows how nuclei shapes at small scales are appropriate proxies for deformation inside and around tissue micro-structures like lumen and cavities, which allows for an accurate description of the local epithelial architecture. At the tissue scale, nuclei shapes can serve as a proxy for tissue mechanical stresses in the context of tissues modeled as visco-elastic liquids with residual elasticity in cell deformation [51].

Quantitative analysis of gene expression at global and cell-scale level

A key interest in in toto imaging is that it allows quantifying expression gradients of genetic markers and proteins of interest across entire organoids. In the case of gastruloids, gradients of differentiation emerge during Wnt activation, prior or concomitant with symmetry breaking [16, 31, 52, 53]. Due to the previously mentioned limitations of light-sheet microscopy for in toto coverage and quantitative intensity analyses for gastruloids beyond 100 - 200 µm diameter, investigations of such gradients have thus far been restricted to analyses of low resolution 2D images [16, 52]. Because of this limitation, it is still unknown when gradients of differentiation are established during gastruloid development. To explore this, we analyzed images of gastruloids fixed at 77 h, labeled with Hoechst and additionally immunostained with T-Bra, a mesendoderm marker expressed upon Wnt mediated cell differentiation. Previous investigations with 2D imaging indicated a gradient of higher T-Bra expression in outer cell layers and lower expression in the interior [31].

Although our dual-view imaging pipeline delivers in toto coverage, images of the ubiquitous Hoechst nuclei stain still exhibited intensity gradients along the z direction and between outer and inner cell layers (Fig.5a). We attributed these gradients to optically induced heterogeneities and developed a correction method that re-normalizes the channels of interest based on the local nuclei intensity. To this aim, we first corrected the intensity decrease, which is dependent on the emission wavelength and the depth (see Fig.S3). We applied on this corrected signal a convolution with a 3D Gaussian kernel masked with the nuclei segmentation. The result of that convolution is a coarse-grained map of nuclei intensity, that we use to locally normalize the original image (Fig.5a). Since this normalization was done in 3D, the intensity values in the mid-plane of a gastruloid were enhanced while the intensities in the first and last z-planes remained constant (Fig.S8b). Furthermore, the method effectively and robustly homogenized nuclei intensities across gastruloids from different acquisitions and imaging conditions (Fig.S8c and d). We evaluated different sizes for the Gaussian kernel of the convolution and chose an optimal one, with a typical value around the diameter of the nuclei (10 - 15 µm), that lead to the most homogeneous Hoechst signal after re-normalization (see Methods III I and Fig.S8e). This normalization scheme accounts for heterogeneities in gastruloids of arbitrary shape as well as local cell density and thus variations in the optical paths of excitation and fluorescence light rays that cannot be captured by a simple geometrical model.

Fig. 5.

When we applied the correction method to the T-Bra signal and evaluated re-normalized T-Bra intensity in a lateral section across the mid-plane of a 77 h gastruloid, we still observed a gradient with higher intensity in outer cell layers (Fig.5a), while the nuclei signal recovered a flat profile. Besides reflecting a true T-Bra expression gradient, a possible explanation for this observation could be hindered penetration of primary and secondary antibodies being considerably larger than Hoechst stain. To evaluate this, we compared immunostained T-Bra signal with the intensity 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. We found a strong correlation between the two channels at all depths, excluding penetration artifacts (Fig.S7).

In addition to optical artifacts, our re-normalization method also corrects for intensity variations effectively induced by residual rotations of aggregates after flipping, which tilt the optical axis with respect to the axis of the first acquisition (Fig.2g). This is illustrated in Fig.5b, providing an example where the nuclei signal exhibits an asymmetrical gradient due to the fusion of two opposite views that were subject to residual rotations. After re-normalization, the nuclei profile is homogeneous, whereas the profile in T-Bra drastically changed, showing a polarization in the opposite direction with respect to the initial gradient that was also visible both in the nuclei stain and T-Bra channel. This emphasizes the importance of our method to properly quantify gene expression gradients in 3D. The observed polarized T-Bra appearance is in agreement with our previous study, in which we have reported that T-Bra expression already commences to be polarized at 77 h in some gastruloids [31].

Because the intensity decay in depth is wavelength-dependent (Fig.S3b), normalizing the different channels only based on the Hoechst intensity in the blue channel would have induced over-corrections in depth. We applied a wavelength-dependent correction to the different channels, by first measuring the characteristic length of exponential decay for the four colors independantly (Fig.S3a and c), and second by computing a map of Hoechst intensity corrected by this characteristic length. This gave us the equivalent of coarse-grained maps of nuclei intensity in each color, serving as a reference to correct the different channels (details in Methods III I).

We then applied the correction scheme to four-color 3D acquisitions of multiple gastruloids of different size (obtained by varying the initial cell number of the culture) imaged in the same slide, immunostained with T-Bra, the pluripotency marker Oct4, and the early endoderm marker FoxA2 (Fig.5c). We sorted these gastruloids into two groups and observed an approximately constant T-Bra expression for small aggregates (193 ± 20 µm average diameter, n=8), whereas big aggregates (291 ± 15 µm average diameter, n=5) showed a clear decay of T-Bra expression towards the interior (Fig.5c). Notably, we quantified T-Bra intensity as a function of distance to the gastruloid border using the in toto 3D image stack, thus reflecting expression gradients in 3D. For FoxA2, a similar pattern was observed, but expression appeared in much fewer cells (Fig.5d). Because FoxA2 signal was sparse and noisy, we could not compute simply the map of expression as done previously for T-Bra. Instead, we computed a map of the fraction of FoxA2 cells in a given neighborhood. We first thresholded the signal in segmented nuclei, using the interactive visualization tools developed in napari (see Fig.6 and Methods). This allowed us to identify the correct intensity threshold and generate a binary mask of FoxA2 expression. We then computed the number of FoxA2 positive cells rescaled by the local nuclei density, using the same method as applied for proliferation analysis, except that we consider FoxA2 positive cells instead of ph3 positive cells. The coarse-grained map of FoxA2 positive cells was then analyzed radially for gastruloids of various sizes, small gastruloids of 193 ± 20 µm average diameter (n=8) and bigger gastruloids of 278 ± 16 µm average diameter (n=3). We observed slightly higher FoxA2 levels in outermost cell layers, and in general more cells expressing FoxA2 in small gastruloids proportionally to the total number of cells. We applied the same method of radial analysis on 3D maps of Sox2 fraction, see Fig.S8a on 5 gastruloids at 60h. Cells expressing Sox2 are heterogeneously distributed but there is no radial gradient. Taken together, these analyses demonstrate that coarse-grained maps coupled with radial profiling can quantify large-scale patterns and gradients of gene expression across gastruloids.

Fig. 6.

Our ability to image multiple expression markers in the same acquisition without cross-talk and to reliably segment all cells allows correlating expression of different genes at the single cell level. To demonstrate this, we analyzed correlations between T-Bra, Oct4, and FoxA2 in the data described above, after spectral filtering and intensity re-normalization. Both processing steps were crucial to obtain unbiased results (Fig.S8f,g). To discriminate marker positive from cells only showing background levels, we applied Otsu thresholding on the cells intensity distribution for each marker and confirmed the distinction by visible inspection.

We expected that FoxA2 and Oct4 are mutually exclusive (since cells are either in the process of endoderm differentiation or still pluripotent). A similar behavior was expected for T-Bra and Oct4/ T-Bra and FoxA2, although some co-expression could be present since T-Bra marks very early mesendoderm differentiation en route to meso- or endoderm and Oct4 is a relatively long-lived pluripotency marker [54].

We indeed observed that Oct4 and FoxA2 appeared as two mutually exclusive populations (Fig.5f, left). On the contrary, while the brightest Oct4 cells also showed low T-Bra signal and vice versa, a large pool of T-Bra and Oct4 positive cells appeared, not showing a striking correlation pattern for most cells (Fig.5f, center). For FoxA2 and T-Bra, a clear anti-correlation was observed, with cells appearing in all four quadrants (Fig.5f, right). Interestingly, we overall detected fewer FoxA2 positive cells than T-Bra positive cells. While we found many cells that were T-Bra positive and FoxA2 negative, most of the FoxA2 positive cells showed at least some T-Bra signal. This indicates that T-Bra expression might be required to activate FoxA2. However, cells expressing both markers still showed a negative correlation, suggesting that once FoxA2 is activated, T-Bra expression is repressed. Overall, these data suggest a sequence of cellular transitions from Oct4 progressively to T-Bra, from which a fraction of cells, presumably the ones expressing low T-Bra levels, will commit to FoxA2. To further explore such transitions and localize specific cell populations, tools that allow interactive mapping of cell populations selected in the correlation plot to the actual gastruloid are needed. We have, therefore, integrated this functionality into our napari plugin (see Fig6). We believe that our integrated imaging and analysis pipeline will allow researchers to study cell differentiation and cell state transitions at different stages of gastruloid development for the first time in full spatial context.

Community-driven tools

We packaged our preprocessing and analysis pipeline into a Python library called Tapenade (for “Thorough Analysis PipEliNe for Advanced DEep imaging”). The code is freely available on Github at github.com/GuignardLab/tapenade/. The library comes with an extensive documentation and Jupyter notebooks accessible to non-specialist users. All datasets used to train our custom StarDist models to detect nuclei and mitotic events, along with the optimized weights for nuclei and mitotic events detection are currently available in section IV. A non-specialist user can quickly segment a custom dataset by downloading the weights and loading them into one of the user-friendly implementations of Stardist like the Stardist Fiji plugin, the napari plugin stardist-napari, or the ZeroCostDL4Mic notebook [46, 55].

While working with large and dense 3D and 3D+time gastruloid datasets, we found that being able to visualize and interact with the data dynamically greatly helped processing it. During the preprocessing stage, dynamical exploration and interaction led to faster tuning of the parameters by allowing direct visual feedback, and gave key biophysical insight during the analysis stage. We thus created four user-friendly napari plugins designed around facilitating such interactions (Fig.6):

a. napari-file2folder

Efficient handling of large nD bioimaging datasets is a common challenge in bioimage analysis. To address this, we developed a napari plugin for lazy loading and efficient exploration of heavy bioimages using the zarr library. The plugin supports a wide range of standard bioimage file formats, including lsm, tif, ome.tiff, zarr, and czi. Upon lazy loading of an image, the user can immediately inspect its total shape across dimensions, even for datasets that cannot be fully loaded into memory. The plugin enables interactive selection of a single element along a specified dimension for focused exploration (e.g., viewing a single z-slice or a timepoint). Additionally, it offers a feature to save all elements along a selected dimension as separate tif files in a specified folder, which is particularly useful for downstream processing workflows with other plugins (Fig.6a). The plugin is freely available on the napari hub or at github.com/GuignardLab/napari-file2folder.

b. napari-manual-registration

When using our automatic registration tool to spatially register two views of the same organoid, we were sometimes faced with the issue that the tool would not converge to the true registration transformation. This happens when the initial position and orientation of the floating view are too far from their target values. We thus designed a napari plugin to quickly find a transformation that can be used to initialize our registration tool close to the optimal transformation. From two images loaded in napari representing two views of the same organoid, the plugin allows the user to either (i) annotate matching salient landmarks (e.g bright dead cells or lumen-like structures) in both the reference and floating views, from which an optimal rigid transformation can be found automatically using singular value decomposition [56] (Fig.6b), or (ii) manually define a rigid transformation by continually varying 3D rotations and translations while observing the results until a satisfying fit is found (Fig.6c). The plugin is freely available on the napari hub or at github.com/GuignardLab/napari-manual-registration.

c. napari-tapenade-processing

From a given set of raw images, segmented object instances, and object mask, the plugin allows the user to quickly run all preprocessing functions from our main pipeline with custom parameters while being able to see and interact with the result of each step (Fig.6d). For large datasets that are cumbersome to manipulate or cannot be loaded in napari, the plugin provides a macro recording feature: the users can experiment and design their own pipeline on a smaller subset of the dataset, then run it on the full dataset without having to load it in napari. The plugin is freely available on the napari hub or at https://github.com/GuignardLab/napari-tapenade-processing.

d. napari-spatial-correlation-plotter

This plugins allows the user to analyze the spatial correlations of two 3D fields loaded in napari (e.g. two fluorescent markers). The user can dynamically vary the analysis length scale, which corresponds to the standard deviation of the Gaussian kernel used for smoothing the 3D fields. If a layer of segmented nuclei instances is additionally specified, the histogram is constructed by binning values at the nuclei level (each point corresponds to an individual nucleus). Otherwise, individual voxel values are used. We took inspiration from the plugin napari-clusters-plotter [57] in letting the user dynamically interact with the correlation heatmap by manually selecting a region in the plot. The corresponding cells (or voxels) that contributed to the region’s statistics will be displayed in 3D on an independent napari layer for the user to interact with and gain biological insight (Fig.6e). The plugin is freely available on the napari hub or at github.com/GuignardLab/napari-spatial-correlation-plotter.

II. Discussion

Multiscale imaging and analysis pipeline

We have developed an integrated pipeline for in toto imaging of whole-mount gastruloids and automated image processing. A key innovation of our approach is its specific optimization for dense 3D tissues, such as gastruloids, which present similar imaging challenges as tumors and spheroids, unlike many organoids that are shallow and less light-diffusive. Unlocking deep whole-mount imaging is crucial for researchers working with such dense tissues, as cryosectioning strategies cannot fully substitute for comprehensive three-dimensional (3D) in toto analysis [18, 21]. Reconstructing 3D tissue microstructure from cryosectioned data is complex and challenging, involving slicing tissue into thin sections and imaging them individually. To achieve a 3D reconstruction, these 2D images must be meticulously aligned and stitched together, a process complicated by frequent issues such as tissue deformation and section loss. In contrast, whole-mount imaging preserves spatial relationships and tissue integrity, making it the method of choice for studying cellular interactions and tissue architecture in their native context. Additionally, whole-mount imaging reduces sampling bias and enhances the detection of rare cell populations or phenotypic variations that might be missed in thin tissue sections. This limitation also apply to state-of-the art spatial transcriptomic techniques, due to the fact that only few RNAs are typically captured per cell for each individual gene and many cells do not show any counts for a gene of interest expressed at low level [58].

A Python integrated pipeline for non-specialists, with napari-enabled 3D exploration

We built our pipeline in Python to leverage its rich ecosystem of scientific computing libraries while keeping the code as accessible and user-friendly as possible. We provide detailed installation instructions and Jupyter notebooks with step-by-step instructions designed to be used by non-specialists, with a focus on ease of use and reproducibility. To tackle the challenges that arise from manipulating deep and dense 3D data, we developped napari plugins for manual registration, preprocessing, and analysis. Napari provides a user-friendly graphical interface which enables real-time visualization and exploration of 3D data during each step of the preprocessing and analysis. This was essential for the development of the pipeline, as the direct visual feedback allowed us to quickly iterate and optimize the processing steps, while gaining insights into the 3D organization of gastruloids. We believe our pipeline offers a distinctive and effective integration of 3D data exploration and processing, providing a valuable tool for the organoid and tumoroids research community.

In this work we have not used refined masks of the gastruloids that would account for and exclude hollow structures like lumen. In Fig.4(a) and (c), we noted that this lead to noticeable patterns in the packing-related coarse-grained fields and in the local alignment of the true strain tensor major axis. Furthermore, obtaining precise masks is essential to extend our analysis to inherently hollow datasets, like intestinal gastruloids. To improve our current implementation, our masking approach could be augmented by a second method used to specifically detect these structures. Pixel classification based on Random-Forests could be used to get accurate segmentation using only a few manually annotated lines. These approaches are now widely accessible to non-specialist users in softwares like Fiji (with Labkit [59]), Ilastik [60], or napari (with napari-APOC [61]).

We did not explore continuously adjusting laser power based on excitation depth to enhance signals in deeper tissue layers, which could further extend the depth limits of our imaging method. In principle, incorporating such adaptation should not affect our pipeline, as long as the laser used for exciting the nuclei channel (which normalizes all other signals) is the same or maintains the same power-to-depth relationship as the lasers exciting the other fluorophores to be normalized. Here, we developed a wavelength-dependent normalization scheme based on a simplified columnar model. Future work could focus on creating a complete model that accounts for photon paths in a complex geometry, for example in elongated gastruloids, and integrate the spectral filtering approach to the wavelength-dependent normalization scheme.

In this work, we assume that the nuclear signal is the most reliable for normalizing other signals due to its extensive coverage of the tissue (over 60% in gastruloids). Instead of normalizing signals at the level of individual nuclei, we developed a correction scheme that identifies the optimal scale of characteristic gradients in nuclear intensities to generate a continuous field for re-normalization. This approach has the added advantage of being applicable to markers with different subcellular localizations (e.g., nuclei, cell membrane, intracellular organelles).

A necessary component of our quantitative imaging pipeline was the use of spectral unmixing to address spectral cross-talk, an inherent challenge in multi-color acquisitions, particularly with two-photon excitation. This allowed us to image four different fluorescent targets simultaneously and was required to capture and correlate multiple genetic markers on the single-cell level. This can be further expanded in the future by using state-of-the-art hyperspectral detection, which allows capturing up to 32 channels simultaneously and thus decompose even more spectrally over-lapping fluorescent targets [62]. Two-photon excitation is advantageous in this context, allowing for the excitation of multiple targets with a single laser line. In our work, three of the four fluorophore species (Hoechst, AF488, and AF568) were excited using a 920 nm laser line. Additionally, iterative rounds of immunostaining and washing can be applied to maximize the information content [24, 63]. In this context, the interactive registration tools that we have implemented here will be key ingredients to correct for sample movements between subsequent staining steps.

Exploring gastruloid development in 3D

In the context of gastruloid development, which serves as the central biological theme of this paper, it is crucial to understand how spatial heterogeneities in gene expression arise at both the cell microenvironment scale and the larger gastruloid scale. During early gastruloid symmetry breaking and elongation, gradients of T-Bra emerge. These radial gradients can be captured with our coarse-grained analysis (see Fig.5c). However, while both the endoderm marker FoxA2 and T-Bra exhibit radial gradients at the coarse-grained scale, their anti-correlated expression is only evident at the single-cell level (see Fig.5c-f). This example highlights the importance of a multiscale analysis pipeline that allows simultaneous visualization of tissue-scale properties and local cellular events. Our methodology utilizes segmentation of the Hoechst-stained nuclei channel to (1) identify individual cells for cell-scale analysis and (2) normalize other channels using the nuclei signal. By developing a diverse training set for nuclei segmentation, we achieved precise determination of nuclei shapes. This precision enables systematic studies of the relationship between nuclear morphology and gene expression patterns, and can be extended to cell morphology using membrane stains such as ZO-1 and cell shape segmentation methods like CellPose [64]. Nuclei deformation patterns show large-scale alignment at the gastruloid level (see Fig.4c-d) and localized alignment near lumens and cavities, providing potential insights into the structuring of epithelial and mesenchymal tissues during gastruloid development. In older gastruloids with cavity-containing structures, such as somitic structures [38], automated cavity detection will be crucial for accurately quantifying coarse-grained fields of cell or division densities.

While this work focuses on gastruloid datasets, we believe that the pipeline effectively extends to organoid, tumoroid, and other datasets that encounter challenges related to whole-mount deep imaging. We designed the library in the prospect of long-term maintenance, and it will receive continuous improvements in the future.

III. Materials and methods

A. Sample preparation

1. Gastruloid culture

Gastruloids were generated using the protocol described previously in [31], from four different cell lines. Briefly, cells were seeded and aggregated for 48 h in low-adherence 96-well plates (Costar ref: 7007) and subsequently pulsed with the Wnt agonist Chiron, which was washed out after 24 h, i.e. at 72 h of aggregate culture. Time points indicated throughout the paper refer to the time after seeding. The two cell lines T-Bra-GFP/NE-mKate2 and E-cad-GFP/Oct4-mCherry were used for gene expression analyses. The E14Tg2a.4 line (ATCC) was used for spectral calibration and benchmarking as it is devoid of any fluorescent construct. The H2B-GFP (a generous gift from Kat Hadjantonakis) was used for live imaging experiments and morphometric analyses. In brief, between 50 and 400 cells were seeded in low adhesive 96 wells plates in a neural differentiation medium. In order to promote endoderm formation and to reduce symmetry breaking variability compared to the standard gastruloids protocol [53], we enriched the medium with Fibroblast Growth Factor and Activin as in [31]. All cell lines were tested to be free from mycoplasma contamination using qPCR.

2. Immunostaining

Samples were immunostained using the same protocol as in [31] apart from the fact that samples were incubated in a 1 mM glycin solution, and that primary antibodies (see table below) were incubated during 3 consecutive days. Secondary antibodies AF488, A568, AF647 (ThermoFisher) and Hoechst 33342 (ThermoFisher) were chosen to be compatible with four-color two-photon imaging. They were diluted 500 times for incubation, except of Hoechst which was diluted 10000 times in order to minimize the crosstalk of the Hoechst signal into other channels.

3. Mounting

Samples were mounted between two glass coverslips of standard thickness (0.13 to 0.17 mm) using 250 µm or 500 µm spacers (SUNJin Lab, Taiwan, R.O.C.) in 20 µl or 40 µl mounting medium. For the clearing benchmark, gastruloids were incubated overnight and then mounted in either PBS, 80% glycerol solution in PBS (v/v) [41], 20% optiprep solution in PBS (v/v) [42], or antifade mounting medium (ProLong Gold Antifade Mountant, Thermo Fisher Scientific). For all other experiments with fixed samples, gastruloids were mounted in 80% glycerol solution in PBS (v/v). For live imaging, gastruloids were transferred from 96-well plates to either MatTek dishes (MatTek corporation, ref: P35G-1.5-14 C) or to micro-well plates (500 µm wells, SUNBIOSCIENCE ref: Gri3D) in culture medium and imaged in a chamber maintained at 37 ^°C, 5% CO₂ with a humidifier.

B. Two-photon imaging

Four-colors two-photon imaging of immunostained samples was performed on an upright Nikon A1 R MP microscope using an Apo LWD 40 × /1.1 WI λS DIC N2 objective. Fluorescence was simultaneously excited at 920 nm (to excite simultaneously Hoechst, AF488 and AF568) with a tunable femtoseconde laser (Chameleon Discovery NX from Coherent) and at 1040 nm (to excite simultaneously AF568 and AF647) with the same laser secondary output at fixed wavelength. Fluorescence was detected on four non-descanned GaAsP detectors with the following emission filters: 450/70 (Ch1), 530/55 (Ch2), 607/70 (Ch3), 700/50 (Ch4). To spectrally separate the channels, three dichroic mirrors (488 LP, 562 LP, 650 LP) were used. Multiposition imaging was used to automatically acquire image stacks on multiple gastruloids mounted in the same sample slide.

Live imaging of gastruloids was performed on a Zeiss 510 NLO (Inverse - LSM) with a femtosecond laser (Laser Mai Tai DeepSee HP) with a 40 ×/1.2 C Apochromat objective. The Histone-GFP signal was excited at 920 nm and detected with a non-descanned GaAsP detector with a 560 LP. The images were acquired with the full field-of-view (318 µm, pixel size 0.62 µm), and a z-spacing of 1 µm. Depth in z varies from 150 to 350 µm, depending on the acquisitions and the sample size, but always past the midplane of the gastruloid.

C. Confocal imaging

Confocal imaging was performed on a Zeiss LSM880 system (Carl Zeiss) using a Plan-Apochromat 40 ×/1.2 NA water immersion objective. Fluorescence was excited with the 488 nm line of an Argon laser. To split excitation and emission light, a 488 dichroic mirror was used. Fluorescence was detected between 491 and 700 nm on a 32-channel GaAsP array detector, with a pinhole set to one airy unit.

D. Image quality estimation

To quantitatively assess clearing efficiency and compare the information content of images at depth, we used the Fourier ring correlation quality estimate (FRC-QE) [43]. We selected a column of 400 ×400 pixels for the clearing benchmark and 100 ×100 pixels for all other analyses in the center of aggregates and computed FRC-QE using the available Fiji plug-in.

E. Preprocessing: Spectral unmixing

To remove spectral cross-talk, four-channel images were spectrally unmixed into n images corresponding to the contribution from each of the n fluorophore species present in the sample.

1. Spectral calibration

To this aim, reference spectra p_i of the utilized fluorophores were determined from calibration images acquired on aggregates stained with only one fluorophore species at the same excitation and detection settings. To account for wavelength dependent signal decay with imaging depths, spectral patterns were computed as function of imaging depth. It is important to note that spectral patterns depend on the gain settings of each detector. However, once calibrated for a certain set of gain values, reference patterns for different gains can be easily interpolated by simply measuring the count ratios for both sets of gains on a sample providing sufficient signal (i.e. well above background counts) on each detector for all gains.

2. Spectral decomposition

Four-channel images I^k(x, y, z) were decomposed using a previously published unmixing algorithm originally developed in the context of fluorescence lifetime correlation spectroscopy [65], image correlation [66] and fluorescence correlation spectroscopy [45, 67]. First, spectral filter functions were calculated using the calibrated reference spectra :

Here, P is a matrix with elements and D a diagonal matrix computed from the average image intensity in each channel k, D = diag[1/⟨I^k(x, y, z)⟩_x,y].

The spectral filter functions represent intensity weights per spectral channel. To obtain the filtered images I_i(x, y, z) for each fluorophore species i, the channel images are multiplied with the filter functions and summed over all channels (m = 4 here):

It should be noted that unmixing successfully works as long as signal levels are well above background counts and all images, acquired for spectral calibration and in actual imaging experiments, are devoid of detector saturation.

This way, even closely overlapping fluorophore species can be discriminated with only few detector channels, at the expense of increased image noise [66]. In this context, pure detector channels, i.e. channels that predominantly detect one of the fluorophore species, are beneficial. This is the case in the example of GFP and AF488 presented in Fig.S1e,f, in which the first channel detects mostly GFP fluorescence. In addition, samples in which signal from different fluorophore species is detected in different pixels (e.g. pixels corresponding to nuclei and membranes) are generally easier to unmix and suffer less from noise increase than samples were different fluorophore species contribute to the signal in the same pixels.

F. Preprocessing: Dual-view registration and fusion

Multiple organoids were acquired from the same slide, first from one side for all of them, then from the other side. Therefore, the first step of the registration and fusion algorithm was to pair the two sides of each organoid together. This pairing was done by extracting the sample positions stored in the metadata of the acquisitions and finding the pairing that minimizes the sum of the distances between paired samples using a linear assignment algorithm.

Registration of these dual-view image stacks was achieved using an existing intensity based blockmatching algorithm [44], adapted for 3D microscopy images and previously used in [6, 7]. From the two distinct images, one is arbitrarily chosen as the reference image and the other, named the floating image, is registered onto the reference one. Starting from an initial transformation, the blockmatching algorithm compares, on multiple hierarchical levels, each block in the floating image with the corresponding neighboring blocks in the reference image, and finds the rigid 3D transformation (specified by rotation around axes X, Y, Z and translation along axes X, Y, Z) from these sets of blocks. Then, we apply the transformation to the floating image, so that the two views are in the same referential, and fuse them into one image stack.

Because of the computational cost of comparing blocks in large 3D images, the size of the neighborhood that the registration algorithm explores is restricted. Therefore, the magnitude of the transformation between the two image stacks has to be relatively small (typically less than 30^° rotations and 100 units in translation). As a consequence, the user has to provide an initial transformation to the algorithm, which approximates the actual transformation. In the case of dual-view imaging with opposite views, the initial transformation would correspond to 180^° rotation around the X axis and some additional translation. To explore the initial transformation, we developed the napari plugin napari-manual-registration detailed above and illustrated in Fig.6. This tools allows to interactively rotate and translate the floating with respect to the reference image stack. The tool also allows users to specify manual landmarks in both image stacks, e.g. lumen, dead or dividing cells, from which an initial transformation can be automatically computed.

While the blockmatching algorithm generally allows for any type of transformation (rigid, affine, non-linear, …), we restrained the search space to rigid transformations. Moreover, since the stacks comprise multiple channels, we used a reference one that was expressed ubiquitously, here the nuclei staining, on which the transformation was computed. The other channels were registered using the same transformation matrix as the reference channel. Because the SNR of each stack decreased roughly with the distance to the objective, we weighted the contribution in the fusion of the two opposing sides of the sample. The weights f₁ and f₂ respectively for the reference view and floating view were computed as a sigmoid:

with z being the depth in the stack in voxels, z₀ the middle (i.e inflection point) of the sigmoid, where f₁(z₀) = 0.5, and p the slope of the sigmoid. The decay length l of the sigmoid and the effective fusion width δ are linked by l = 4/p and δ = 2l.

As illustrated in Fig.2e, the two sides were imaged up to a certain depth depending on the sample size and the image quality, past the mid-plane to ensure some overlap with the other view. Because we acquired the same depth and number of z-planes for both views, the sigmoid is centered in the middle of the overlapping region, which is why we considered z₀ = 0.5. The overlapping region to fuse had a typical width value of w = 70 to 150 µm (depending on the image depth) around the midplane. The coordinate along the overlap region was mapped to the range [0, 1] to allow setting the parameters of the sigmoid fusion function independently from the size of the overlap or the size of the gastruloid, as illustrated on Fig.S5. The slope of the sigmoid was set to p = 15, which leads to l = 0.27 and δ = 0.53. Intuitively, this means that 53% of the overlapping region was used to mix both signals with non-negligible mixing coefficients to create the fused part.

G. Nuclei segmentation

To segment nuclei, we trained a custom Stardist3D [46] model on three annotated datasets acquired with two-photon microscopy in different acquisition modalities: one fixed dateset, two live datasets (one with mosaic fluorescence labeling), see Table I and S2a.

1. Global contrast enhancement

Before training and inference with our Stardist3D model, we applied percentile-based histogram clipping and normalization to the input images. For a given image voxel at position , its intensity is first transformed as follows:

where and are low and high percentile intensity values in the image.

In practice, we choose and to be the first and 99th percentile of the intensity values in the image. As shown Fig.S2d, applying global contrast enhancement was not sufficient to restore sufficient intensity in the deeper regions of some samples, which lead to poor segmentation performances. We thus extended our contrast enhancement method to apply in local sub-regions of the images.

2. Local contrast enhancement

Because the image contrast is not homogeneous within and across samples, we homogenized intensity profiles across the training datasets by applying local percentile-based histogram clipping and normalization.

For a given image voxel at position , its intensity is first transformed as follows:

where and as before, the first and 99th percentile of the intensity values within the window, and w close to the average size of the nuclei in the dataset The transformed image Í^′ is finally clipped to the range [0, 1] to obtain the locally enhanced image.

The training datasets were rescaled to an isotropic voxel scale of 0.62 µm/pix in all directions, from which crops of 64 × 64 × 64 voxels were extracted and augmented during training by 3D rotations, flipping, intensity shifts, and additive Gaussian noise. Applying scaling transformations to crops during training is often regarded as the best way to achieve robustness in the case of a wide object size distribution. This is typically the case when several raw datasets acquired in different experimental conditions are combined, particularly due to different pixel sizes across datasets. However, we found in our case that the size distribution related to a single dataset is narrow enough that robustness is better achieved by first rescaling all training datasets to a common isotropic nuclei size. [64, 68].

Before inference, a new dataset simply needs to be rescaled to reach that common nuclei size.

Another StarDist model was trained to detect the divisions stained by ph3. To this aim, we iteratively trained StarDist3D models with a very small batch of annotated data, using the model trained on all nuclei, and corrected the prediction using napari to train again. At the end, 31 images of whole-mount gastruloids and 5300 divisions were used to train a custom Stardist3D model, using the same parameters as for the nuclei detection model in terms of preprocessing, crops size and data augmentation.

Whereas the model trained on nuclei was not very successful on ph3 stained nuclei, because of their particular shapes, the final model after multiple iterations achieved better results, around F1=82 %. Some over-segmentations around large and deformed nuclei were still present, and using 3D Gaussian blur did not achieved a better score, because they led to misdetections where neighboring nuclei were in mitosis. Therefore, after prediction, post-processing was used to filter out small volumes, events detected from the background or dead cells, and minimal manual curation was done, using napari.

H. Masked Gaussian convolution to probe spatial fields at different scales

To probe volumetric signals at different length scales σ ranging from the average nuclei size to the width of the organoid, we convolve signals with a Gaussian kernel of size σ. We propose two ways to compute convolutions depending on whether the signal is dense (e.g an image representing a signal intensity) or sparse (e.g nuclei volumes defined only locally on nuclei centroids).

a. Dense data

We apply convolution with a normalized Gaussian kernel to dense signals. For a given voxel at position , its value , whether scalar or vector, is transformed according to:

with , and a the convolution window size. In our implementation, we choose a = 3σ.

When a binary mask (e.g defining the inside of the organoid) is available, we apply a masked convolution with a correcting normalization factor. Let M be a binary masking function, e.g for voxels at positions inside the organoid, and 0 everywhere else, we compute

This additional normalization factor prevents the artifactual decay of the convolved signal close to the mask boundaries by “excluding” voxels outside the mask from the convolution.

Note that the masked formulation allows the user to compute values of a quantity outside its domain of expression. For instance, if both a mask M_inst of nuclear instances and a mask M_in/out that marks the inside and outside of the sample are provided, a quantity like the nuclear Hoechst expression can be computed in the full region defined by M_in/out, even though it only appears in the nuclei.

b. Sparse data

For sparse data, i.e quantities defined locally at n specific positions , we compute the convolution and masked convolution in a similar manner, but with a mask M that is non-zero only at the positions at which the sparse signal is defined (e.g nuclei centroids). This formulation allows the user to choose between studying the sparse signal at its original positions, i.e by computing the above convolution at positions , (ii) at resampled positions, e.g on a uniform grid for convenience, or (iii) as a dense continuous field across the whole domain of the gastruloid.

We provide a second implementation of the masked convolution specifically for sparse data, which is more convenient when the smoothed sparse signals should also be smooth (e.g when it should be defined at its initial positions or on a regular grid). This implementation uses the Kd-Tree data structure implemented in Scipy to efficiently find neighbor relationships from the positions at which the sparse signal is defined, which is used to optimize the computation of the convolution (8).

I. Preprocessing: Intensity normalization

Quantifying gene expression in 3D tissue samples from imaging data is hindered by various sources of noise and optical artifacts that induce large scale intensity gradients (e.g. scattering, aberrations). The main sources of artifacts are (local) variations in optical paths and cell density, that usually manifest as decreased intensities towards deep tissue regions. Artifactual gradients can also appear during the registration of two non-collinear views, which further complexifies the expected distribution of artifacts and resulting intensity gradients.

We model the intensity signal from a given channel as

where is a binary mask of the signal instances (e.g nuclei or membranes), i.e if is inside an instance and 0 otherwise, is the gene expression field, which provides a continuous representation of the cell-scale variations of the biological signal, is the field representing the artifactual intensity variations, and is an additive uncorrelated noise term which models both variations in the signal at scales below the cell-scale and acquisition noise. This decomposition is illustrated on Fig.S4.

We hypothesize that the Hoechst signal is globally homogeneous throughout the volume, i.e , and that the field of optical artifacts A has large scale variations (typically superior to a few cell diameters). In particular, direct observations of homogeneous signals like Hoechst with different emission wavelengths, as displayed of FigS3(a-b), indicate that optical artifacts are mostly dominated by exponential decay with depth with a wavelength-dependent decay length. A second order contribution comes from scattering effects related to heterogeneities in photon’s optical paths, cellular density, or to geometrical considerations, e.g proximity to the sample’s boundaries. For every wavelength λ and every signal, we thus model A to first order as

with d the wavelength-dependent decay length. The second term a on the right hand side corresponds to second order corrections to the exponential decay model. Note that, in accordance with Beer-Lambert law, the value d(λ) also likely depends on experimental parameters ℰ, like the nature of the mounting medium or the settings of the experiments, in a multiplicative fashion, i.e

so that in similar imaging scenarios, ratios of decay lengths at two wavelengths

only depend on λ₁ and λ₂.

To estimate the decay lengths across the range of wavelengths in use, we measured the intensity decay in samples with ubiquitous nuclei stainings in four emission wavelengths: Hoechst for the blue channel (λ = 405nm), H2B-GFP amplified for the green channel (λ = 488nm), SPY-DNA555 or the endogenous H2B-Tomato for the red channel (λ = 555nm) and DRAQ5 for the far-red channel (λ = 647nm). Stainings for different wavelengths were done on separated samples to avoid signal crosstalk between the channels. All samples were imaged with the same excitation and detection settings. For each sample, we extracted the intensity profile in a central column along the z-axis to minimize boundary artifacts, and fitted equation (10) to recover the decay lengths d(λ) (Fig.S3(a-c)). We observe a monotonous relationship between median values of d across different samples and λ. To confirm our measurements, we also computed decay lengths via a second method: notice that the model of equation (10) leads to the following identity for two ubiquitous signals I₁ and I₂ at wavelengths λ₁ and λ₂ (neglecting the noise term ε in equation (9) and the second order term from equation (10))

We imaged new samples with both Hoechst (λ = 405nm) and the Far-red nuclei stain (λ = 647nm). There was no significant spectral overlap between the two signals. By fitting equation (13), we obtained another estimate of the ratio r(λ = 405nm, λ = 647nm) ≈0.46, consistent with the value of 0.48 obtained with the previous method (Fig.S3e). In Fig.S3(d), we plot the decay length ratios with the decay length of the Hoechst channel d_blue. In our model, these ratios are independent of imaging conditions, and the monotonicity and smoothness of the trend suggests that the ratios d_blue/d(λ) in the range λ∈ [405nm, 647nm] could be estimated by interpolating our results.

In our experiments, it was not possible to image ubiquitous nuclei signals at the same wavelengths as relevant biological signals, while Hoechst staining was often easily available. We thus developed a correction scheme based on Hoechst staining to separate the contribution of artifactual gradients from true gradients in gene expression. The rationale of the scheme stems from the fact that equation (13) also allows us to relate two ubiquitous signals up to a multiplicative constant:

The Hoechst signal could therefore be used to predict the intensity profile from another “virtual” ubiquitous signal imaged at some wavelength λ ≠ λ_blue that coincides with the wavelength of the biological signal of interest.

First, we compute the Gaussian average of the reference Hoechst signal I_h at a scale σ using equation (8), in which we keep only voxel intensities that belong to the nuclei using a mask M_nuc:

Here, σ is chosen large enough to average out the additive noise ε but small enough that the structure of the field of optical artifacts A is preserved. The practical choice of σ is discussed below.

We then normalize other biologically relevant intensity signals I_b by the Gaussian-averaged nuclei field to get

with ε^′ the normalized version of the additive noise ε. This step effectively removes the field of optical artifacts from the signal. The choice of σ used for the Gaussian averaging of the ubiquitous signal depends on the length scale at which the field of artifacts A varies. Since this is not known in practice, we choose in our implementation the value σ that leads to the most homogeneous normalized ubiquitous signal across all z-planes. Formally,

with MAD the median absolute deviation and Z_max the depth of the last z-slice. Fig S3f shows the decay profile of a sample normalized using different values of σ, as well as the corresponding Gaussian-averaged nuclei fields.

With the aim of comparing normalized signals across different samples acquired with the same imaging setup, we further rescale the normalized signal by a global multiplicative factor designed to mitigate the sample-dependent impact of the normalizing factor in equation (16). We multiply each normalized signal by a multiplicative factor such that the median intensity of the raw Hoechst signal computed in a region of brightest intensity is transmitted to the normalized Hoechst signal. The rationale is that the median intensity in the brightest regions reflects the expected intensity in the absence of artifacts, e.g. in outer cell layers close to the cover glass.

We tested the normalization method on a sample stained with both Hoechst and a Far-red nuclei stain, using Hoechst to normalize the Far-red channel and then look at the intensity profile. Fig.S3e shows the depth decay for the unnormalized signal, as well as the signal normalized using the ratio r(λ_blue, λ_far-red) = 0.48 as found experimentally. On the same plot we show the profile for the ratio r = 1 that would model a situation where we do not correct for wavelength dependency. In the latter, the intensity is over-corrected in depth, whereas with the correct ratio we recover a flat profile in depth. Because there is an error associated with the measurement of r, we tested small variations around the value 0.48. The depth decay changes slightly, as shown on the right panel with a 10% difference above and below. XZ profiles for these two panels are shown on the right. In addition, we correct the Hoechst for its own intensity decay, in that case using r = 1 which recovers a flat profile in depth, this is shown Fig.S3g with the corresponding XZ views. The last sanity check is performed by plotting the profile not only along the Z axis but also along X and Y, to ensure the absence of any other gradients of intensity. This is Fig.S3h, plotted on the Far-red channel normalized with r = 0.48.

J. Density maps

To quantify the local abundance of a target object, e.g. nuclei or specifically mitotic nuclei, we study three spatial fields: the density (number of objects per volume), the volume fraction (percentage of local space filled by the objects), and the map of objects volumes.

All fields are first computed at the voxel scale. The voxel-scale density field and map of objects volumes are computed by replacing each object by a single voxel of value, respectively 1 or the object’s volume placed at its centroid. The voxel-scale volume fraction corresponds to a binary mask of the target objects, and can be obtained directly with semantic segmentation methods or by binarising segmented object instances as output by Stardist3D. The voxel-scale fields are finally convolved with a Gaussian kernel of the desired analysis scale.

These quantities give a complete description of the local object abundance, as the volume fraction and map of objects volumes inform on the impact of volume on the object density. For instance, a large spatial region can misleadingly display a constant volume fraction but have large variations in object density and local objects volumes if the latter two compensate.

K. Morphometric analysis

To quantify object-scale deformation, we fit an ellipse to each segmented object. This can be done efficiently by computing the inertia tensor I. For a given object composed of N voxels at positions, we first compute the positions of the voxels relative to the object’s centroid, with. The inertia tensor is then defined as:

with and δ_i,j the Kronecker delta. The eigenvectors and eigenvalues of the inertia tensor correspond to the principal axes and moments of inertia of the object. Notably, the principal axes of I are the same as the principal axes of the ellipse that best fits the object. The lengths of the ellipse’s semi-axes can be obtained from the moments of inertia: if (I₁, I₂, I₃) are the moments of inertia associated with the principal axes , the lengths of the corresponding semi-axii (L₁, L₂, L₃) are given by:

Computing the inertia tensor for each object leads to a sparse tensor field that can be studied at arbitrary scales using the sparse Gaussian convolution equation defined previously. For instance, it can be used to study local object alignments by focusing only on the largest semi-axis (the main axis of deformation of the object). In practice, quantifying local object alignment by averaging inertia tensors introduces a bias by over-representing large objects. To avoid this bias, we instead compute the true strain tensor from the inertia tensor. The true strain tensor T is defined as:

with V the matrix of eigenvectors of I and . The eigenvalues (T₁, T₂, T₃) of T represent the relative amount of deviation of each principal length from the radius of a sphere with the same volume as the object. Due to its relative nature, the true strain tensor is suited to compare the deformations of two objects with different sizes, and is thus appropriate for averaging in a local neighborhood.

Acknowledgements

This work is supported by the French National Research Agency (“France 2030”, ANR-16-CONV-0001 from Excellence Initiative of Aix-Marseille University - A*MIDEX and generic grant to P.-F.L. ANR-19-CE13-0022), the Fondation de la Recherche Médicale (to P.-F.L. EQU202003010407), a generic grant to S.T. ANR-22-CE30-0021 and an ERC grant (to P.-F.L., ERC SyG 101072123). We also acknowledge the France-Bioimaging Infrastructure (ANR-10-INBS-04). V.DE. acknowledges support by an HFSP long-term postdoctoral fellowship (HFSP LT0058/2022-L).

Additional information

IV. Code and data availablity

All datasets used to reproduce the figures of the articles, the trained Stardist3D weights, and the associated training datasets can be found in the following link: https://zenodo.org/records/14748083

The Tapenade package is available on GitHub: https://github.com/GuignardLab/tapenade.

The napari plugin napari-manual registration is available here: https://github.com/GuignardLab/napari-manual-registration

The napari plugin napari-tapenade-processing is available here: https://github.com/GuignardLab/napari-tapenade-processing

The napari plugin napari-spatial-correlation-plotter is available here: https://github.com/GuignardLab/napari-spatial-correlation-plotter

Contributions

S.T, P-F.L and L.G designed the project. S.T and V.D-E designed the experiments. A.G, V.D-E, S.T and A.R conducted the experiments and collected the data. V.D-E developed the spectral filtering approach, A.G and J.V annotated nuclei data for the training of the StarDist3D network. A.G and L.G adapted the registration algorithm, with input from J.V. and V.D-E.

J.V developed and packaged the pipeline and the napari plugins, with technical input from A.G, and conceptual inputs from A.G, S.T, L.G, V.D-E. A.G, J.V, and V.D-E analyzed the data, with inputs from S.T, L.G, P.R and P-F.L. S.T, V.D-E, A.G and J.V wrote the manuscript with inputs from L.G, P-F.L and P.R. All authors discussed the scientific results and their interpretation. A.G, V.D-E, J.V and L.G created the figures and visualizations, P-F.L, S.T, L.G and P.R acquired funding.

Abbreviations

• SNR: signal-to-noise-ratio

• IoU: intersection over union

• ph3: phospho-histone H3

• PBS: phosphate-buffered saline

• FRC-QE: Fourier ring correlation quality estimate

• Ecad: E-cadherin

• H2B: Histone 2B

Additional files

Supp movie 1

Supp movie 2

Supp movie 3

Supp movie 4

Supplementary Information