key: cord-0325397-njywc40i
authors: Trébeau, Céline; de Monvel, Jacques Boutet; Altay, Gizem; Tinevez, Jean-Yves; Etournay, Raphaël
title: Extracting multiple surfaces from 3D microscopy images in complex biological tissues with the Zellige software tool
date: 2022-04-22
journal: bioRxiv
DOI: 10.1101/2022.04.05.485876
sha: 2b298d80445f69ad240183783ffb47def5bda27a
doc_id: 325397
cord_uid: njywc40i

Efficient tools allowing the extraction of 2D surfaces from 3D-microscopy data are essential for studies aiming to decipher the complex cellular choreography through which epithelium morphogenesis takes place during development. Most existing methods allow for the extraction of a single and smooth manifold of sufficiently high signal intensity and contrast, and usually fail when the surface of interest has a rough topography or when its localization is hampered by other surrounding structures of higher contrast. Multiple surface segmentation entails laborious manual annotations of the various surfaces separately. As automating this task is critical in studies involving tissue-tissue or tissue-matrix interaction, we developed the Zellige software, which allows the extraction of a non-prescribed number of surfaces of varying inclination, contrast, and texture from a 3D image. The tool requires the adjustment of a small set of control parameters, for which we provide an intuitive interface implemented as a Fiji plugin. As a proof of principle of the versatility of Zellige, we demonstrate its performance and robustness on synthetic images and on four different types of biological samples, covering a wide range of biological contexts.

Efficient tools allowing the extraction of 2D surfaces from 3D-microscopy data are essential for studies 11 aiming to decipher the complex cellular choreography through which epithelium morphogenesis takes 12 place during development. Most existing methods allow for the extraction of a single and smooth manifold 13 of sufficiently high signal intensity and contrast, and usually fail when the surface of interest has a rough 14 topography or when its localization is hampered by other surrounding structures of higher contrast. 15

Multiple surface segmentation entails laborious manual annotations of the various surfaces separately. As 16 automating this task is critical in studies involving tissue-tissue or tissue-matrix interaction, we developed 17 the Zellige software, which allows the extraction of a non-prescribed number of surfaces of varying 18 inclination, contrast, and texture from a 3D image. The tool requires the adjustment of a small set of control 19 parameters, for which we provide an intuitive interface implemented as a Fiji plugin. As a proof of principle 20 of the versatility of Zellige, we demonstrate its performance and robustness on synthetic images and on 21 four different types of biological samples, covering a wide range of biological contexts. 22 23 Introduction. 24 The interplay between gene regulatory networks and physical forces in driving collective cell behaviors is 25 key to tissue morphogenesis during development and to tissue homeostasis throughout life. Recent 26 quantitative studies of epithelial morphogenesis have begun to unravel the basic cellular and physical 27

principles of tissue development, by providing the tools to integrate multiple scales of tissue dynamics [1-28 4] . These tools are instrumental to quantify how cell shape changes, cell divisions, cell rearrangements and 29

cell extrusions contribute to tissue remodeling, and to establish data-driven computational models of tissue 30 morphogenesis. 31

Quantitative analysis of an epithelium starts with the extraction of its apical surface from 3D-32 microscopy images (z-stacks of xy-optical sections) encompassing the volume immediately surrounding the 33 epithelium. However, this is a difficult task because this surface is usually not flat (it is best modelled as a 34 curved surface, or 2D submanifold embedded in 3D space), and is often surrounded by other biological 35 structures such as cell layers, acellular membranes, extracellular matrix, and vesicles that hamper its 36 visualization and reconstruction. To make the surface extraction tractable, these studies rely on specific 37

preparations of the specimen, allowing to expose the entire epithelial surface labeled with junctional 38 fluorescent markers to reveal the network formed by epithelial cell-cell contacts. Once the epithelial surface 39 has been extracted, automated cell segmentation and cell contour tracking tools can be used to follow the 40 dynamics of every cell within the epithelium. 41

Another challenging experimental limitation of these studies is that some of the structures 42 surrounding the epithelium exert external physical constraints that are known to critically affect epithelial 43 morphogenesis by directing cellular dynamics and signaling pathways [1, [5] [6] [7] . To understand the physical 44 forces controlling tissue morphogenesis, it is thus essential to also characterize how the dynamics of these 45 extra-epithelial surfaces relate to that of the epithelium (see [8] for review). This calls for the development 46 of dedicated tools allowing the automated extraction of information from several surfaces of interest in a 47

given sample, since the sheer volume of the data precludes any attempt at a manual analysis. 48

Several surface extraction tools have been developed, some of which are available as open access 49 software, such as PreMosa [9], FastSME [10], and LocalZProjector [11] . These tools focus on the extraction 50 of a single, near-horizontal epithelial layer, which is assumed to (i) be sufficiently smooth, (ii) show enough 51 contrast against surrounding background signals, and (iii) cover the entire image field-of-view. Specifically, 52 it is assumed that the fluorescent marker used to label the epithelial cell network should provide the highest 53 contrast in the image and allow to select it out from autofluorescent extracellular structures such as the 54 cuticle in flies, or other acellular membranes in mammalian epithelia. The surface is then localized using 55 heuristic algorithms based on the detection of the pixels of maximum contrast and/or brightness. However, 56

applying these tools on more complex biological images with several epithelia of weaker contrast often 57 leads to incorrect localization of the surface of interest, and its blending with the nearby unwanted 58 biological structures. 59

MinCostZ on the other hand, is the only available open-source tool that allows the extraction of up 60 to two surfaces from a 3D stack, and imposes explicit continuity constraints on the reconstructed surfaces. 61

MinCostZ surface extraction relies on a previously developed formulation of the problem as a graph-cut 62 optimization [12] . It is implemented as an ImageJ plugin [13] , taking as control parameters, the number of 63 surfaces to extract, the maximum slope and the range of distances allowed between the surfaces, as well 64 as some user-defined cost function that should reflect the characteristics of the surfaces in term of signal 65 intensity, contrast and texture. Despite its interest, this approach remains computationally costly and 66 difficult to apply in practice due to the non-trivial choice of the cost function and the need to know 67 beforehand the relative positions of the surfaces to be extracted. 68 with a single set of reconstruction parameters. A sensitivity analysis also reveals a high robustness of Zellige 97 against small variations of these parameters. This will make it a tool of choice in terms of versatility and 98 ease of use for the investigation of biological surfaces. 99

Results & discussion.

Proof of concept of multiple surface extraction on a synthetic image.

The implementation of Zellige is summarized in Figure 1 and in the Methods section ( Figure S1 ). Figure 2  102 shows the results produced by Zellige on a phantom 3D image [29] containing three distinct synthetic 103 surfaces generated as described in Supplementary note 2. The three surfaces are extracted with little errors. 104 We assessed the quality of the reconstruction by comparing each of the height-maps produced by Zellige 105

to the corresponding ground truth (GT) height-map, which is exactly known in this case (Figure 2A -C, and 106

Supplemental note 3). For the three surfaces, the reconstruction has subpixel accuracy over >99% of the 107 GT pixels (Figure 2D- construction rounds, respectively) were set to their default reference values (see Figure S2 and 116

Supplemental note 4) and did not need to be adjusted. 117

Thus, using a single set of control parameters, Zellige can extract multiple surfaces of various 118

shapes and textures with virtually no error, without requiring the user to provide information about their 119 number or relative position, nor about their shape or texture characteristics. 120

Performance of Zellige on biological samples. Over the past few decades, the Drosophila model has been invaluable to decipher the molecular and cellular 123 mechanisms underlying organ embryogenesis [30, 31] . Epithelium morphogenesis studies not only revealed 124 the importance of mechanical stresses (including stress boundary conditions) and planar polarity signaling 125 on cell dynamics to generate tissues of reproducible sizes and shapes, it also highlighted the importance of 126 extracellular matrix attachments in constraining the tissue stresses that guide patterning [1, 32] . At the 127 pupal stage, the fly undergoes dramatic remodeling of its larval organs into adult organs. Large scale tissue 128 flows initiate at a timing that coincides with molting, when the epithelium contracts away from the 129 overlaying cuticular sac, a protective acellular membrane that imposes mechanical boundary conditions to 130 the tissue. 131 surfaces of interest can be identified, with varying signal intensities, noise levels and features (Figures 3A-135  B) . The abdomen is formed of an epithelium (surface S2) overlaid by a cuticle (surface S1). Lying just beneath 136 these two surfaces, one can observe globular structures showing in some places a higher intensity than the 137 signal coming from the surfaces. The wing also consists of an epithelium of low intensity signal (surface S4), 138

and an overlying cuticle (surface S3). These two surfaces are relatively flat, except for surface S3 which is 139 very steep near one of its edges. 140 Figure 3C shows a 3D graphical representation of the height-maps reconstructed by Zellige (green) 141 and those reconstructed by an expert biologist (blue), taken as ground truth (GT). While these height-maps 142 clearly show greater roughness than those of the synthetic surfaces presented earlier, they could again be 143 obtained with a single set of control parameters that were adjusted manually with the Zellige interface (see 144

Supplementary Table S1 ). We observe an excellent match between the four reconstructed and 145

corresponding GT height-maps, despite the rather complex topography of surfaces S1 and S2 (with slopes 146 reaching up to 45°), and the near-vertical inclination of surface S3 at its boundary. Yet, small deviations may 147 be seen in the regions of highest slope of the surfaces. Some of these deviations are likely attributable to 148 uncertainties in the definition of the GT height-maps, whose accuracy depends on the expert. 149 Figure 3D shows the differences between the reconstructed and GT height-maps, plotted as color-150 coded error maps. These differences are <2 (in pixel units) for most pixels, while some regions of higher 151 error values can be seen locally in surfaces S1 and S2, and at the boundary of surface S3. Note that for 152 surfaces S2 and S4, which contain junctional epithelial meshes composed of larger and smaller cells, 153 respectively, the GT height-map encompass not only the mesh but also the interior of the cells, where no 154 junctional signal is detected. The distance calculated inside the cells is thus more subjected to intensity 155 fluctuations, especially for surface S2. Nonetheless, the RMSEs of surfaces S1, S2 and S4 are less than 1, 156 showing that on average the reconstructed height-maps match the corresponding ground truth with 157 subpixel accuracy. The higher RMSE (1.25) of surface S3 is largely due to the region of steep region at the 158 edge of this surface (yellow region on the error map for this surface, Figure 3D ). The coverage of the 159 reconstructed height-maps is excellent (≥ 96%) for surfaces S1 and S2, and slightly lower, but still very good 160 (≥ 93%) for the smaller surfaces S3 and S4. Figure 3E shows the 2D projections of the 3D image obtained 161

for each of the reconstructed surfaces and for the corresponding ground-truths. The inaccuracies visible on 162 the error maps (see Figure 3D ) do not significantly impact these projections, which appear very similar to 163 the projections obtained with the corresponding ground-truths. Thus, while the biological sample contains 164 significant noise and shows a much more variable contrast (especially with the presence of high intensity 165 globular structures near surfaces S1 and S2), Zellige makes it possible to segment these surfaces selectively, 166

with a quality of segmentation comparable to that obtained by manual expert segmentation. 167

This possibility brings several perspectives that are not offered by single surface extraction 168 algorithms. First it opens the possibility to systematically study the tissue axial movements (along z) relative 169

to the cuticle during molting, allowing for example to gain insights into the early tissue contraction of the 170 wing hinge that acts as a mechanical inducer over the wing blade [8]. Second, Zellige makes it possible to 171

automatically extract structures such as the abdomen epithelium, which is usually segmented manually 172

[34], due to the difficulty to separate the large larval cells from the cuticle mesh and from other globular 173 structures (such as fat bodies or macrophages) present underneath the epithelium. All these structures 174 become intertwined when using other extracting tools. In this context, Zellige opens new opportunities to 175 study collective cell behavior during epithelial morphogenesis in vivo, and to integrate in the analysis the 176 surrounding surface-like structures involved in the mechanics of the system. 177

Example 2: Extracting a thin cochlear epithelium surface from a multilayer dataset. tightly controlled patterning processes during which the cochlear sensory epithelium extends and develops 182 its characteristic coiled snail shape, while adopting a striking cellular mosaic organization, with graded 183 changes of morphogenetic parameters along the cochlea [38, 39] . These morphogenetic processes are well 184 recapitulated in organotypic cultures, on the condition that the mesenchyme that underlies the epithelium 185 6 be preserved. The cultures are then amenable to live imaging [37, 40], pharmacological [41] and genetic 186

manipulations. 187 Figure 4A shows a confocal swept field microscope acquisition of an embryonic mouse cochlea 188 [42] . The sample contains only one surface of interest, the cochlear epithelium, but this surface lays on top 189 of a thick tangled mesh of non-epithelial cells originating from the mesenchyme. The whole biological tissue 190 is stained for filamentous actin (F-actin) using phalloidin. The epithelium surface presents a non-uniform 191 signal included in a small z-range (6 ≤ z ≤ 10), and a mesh of very heterogeneous size. Between sections 192 z = 10 and z = 14 one can observe the basolateral region of the epithelial cells, also stained for F-actin. he 193 particularity of this sample is that the mesenchyme presents an intense and contrasted signal over a wide 194 range of z-values (14 ≤ z ≤ 43). This makes it challenging to extract the surface of the epithelium, which is 195 characterized by low intensity and low contrast. 196 surface is also reconstructed with an excellent coverage (> 99%, Supplementary Table S1). 203

As the sample contains a single epithelial surface of interest, we compared the performance of 204

Zellige with three other software that can extract only a single surface ( Figure 4D ). The projections of 205

PreMosa, FastSME and LocalZProjector completely miss the epithelium. Only regions of high contrast 206

corresponding to the mesenchyme are projected. In contrast, Zellige generates a projection very close to 207 the ground truth. This demonstrates the efficiency of Zellige to selectively extract a low contrast surface, 208

despite the presence of several structures of higher contrast. Indeed, Zellige detects every structure as a 209 possible surface seed without any assumption on its contrast, and only extends this seed into a surface if 210 enough spatial continuity is found in the surrounding signal. This feature allows to separate individual 211 surfaces from other structures spatially, which should greatly facilitate the analysis of live imaging 212 experiments. 213

Example 3: Extracting a single bronchial epithelial surface rendered abnormally rough by SARS-CoV-2. 214

Recently, we used Zellige to extract the surface of a primary culture of bronchial epithelial cells following 215 infection by the SARS-CoV-2 virus [43] . The infection causes the surface of the epithelium to become 216

abnormally rough due to cell damages as seen from discontinuities within the cell layer. The sample we 217 chose from this study is a 3D confocal image of the epithelium responding to SARS-CoV-2 infection [44] 218

( Figure 5A ). The surface of interest in this image corresponds to the layer of epithelial cells stained for the 219 tight junction protein Zona Occludens-1 (ZO-1). The surface roughness causes the network of junctions to 220 extend over the height of the z-stack, with a signal of varying intensity ( Figure 5A ). In addition, the junctional 221 network remains non-planar even at the level of a single cell, hence violating the smoothness condition 222 commonly assumed to hold in the context of epithelial surface extraction. We also observe the presence of 223 nearby punctiform structures of high contrast that are mainly located outside of the epithelium surface. 224

This sample therefore provides an example of a surface with a complex landscape, interspersed with a 225

constellation of signals which may interfere with the surface segmentation. 226

The 3D representation of the reconstructed and corresponding GT height-maps ( Figure 5B ) makes 227 it possible to appreciate the roughness of the surface of interest. The two height-maps overlap quite 228

satisfactorily. As shown on the error map ( Figure 5C ), a large majority (71.1%) of pixels of the reconstructed 229 height-map show errors smaller than 1 pixel (~ 96% of them showing errors smaller than 2 pixels). The error 230

is however larger in regions where the cell size is larger, as well as in areas where the ZO-1 signal intensity 231 is very low, preventing a complete reconstruction of the junctions. Nevertheless, the overall RMSE remains 232 small, with a value of 0.81 (Figure 5) . The coverage of the reconstructed surface is also excellent (98% of 233 the GT height-map), despite the above-mentioned discontinuities. Figure 5D shows the comparison of the 234 projection generated by Zellige to those produced by PreMosa, FastSME, and LocalZProjector. The 235

projection generated by PreMosa misses many junctions of the epithelial network, but it removes quite 236 well the punctiform signal originating from other optical sections. FastSME performs better than PreMosa 237 in reconstructing the junctions, but they produce a projection where the punctiform signal remains strong. 238

In contrast, Zellige and LocalZProjector manage to both reconstruct the surface well and to filter out the 239 punctiform signal quite effectively. This result demonstrates the efficiency of Zellige to extract a surface 240 with complex topography by excluding intense and contrasted spurious signals away from the epithelium 241

surface. 242 apical junctional network of the epithelium is difficult to visualize in microscopy images as it is seen from 252 below, through the basal layer. Another difficulty is the spherical geometry of the vesicle system, which 253 makes the epithelial surface of interest difficult to extract in regions of high inclination relative to the focal 254

plane. 255 Figure 6 shows the result of applying Zellige on a 3D confocal microscopy image of half of a 256 developing inner ear organoid at 14 days of culture, a stage at which markers characteristic of the mouse 257 otic vesicle can be detected [46] . The sample was fixed and stained for F-actin to visualize all cellular 258 structures including the epithelium. Two surfaces of interest can be identified, namely the basal side of the 259 epithelium and the apical junctional network (Figure 6A-B) . Both surfaces are mesh-like structures of high 260 inclination, high signal intensity and high contrast. The vesicle lumen also contains cell debris of high 261 intensity and contrast that are not part of any surface of interest. 262 Figure 6C shows a 3D graphical representation of the height-maps reconstructed by Zellige and 263 those reconstructed by an expert biologist, taken as ground truth (GT). Due to their dome-shaped 264 topography, the manual segmentation of these surfaces was rather laborious, and is more likely prone to 265 errors in the regions of high inclination. Despite this, we observe an excellent match between the two 266 reconstructed and corresponding GT height-maps ( Figure 6D ). The distance between the two height maps 267 is <2 for the large majority (> 92%) of pixels, while larger error values occur locally in regions of near-vertical 268

inclination of the surfaces. The RMSEs of both apical and basal surfaces are close to 1, showing that on 269 average the reconstructed height-maps match the corresponding GT with about pixel accuracy. The 270

coverage of the reconstructed height-map is nearly 100% for the basal surface, and ~88% for the apical 271

surface. Figure 6E shows the 2D projections of the 3D image obtained for each of the reconstructed 272 surfaces, as compared to the projected GT height maps. Thus, for the extraction for these highly inclined 273 surfaces, Zellige produces height maps of a quality comparable to those obtained from manual expert 274 segmentation. 275

In this example, Zellige could be combined with a 2D cell tracking framework such as TissueMiner 276

[3] to perform cell dynamics analysis. Note that geometric distortions introduced in the projected surfaces 277 by the epithelium inclination could be corrected for, using complementary tools such as DProj [11] . This 278 approach could provide a means to quantitatively address how an inner ear organoid epithelium patterns 279

at the cellular and organoid scales, while quantifying the epithelial thickness changes due to cellular 280 intercalation or cell shape changes in the depth of the epithelium. This would also permit to better 281 characterize the variability of inner ear organoids within in a given aggregate, and it could allow one to 282 explore how the organoid interacts with surrounding tissues and how these interactions influence the 283 differentiation of their constituent sensory cells. Table S1 ). The RMSE and coverage of each of the 290 reconstructed surfaces were evaluated and plotted as a function of the value of the parameter that was 291

varied. 292 Figure S3 shows the results of the sensitivity analysis carried out on the image of surfaces S1 and S2, and > 85% for surfaces S3 and S4. The surfaces S1 and S3 show lower signal intensity 300 and lower contrast than surfaces S2 and S3, making them more difficult to extract. Surface S4 has the lowest 301 contrast, and fails to be reconstructed if the classification threshold values are too stringent (namely for 302

Totsu > 12, or TA > 8). Nevertheless, the intervals of stability of TA and Totsu (that is, the intervals of values 303 over which a high-quality extraction of all surfaces is obtained) remain relatively wide (cf. Figure 7) . The 304

smoothing parameters σxy and σz also have some effect on the reconstruction of the surfaces. When σxy is 305 less than 3, the RMSE is higher for the mesh-like epithelial surfaces S2 and S4 (formed by the junctional 306 network of the epithelium). A minimal smoothing along the axial direction is also important to ensure that 307 the reconstructed surfaces are not too fragmented, preventing their complete reconstruction. Yet, σz 308

should not be chosen too large either, to avoid merging nearby surfaces along the z axis. In this case, the 309 closely positioned surfaces S1 and S2 are well-separated if setting σz close to 1, but they become merged 310

when σz > 2. In general, the surface construction parameters have little effect on RMSE and coverage 311

( Figure S3B) . The sensitivity analysis on this challenging specimen shows a good robustness of Zellige to 312 extract the four surfaces with a single set of parameters, each of which can be chosen in a reasonably wide 313

interval considering the other fixed. 314

The results of the sensitivity analyses performed with the other biological image stacks 315 (examples 2, 3 and 4 described above) are shown in Figures S4, S5 and S6 (Supplemental note 4) . These 316 results are summarized in Figure 7 , which show the stability intervals over which the extracted height-maps 317 satisfy the criteria RMSE ≤ 1.5 and coverage ≥ 85%, for which the reconstruction can be considered of high 318 quality. Overall, the stability intervals for the two classification threshold parameters (TA and Totsu) are 319 narrower for specimens containing a surface of low signal intensity and low contrast, while still covering 320 about 1/4 of their respective width. The graphical user interface of Zellige allows the user to adjust TA and 321

Totsu interactively, making it intuitive to search for reasonable values. We found that 2 ≤ σxy ≤ 3 and σz = 1 322

generally give high quality results for all tested specimens ( Figure 7A) . We therefore expect only little 323 adjustment to be required by the user on the smoothing parameters from their default permissive values 324 (set to σxy = 1 and σz = 1). Values of σz that are too large may lead to the merging of a surfaces with a nearby 325 structure of high contrast (surface or else), as it happens for the epithelium surface in the cochlea specimen, 326

which merges with the underlying mesenchyme signal when σz is greater than 2 (Figure 4 and Figure S4 ). 327

The effect is even more pronounced for surfaces of high inclination (Figures 2, 6 and Figures S2, S6) , or 328 presenting a particularly rough texture (Figures 5 and Figures S5) . Finally, the computation times for running Zellige on a given dataset ranged between a few seconds 337 and a minute on a standard PC computer (see Figure S7 and Supplementary note 5), except in a few 338 exceptional cases corresponding to extreme values of the control parameters. As a safeguard, a stopping 339 criterion could be implemented so as to exit the run (declaring the current parameter values invalid) if the 340 surface assembly computation exceeds a user-prescribed duration. 341

Overall, the sensitivity analysis indicates that the surface extraction performed by Zellige is robust 342 to variations of the control parameters of step 1 (surface pixel selection step). In general, the reconstruction 343 is more sensitive to the amplitude threshold parameter (TA), and the Otsu threshold parameter (Totsu)  344 should be kept sufficiently low for samples containing surfaces of intensity close to the background. 345

However, in some cases such as in our example 3, the opposite is true. Thus, the two threshold parameters 346 play somewhat complementary roles, and the possibility to adjust them independently is useful in practice 347

to be able to cover as many cases as possible. A smoothing along xy appears necessary to correctly 348 reconstruct the surfaces supported by a junctional mesh. Not surprisingly, best results are obtained when 349 the radius of the gaussian filter used for this (parameter σxy) is adapted to the mesh (or cell) size. Likewise, 350 a smoothing along z is beneficial, but the extent of this smoothing (parameter σz) should not be too large 351

to avoid causing the fusion of nearby surfaces. With a few exceptions, the values of the RMSE and coverage 352

show little sensitivity to the values of the parameters controlling step 2 (surface assembly step), at least 353 once putative surface pixels have been properly selected. In the presence of several surfaces of potentially 354 very different sizes, the parameter controlling the fraction of OSE sizes allowed for OSE seeds (parameter 355 TOSE1) should be relatively large (≥ 0.5 or greater, i.e. allowing more than 50% of the largest OSE sizes for 356 seeds) to allow the extraction of a surface of small size (for example, to extract the surface n° 4 of 357 example 1, which covers less than 20% of the xy-field of view, TOSE1 must be larger than 0.6). Extreme values 358 (close to 1) for the connectivity rates (C1 and C2) are too stringent and lead to a drop in the coverage of the 359 reconstructed surfaces. To sum up, we see that the most critical parameters for a satisfactory extraction of 360 the different surfaces are those controlling step 1. In most cases the parameters controlling step 2 do not 361 need to be adjusted and can be fixed to their default reference values. 362 363 Conclusion. 364 We have developed Zellige, a new tool to extract multiple surfaces from 3D fluorescence microscopy 365

images. Zellige automatically finds surfaces by first identifying putative pixels that are likely to belong to a 366 biological surface, and second by assembling a surface through connection of adjacent pixels satisfying 367 natural proximity constraints. By using Zellige on synthetic epithelium images we have shown that it 368

accurately reconstructs a surface with excellent performances in terms of both the distance to the ground-369 truth height-map and the surface coverage (Figure 2) . Zellige can deal with complex images containing 370 multiple surfaces, with computation times not exceeding a few tens of seconds on a standard computer. 371

Importantly, the user is not required to specify the number of surfaces to be extracted. In the Drosophila 372 specimen (Figure 3) , the software readily extracts the 4 surfaces of interest that could be identified. Since 373

Zellige detects putative surface pixels in the first step by combining local and global thresholds, it can deal 374 with images where the multiple surfaces display different features, such as in the mouse cochlear embryo 375 (Figure 3) . With this difficult dataset, we could also confirm Zellige's robustness against very low signal-to-376 noise levels. The constructive approach of surface region growing used by Zellige in its second step enables 377 it to circumvent the surface smoothness requirement, that is classically assumed by other surface extraction 378 tools. For instance, it could reconstruct the highly irregular surface of a bronchial tissue infected by SARS-379

CoV-2 ( Figure 5) . 380

The robustness and flexibility of Zellige come at a price, namely, the requirement to specify 12 381 parameters when running the surface extraction. However, the sensitivity analysis we performed shows 382 that adjusting only 4 of these parameters is enough in practice to handle a wide range of image types. These 383 parameters correspond to intuitive notions (e.g., thresholding and smoothing), which makes Zellige 384 particularly easy to use. The Fiji interface that we implemented to perform this adjustment should make 385

Zellige even more user-friendly and effective for biological applications. 386

To our knowledge, Zellige is the only open-source tool that can extract an unspecified number of 387 epithelial surfaces from a 3D volume, possibly larger than two. This unique feature is especially useful in 388 complex images that could be processed only by specialized tools before. For instance, Zellige can extract 389 surfaces with projections on the xy plane that completely overlap, such as the basal and apical epithelia in 390 the organoid image of Figure 6 . Previously, such an image could be processed only by tools that relied on 391 segmenting a mesh around the object surface, such as MorphoGraphX or ImSAnE [16, 47] . 392

The flexibility and robustness of Zellige should allow to considerably relax the constraints that were 393

previously imposed on the sample preparation and the image acquisition steps by the subsequent analysis. In the first step, or surface pixel selection step, a segmentation is applied to the 3D image to select 420 pixels that likely belong to a surface of interest (Figure 1 and Figure S1A ). These putative surface pixels are 421 detected as local maxima of image intensity along the z-axis, after using two independent binary classifiers, 422 one based on pixel contrast and the other one on pixel intensity. Five adjustable parameters control the 423 selection step: two threshold parameters (TA and Totsu) control the strength of the binary classifiers applied 424 on contrast and intensity, respectively, and three parameters (Smin, σxy and σz) control clean-up operations 425 applied at the end of the classification (removal of small isolated spots, and local averaging along the xy 426 plane and the z axis, respectively). 427

In the second step, or surface assembly step, an iterative algorithm is used to extract the height-428 maps of each of the surfaces present in the image (Figure 1 and Figure S1A ). The assembly starts by 429 grouping neighboring putative surface pixels together within each orthogonal (xz or yz) section of the 3D 430 image, in order to form a set of building blocks referred to as orthogonal surface elements (OSEs). These 431 building blocks are then used to assemble the surfaces, in a process analogous to jigsaw puzzles, where 432

OSEs adjacent to the surface boundary are added if they match this boundary, and rejected otherwise, until 433 no matching OSE can be found ( Figure S1 ). In order to increase the robustness of the assembly step, Zellige 434 applies it in two rounds, proceeding along different axes during the first and second rounds. Each round is 435 controlled by 3 adjustable parameters: a threshold parameter (0 ≤ TOSE ≤1) sets a minimum size for the 436 building blocks that can be used as seeds to initiate the assembly of a surface; and two other parameters 437 (0 ≤ R ≤ 50 and 0 ≤ C ≤ 1) set the matching constraints used to accept or reject the addition of OSEs to a 438

surface. The assembly step is thus controlled overall by 6 parameters, i.e. two groups of 3 parameters (TOSE1, 439 R1, C1) and (TOSE2, R2, C2) controlling the first and second assembly rounds, respectively. 440

Finally, the height-maps of each of the reconstructed surfaces are used to obtain a corresponding 441 2D projection (Figure 1) . In practice a maximum projection restricted to a subvolume of width δz (δz being 442 a user-defined parameter) centered around the surface of interest is performed. 443

Availability of data and code. Human bronchial epithelium imaging.

The data used in Figure 5 were taken from the recent study [43] to which we refer for the preparation and 468

imaging of human bronchial epithelium cultures. In brief, MucilAirTM were purchased from Epithelix (Saint-469 Julien-en-Genevois, France) and cultured for at least 4 weeks to reconstruct a differentiated human 470 bronchial epithelium in vitro and stained as previously described. Drosophila imaging.

Flies were raised at 25°C under standard conditions. Pupae were collected for imaging as described 476

previously [49] . Ecad::GFP flies [50] were used for live imaging as previously described [1] . In brief, images 477

were acquired with a spinning disk microscope from Gataca Systems driven by the MetaMorph software. 478

The system is equipped with an inverted Nikon TI2E stand, a motorized XYZ stage, and a Nikon Plan Apo 60x 479 oil immersion (NA=1.4) lens and with a Prime95B camera. 480

Cochlea imaging.

The inner ears from wild-type (C57BL/6) mice were rapidly dissected from temporal bones at embryonic 482 stages E14.5 in HEPES-buffered (10 mM, pH 7.4) Hanks' balanced salts solution and fixed in 4% 483 paraformaldehyde, 1 hour at room temperature. Specimen were permeabilized and stained for phalloidin-484

Atto 565 (Sigma) as previously described [39] . Fluorescence images were obtained with a swept-field 485 confocal microscope (Opterra II) from Brucker. This system is equipped with a Nikon Plan Fluor 60x oil 486 immersion lens (NA=1.4) . 487

Biologists in this study hold a designer certificate of animal experimentation (level 1), allowing them 488 to perform experimental work on animals in strict accordance with the European directive 2010/60/EU, and 489

French regulations. The Ethics Committee of the Institut Pasteur (Comité d'Ethique en Experimentation 490

Animale -CETEA) has approved this study with the project identifier dha170006. This approval is based on 491 careful compliance to the 3Rs principle in the care and use of animals (Annex IV -2010/60/EU). 492

Inner ear organoid imaging. Triton X-100, for 3-5 days to clarify the tissue. Finally, the aggregates were whole-mounted using the ScaleS 502 solution and imaged using a confocal laser scanning microscope (A1R HD25, Nikon) equipped with a Nikon 503 25x silicon oil immersion lens (NA=1.05). 504

Acknowledgements. 505 We thank Maia Brünstein of the Hearing Institute Bioimaging Core Facility -C2RT/C2RA and Florian Rückerl 506 of the Photonic BioImaging -C2RT/C2RA, for sharing their expertise on light microscopy. We thank Romain 507

Levayer for sharing flies and lab space with us. We thank Gilles Trébeau (https://fr.linkedin.com/in/gilles-508 trebeau-427601212) for helpful discussions regarding the Java frontend deployment of Zellige. We thank 509

Lisa Chakrabarti, Rémy Robinot and Vincent Michel for sharing data. This work was supported by Step 1

Step 2 n°1 n°2 n°3 n°4 n°1 n°3 n°2 n°4 in the dataset (of dimensions 1200 × 1200 × 51 pixels): surfaces S1 and S2 are relatively close to one another and located within overlapping z-ranges (8 ≤ z ≤ 50 and 20 ≤ z ≤ 50, respectively). Surfaces S3 and S4

(located in the z-ranges 42 ≤ z ≤ 50 and 9 ≤ z ≤ 50, respectively) are relatively far from each other and can nearly be separated by a plane. (C) 3D representations of the height maps extracted by Zellige (in green) and

of the ground truth height maps (GT, in blue) of surfaces S1-S4. The reconstructed height-maps of all surfaces S1-S4 cover >93% of the area of the corresponding GT (cf. Figure S2 and Table S1 ). To reduce the staircase artifacts (more or less visible depending on the surface) due to the digitization of the GT and reconstructed height-maps, all height-maps were smoothed with a 2D gaussian filter with a standard radius of 5 pixels (cf.

Supplemental note 1). (D) Error maps (color-coded distance along the z-axis between the reconstructed and the GT height-maps) plotted for each of the reconstructed surfaces. The large majority of pixels on the reconstructed height-maps (98%, 96%, 91%, and 99% for surfaces S1 to S4, respectively) display errors of <2 pixels. The height-maps of surfaces S1, S2, S4 show subpixel accuracy on average (RMSE < 1), while that of surface S3 is slightly less accurate (RMSE = 1.25). (E) Projections of the 3D image localized to the different surfaces S1-S4 (in this and all subsequent figures, these are maximum intensity projections over a subvolume of width δz=±1 pixel above or below the corresponding height-maps). Upper and lower panels show the projections based on the GT and the reconstructed height-maps, respectively. Scale bar 50 µm. 

Interplay of cell dynamics and 515 epithelial tension during morphogenesis of the Drosophila pupal wing

Unified quantitative 518 characterization of epithelial tissue development

Cellular growth and rearrangement during the development of the 607 mammalian organ of Corti

Cochlear outer 609 hair cells undergo an apical circumference remodeling constrained by the hair bundle shape

Cell migration, intercalation and growth regulate mammalian 612 cochlear extension

Wnt5a functions in planar cell polarity 614 regulation in mice

Zellige example dataset: mouse 616 embryonic cochlea stained for F-actin

SARS-CoV-2 infection induces 618 the dedifferentiation of multiciliated cells and impairs mucociliary clearance

Zellige example 621 dataset: primary culture of human bronchial cells infected by SARS-CoV-2. 2022

Generation of inner ear sensory epithelia from 624 pluripotent stem cells in 3D culture

Zellige example dataset: inner ear 626 organoid epithelium

MorphoGraphX: A platform for quantifying morphogenesis in 4D

ImgLib2-generic image processing in Java

Imaging Drosophila Pupal Wing Morphogenesis

Methods and Protocols

Directed, efficient, and versatile modifications of the 636 Drosophila genome by genomic engineering

Directed Differentiation of Mouse Embryonic Stem Cells Into Inner Ear 639

Sensory Epithelia in 3D Culture

Flowchart of Zellige's algorithmic steps. Surface pixel selection (step 1), surface assembly in the form of a height map (step 2), and subsequent projection localized to the height map, are schematically depicted in the case of a 3D image containing 4 surfaces of interest.

Multiple surface extraction on a synthetic 3D image. (A) The image contains 3 phantom surfaces (S1, S2, S3) of different shapes (sinusoidal, flat, and paraboloidal, respectively), and different textures (surface S1 has constant intensity, while surfaces S2 and S3 are supported by Voronoi meshes of different cell-sizes). (B) 3D representations of the height maps extracted by Zellige (in green) and of the ground truth (GT, in blue) height maps of surfaces S1, S2, and S3. (C) Error maps displaying the distance along the z-axis between the reconstructed and GT height maps for surfaces S1, S2, and S3. (E) Projections of the 3D image localized to the different surfaces S1-S3 (maximum intensity projections over a subvolume of a width z=1 pixel above or below the corresponding height-maps). Upper and lower panels show the projections based on the GT and the reconstructed height-maps, respectively. , and a portion of the developing wing. Scale bar 50 µm. Four surfaces of interest may be identified in the dataset (of dimensions 1200 × 1200 × 51 pixels): surfaces S1 and S2 are relatively close to one another and located within overlapping z-ranges (8 ≤ z ≤ 50 and 20 ≤ z ≤ 50, respectively). Surfaces S3 and S4 (located in the z-ranges 42 ≤ z ≤ 50 and 9 ≤ z ≤ 50, respectively) are relatively far from each other and can nearly be separated by a plane. (C) 3D representations of the height maps extracted by Zellige (in green) and of the ground truth height maps (GT, in blue) of surfaces S1-S4. The reconstructed height-maps of all surfaces S1-S4 cover >93% of the area of the corresponding GT (cf. Figure S2 and Table S1 ). To reduce the staircase artifacts (more or less visible depending on the surface) due to the digitization of the GT and reconstructed height-maps, all height-maps were smoothed with a 2D gaussian filter with a standard radius of 5 pixels (cf. Supplemental note 1). (D) Error maps (color-coded distance along the z-axis between the reconstructed and the GT height-maps) plotted for each of the reconstructed surfaces. The large majority of pixels on the reconstructed height-maps (98%, 96%, 91%, and 99% for surfaces S1 to S4, respectively) display errors of <2 pixels. The height-maps of surfaces S1, S2, S4 show subpixel accuracy on average (RMSE < 1), while that of surface S3 is slightly less accurate (RMSE = 1.25). (E) Projections of the 3D image localized to the different surfaces S1-S4 (in this and all subsequent figures, these are maximum intensity projections over a subvolume of width z= 1 pixel above or below the corresponding height-maps). Upper and lower panels show the projections based on the GT and the reconstructed height-maps, respectively. (D) Color-coded error maps of the reconstructed height-maps for the apical (left) and basal (right) epithelial surfaces of the organoid. The surfaces of interest are reconstructed with an error of < 2 pixels over a large majority (96% and 93% for the apical and basal surfaces, respectively) of pixels, as well as on average (RMSE ~ 0.8 and 1.1 for the apical and the basal surfaces, respectively).(E) Projections localized to the GT height-maps of the epithelium surface (panels on the left), and the height-maps extracted by Zellige (panels on the right). Scale bar 100 µm.

Summary of the sensitivity analysis. The intervals indicated in grey for each parameter and each of the images tested correspond to the parameter values for which the reconstruction satisfies high quality criteria defined by RMSE ≤ 1.5 and coverage ≥ 85%. Black marks indicate the reference value obtained by manual adjustment for each image (cf. Supplemental note 2). (A) Parameters of the surface selection step. (B) Parameters of the surface assembly step.