CartoCell, a high-content pipeline for 3D image analysis, unveils cell morphology patterns in epithelia

Article CartoCell, a high-content pipeline for 3D image analysis, unveils cell morphology patterns in epithelia Graphical abstract Highlights dCartoCell exploits a small training dataset for accurate highcontent 3D segmentation dCartoCell enables quantification of morphological features at the cellular level dUsing CartoCell, we generate a ‘‘cartography’’ of morphological patterns in 3D dCartoCell applies to diverse epithelial models Authors Jesu ´s A. Andre ´s-San Roma ´n, Carmen Gordillo-Va ´zquez, Daniel Franco-Barranco, ..., Pedro Go ´mez-Ga ´lvez, Ignacio Arganda-Carreras, Luis M. Escudero Correspondence [email protected].ac.uk (P.G.-G.), ignacio.argand[email protected] (I.A.-C.), [email protected] (L.M.E.) In brief Analyzing 3D epithelia at the cell level becomes challenging when dealing with hundreds of samples without a training dataset. Andre ´s-San Roma ´netal. introduce CartoCell, a deep-learning pipeline that employs a small annotated dataset to achieve high-content cyst segmentation coupled with detailed analysis and cartography of cell features. Andre ´s-San Roma ´n et al., 2023, Cell Reports Methods 3, 100597 October 23, 2023 ª2023 The Authors. https://doi.org/10.1016/j.crmeth.2023.100597 ll Article CartoCell, a high-content pipeline for 3D image analysis, unveils cell morphology patterns in epithelia Jesu ´s A. Andre ´s-San Roma ´n, 1,12 Carmen Gordillo-Va ´zquez, 1,12 Daniel Franco-Barranco, 2,3,12 Laura Morato, 1 Cecilia H. Ferna ´ndez-Espartero, 1 Gabriel Baonza, 4 Antonio Tagua, 1 Pablo Vicente-Munuera, 5 Ana M. Palacios, 1 Marı ´a P. Gavila ´n, 6 Fernando Martı ´n-Belmonte, 4 Valentina Annese, 1 Pedro Go ´mez-Ga ´lvez, 1,7,8, * Ignacio Arganda-Carreras, 2,3,9,10, *and Luis M. Escudero 1,11,13, * 1 Instituto de Biomedicina de Sevilla (IBiS), Hospital Universitario Virgen del Rocı ´o/CSIC/Universidad de Sevilla and Departamento de Biologı ´a Celular, Facultad de Biologı ´a, Universidad de Sevilla, 41013 Seville, Spain 2 Department of Computer Science and Artificial Intelligence, University of the Basque Country (UPV/EHU), 20018 San Sebastian, Spain 3 Donostia International Physics Center (DIPC), 20018 San Sebastian, Spain 4 Program of Tissue and Organ Homeostasis, Centro de Biologı ´a Molecular Severo Ochoa, CSIC-UAM and Ramo ´n & Cajal Health Research Institute (IRYCIS), Hospital Universitario Ramo ´n y Cajal, 28034 Madrid, Spain 5 Laboratory for Molecular Cell Biology, University College London, London, UK 6 Centro Andaluz de Biologı ´a Molecular y Medicina Regenerativa (CABIMER), JA/CSIC/Universidad de Sevilla/Universidad Pablo de Olavide and Departamento de Citologı ´a e Histologı ´a Normal y Patolo ´gica, Facultad de Medicina, Universidad de Sevilla, 41009 Seville, Spain 7 MRC Laboratory of Molecular Biology, Cambridge Biomedical Campus, Francis Crick Avenue, Trumpington, Cambridge CB2 0QH, UK 8 Department of Physiology, Development and Neuroscience, University of Cambridge, Cambridge CB2 3EG, UK 9 Ikerbasque, Basque Foundation for Science, 48009 Bilbao, Spain 10 Biofisika Institute, 48940 Leioa, Spain 11 Biomedical Network Research Centre on Neurodegenerative Diseases (CIBERNED), 28029 Madrid, Spain 12 These authors contributed equally 13 Lead contact *Correspondence: [email protected] (P.G.-G.), [email protected] (I.A.-C.), lmescudero[email protected] (L.M.E.) https://doi.org/10.1016/j.crmeth.2023.100597 SUMMARY Decades of research have not yet fully explained the mechanisms of epithelial self-organization and 3D packing. Single-cell analysis of large 3D epithelial libraries is crucial for understanding the assembly and function of whole tissues. Combining 3D epithelial imaging with advanced deep-learning segmentation methods is essential for enabling this high-content analysis. We introduce CartoCell, a deep-learning-based pipeline that uses small datasets to generate accurate labels for hundreds of whole 3D epithelial cysts. Our method detects the realistic morphology of epithelial cells and their contacts in the 3D structure of the tissue. CartoCell enables the quantification of geometric and packing features at the cellular level. Our single-cell cartography approach then maps the distribution of these features on 2D plots and 3D surface maps, revealing cell morphology patterns in epithelial cysts. Additionally, we show that CartoCell can be adapted to other types of epithelial tissues. INTRODUCTION Analysis of epithelial tissue properties at the cellular level has enabled advancement in the understanding of different cellular phenomena during morphogenesis. Traditionally, most approaches were based on two-dimensional (2D) analysis of the apical surfaces of monolayer epithelia. However, the need to understand the 3D morphology of epithelial cells to study MOTIVATION A major bottleneck in developing neural networks for cell segmentation is the need for laborintensive manual curation to develop a training dataset. The present work addresses this limitation by developing an automated image-analysis pipeline that utilizes small datasets to generate accurate labels of cells in complex 3D epithelial contexts. The overall goal is to provide an automatic and feasible method to achieve high-quality epithelial reconstructions and to enable high-content analysis of morphological features, which can improve our understanding of how these tissues self-organize. Cell Reports Methods 3, 100597, October 23, 2023 ª2023 The Authors. 1 This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). ll OPEN ACCESS organogenesis, 1–4 cell migration, 5 branching formation, 6 tumorigenesis, 7 or wound healing 8 has become evident in recent years. A major breakthrough was the discovery that epithelial cells can present very complex geometries due to the exchange of neighbors along the apico-basal axis. These cell shapes have been called scutoids, and it has been shown that they have a role in morphogenesis as well as in the connectivity and biophysical properties of tissues, cushioning and minimizing cell surface tension and leading to a balanced energetic state. 9–16 Scutoids represent a new paradigm, but they set a challenge for the quantitative analysis of the complex epithelial 3D packing requiring a very accurate reconstruction of 3D epithelial tissues from microscopy images to allow the capture of precise cell shapes and neighboring relationships. In the last few years, deep learning has become the state-ofthe-art solution for the analysis of biomedical images. 17–19 Deep learning is a subdomain of machine learning that makes use of large (or so-called deep) artificial neural networks to solve a wide variety of tasks. As opposed to conventional algorithms, before they can be used, deep-learning methods (or models) need to be ‘‘trained.’’ In other words, models can ‘‘learn’’ from a set of examples how to solve a specific task. Once trained, the models can be directly applied to new samples, usually called prediction or ‘‘inference.’’ In the particular case of image segmentation, the training dataset is commonly formed by a set of raw images and their corresponding ground-truth annotations or ‘‘labels.’’ This type of learning framework, with both raw and label images available, is known as supervised learning. Furthermore, realistic 3D reconstruction of epithelial cells requires assigning each individual cell a unique label, in a process called ‘‘instance segmentation.’’ The 3D instance segmentation of microscopy data is a difficult task, especially in the presence of a dense concentration of cells and anisotropic voxel resolution, as it is common in volumetric images of epithelial tissue. State-of-the-art learning-based methods tackle these challenges using a top-down strategy, by first training a deep neural network (DNN) to predict representations of the objects of interest (cells in our case) and then extracting individual instances from those representations using different post-processing methods. Common representations include cell masks or boundaries, 20–23 distance or flow maps, 24–26 or a combination of some of the latter. 27,28 On top of those representations, cell instances are then calculated usually by means of watershed 29,30 or graph-partitioning methods. 31–33 Other techniques have shown success in segmenting cell nuclei and tracking cell lineage, 34,35 but they do not have the high level of accuracy in cell shape required to obtain detailed geometric and topological information at the cellular level. Despite the benefits observed from these supervised approaches, their main drawback is the large number of annotated samples needed to establish a training dataset and obtain reliable performance. 36 Preparing and processing such a large amount of data manually or semi-automatically is usually tedious and time expensive. This problem arises from the acquisition time of high-resolution images as well as from the labeling of raw images performed by experts in the field, which is usually the main bottleneck of the protocol. To address this issue, a common strategy consists in using data augmentation, i.e., synthetically increasing the size of the training data by morphological and intensity transformations or noise addition. 37,38 However, data augmentation may not be sufficient to realistically recreate the diversity of image data to be processed. A much less exploited alternative to speed up the segmentation protocol would consist in the use of low-resolution images instead, which are acquired and annotated at a considerably faster pace. Certainly, this option would be ideal if the quality of the output segmented cells remains comparable to that obtained with high-resolution images. In this article, we image, process, and analyze whole MadinDarby canine kidney (MDCK) 3D epithelial cell cultures. 39 Despite their simplicity, these cysts have previously been used as a suitable model system to study the establishment of cell polarity and cell junctions, 3,40–43 epithelial morphogenesis and physiology, 44–51 and tumor progression, 52–55 and for exploring the constraints on epithelial tissue morphogenesis. 56 MDCK cysts have provided valuable insights into the study of more complex systems, helping us to understand self-organization in organoids, embryoids, and the early stages of mammal development. 56–58 In our approach, we use a small training dataset of high-resolution images to subsequently produce a large training dataset of low-resolution images that are automatically segmented. Our method follows a top-down pipeline that makes use of DNN architecture with multiple cell representations and watershed post-processing to initially segment the epithelial cells as 3D instances. These instances are refined by a second post-processing step: a 3D Voronoi algorithm that provides realistic 3D epithelial boundaries where cells are in close contact with each other. The algorithm is based on tiling the space between a set of (Voronoi) seeds by proximity, without leaving any gaps among the generated compartments. 59 These compartments are called Voronoi cells. Honda showed that the Voronoi compartmentalization of a 2D space, after using the cell nuclei as seeds, fitted the pattern of cellular contacts found in epithelial surfaces. 60 In 3D, this approach has been previously used to simulate the shapes of globular cells. 10 Here, it is key to increasing the quality of our cell-segmentation results. In short, we have developed an accessible and fast tool to investigate the complex organization of epithelial tissues. The production of a large number of samples with accurate segmentation has opened up a new way of 3D high-content analysis. The representation of the extracted feature values in each cell provides maps of the cysts at single-cell resolution. Given the similarities with the practices of making and using maps, we called our approach single-cell cartography, and named our high-content segmentation method CartoCell. The simple observation of these maps reveals the presence of cell morphology patterns where cells are distributed following geometric cues. These patterns illustrate how different the cells within the same cyst really are, and how cells with similar characteristics have the tendency to cluster together in specific zones of the cysts. Importantly, the large number of processed individual cells permits us to quantify the frequency of the patterns and even to find hidden traits of organizational features within the 3D structure of the tissue. 2Cell Reports Methods 3, 100597, October 23, 2023 Article ll OPEN ACCESS RESULTS CartoCell, a high-throughput pipeline for segmentation of 3D epithelial cysts The realistic analysis of whole epithelial tissues at the cell level is a critical point to a bottom-up understanding of how tissues self-organize during development. In this work, by means of deep-learning and image-processing strategies, we have developed CartoCell, an automated pipeline (Figure 1) to segment and analyze hundreds of epithelial cysts at different stages (Table S1) with minimal human intervention. CartoCell is subdivided into five consecutive phases. In phase 1, a small dataset of 21 cysts, stained with cell outline markers, was acquired at high resolution in a confocal microscope (Figure 1, 8 cysts of 4 days and 13 cysts of 7 days, STAR Methods). Next, the individual cell instances were segmented using LimeSeg, 61 STAR Methods), a semi-automatic segmentation plugin of Fiji. 62 The final high-resolution label images were the output of a curation process aided by a custom MATLAB code (STAR Methods). In particular, we implemented a MATLAB graphic user interface to facilitate manual deletion, insertion, fusion, and proper profiling of cell instances and lumen segmentations (STAR Methods). The use of high-resolution images simplified and improved the accuracy of manual annotations by providing a clearer visualization of the cysts and their structures. On average, we estimate that the segmentation and curation process took 3–5 complete working days of one person per cyst. The high-resolution images from phase 1 provide the accurate and realistic set of data necessary for the following steps (see discussion). In phase 2, both high-resolution raw and label images were down-sampled to create our initial training dataset (Figure 1). The logic of this step is to take advantage of lower storage requirements and faster acquisition and processing time of low-resolution images. Specifically, the image volumes were reduced to match the resolution of the images acquired in phase 3(STAR Methods). Using that dataset, a first DNN was trained. The DNN employed was a custom stable 3D residual U-Net (3D ResU-Net) 63 (Figures 1 and S1). We will refer to this first model as ‘‘model M1’’ (Figure 1 and STAR Methods). In phase 3, a large number of low-resolution stacks of multiple epithelial cysts was acquired (Figure 1). This was a key step in allowing high-content analysis of samples, since it greatly reduces the acquisition time (STAR Methods). Here, we extracted the single-layer and single-lumen cysts by cropping them from the complete stack (Figure 1 and STAR Methods). In this way we obtained a set of 293 low-resolution images, composed of 84 cysts at 4 days, 113 cysts at 7 days, and 96 cysts at 10 days (Figure S2). Next, we applied our trained model M1 to those images and post-processed their output to produce (1) a prediction of individual cell instances (obtained by marker-controlled watershed) and (2) a prediction of the mask of the full cellular regions (Figures 1 and S1;STAR Methods). At this stage, the output cell instances were generally not touching each other, which is a problem when studying cell connectivity in epithelia. Therefore, we applied a 3D Voronoi algorithm to correctly mimic the epithelial packing. 59,60 More specifically, each prediction of cell instances was used as a Voronoi seed, while the prediction of the mask of the cellular region defined the bounding territory that the cells could occupy (Figures 1 and S3;STAR Methods). In previous works, cell nuclei were used as Voronoi seeds, leading to less reliable cell outlines, since nuclei may not be exactly located at Voronoi centers. 64 In our approach, the predicted full cell instance masks are used as seeds, producing more accurate results, since only the inter-instance space needs to be filled by the algorithm. The output of this phase was a large dataset of lowresolution images and their corresponding accurate labels. In phase 4, a new 3D ResU-Net model (‘‘model M2’’ henceforth) was trained on the ‘‘training low-resolution dataset,’’ composed of the newly produced large dataset of low-resolution raw images and its paired label images, along with the ‘‘training down-sampled dataset’’ (Figure 1 and STAR Methods). This was a crucial step, since the performance of deep-learning models is highly dependent on the number of training samples. In phase 5, model M2 was applied to new low-resolution cysts and their output was post-processed as in phase 3, thus achieving high-content segmentation of the desired cysts (Figures 1 and S2). Optimization of the method Once the CartoCell pipeline was defined, we performed an automatic screening of parameters of the 3D ResU-Net to optimize the quality of the prediction of models M1 and M2 (STAR Methods). To this aim, we elaborated a test set with 60 new low-resolution cysts (20 cysts at 4 days, 20 at 7 days, and 20 at 10 days of development) not used in any of the previous training steps (STAR Methods). These cysts were semi-automatically segmented and manually curated to obtain their groundtruth labels (STAR Methods). The parameter search for the M1 and M2 models aimed to ensure the highest quality of segmentation, based on the comparison between the prediction of the models and the ground-truth labels of the test set (Table S2 and STAR Methods). In addition, this optimization also demonstrated that M2, a model trained with a large number of low-resolution images, outperformed M1, a model trained with a small but perfectly segmented dataset (Table S3). Performance evaluation and comparison with other methods To test CartoCell against current alternatives, we compared the performance of our segmentation pipeline with that provided by the state-of-the-art approaches StarDist 3D, 25 Cellpose, 26 and PlantSeg 23 (Figure S3 and STAR Methods). In all cases, the small down-sampled dataset from phase 2 was used as the training set, and the same 60 new cysts were used as the test set (STAR Methods). Moreover, for the sake of analyzing method robustness and stability, each method was trained ten times under the same conditions. Segmentation metrics were thus provided on average over those ten repetitions (Table S3). In summary, CartoCell compares favorably with state-of-the-art alternatives, especially after retraining our model with the new dataset in phase 4. As demonstrated by the segmentation metric values of model M1 and M2 (Table S3), one of the keys of the success of CartoCell relies on generating a large yet imperfect training dataset, which greatly enhances the prediction accuracy when Cell Reports Methods 3, 100597, October 23, 2023 3 Article ll OPEN ACCESS training a DNN. To measure this improvement, we quantified the number of perfectly annotated cysts that would be required in a single phase to achieve results equivalent to those of the whole pipeline of CartoCell. In particular, we trained our 3D ResU-Net with 50, 100, 150, 200, and 250 perfectly annotated cysts (STAR Methods). As a result, we observed that as many as 90 cysts would be needed to match the performance of CartoCell with only 21 cysts as input (Figure S3). Additionally, we demonstrated that other state-of-the-art DNNs can be used in conjunction with CartoCell. By leveraging on the long imperfect training dataset resulting from phase 2, these methods can be integrated into CartoCell by replacing the model M2 in phase 4 (Figure 1). In line with the findings in our original CartoCell pipeline, the use of the new M2 models improved the segmentation results with respect to those of model M1 (Table S3). Figure 1. CartoCell pipeline for high-content epithelial cysts segmentation Phase 1: the ‘‘high-resolution raw images’’ consist of confocal z-stack images, where the cell membrane is stained. These images are segmented and proofread using LimeSeg and a custom MATLAB code for curation to obtain the ‘‘high-resolution label images’’ (STAR Methods). Together, the raw and the label images encompass the ‘‘training high-resolution dataset.’’ Number of samples = 21. Phase 2: the ‘‘training high-resolution dataset’’ is down-sampled to obtain the ‘‘training down-sampled dataset,’’ which is the training set for the ‘‘model M1.’’ Phase 3: low-resolution images are obtained from confocal z-stack images, stained in a similar way to phase 1. Number of samples = 293. Scale bar, 100 mm. Next, the ‘‘work-flow M1’’ is applied: inference using ‘‘model M1’’ and subsequent post-processing to obtain individual cell instance predictions and cell masks, followed by the 3D Voronoi algorithm to guarantee that predicted cells remain in close contact. As a result, the ‘‘low-resolution label images’’ are generated. Phase 4: training of the ‘‘model M2’’ on the large ‘‘training low-resolution dataset.’’ Number of samples = 314. Phase 5: high-content segmentation of new low-resolution images (unseen by the pipeline) using the ‘‘work-flow M2,’’ which is equivalent to the ‘‘work-flow M1’’ but using the ‘‘model M2.’’ See also Figures S1–S3;Tables S2,S3, and S4. 4Cell Reports Methods 3, 100597, October 23, 2023 Article ll OPEN ACCESS (legend on next page) Cell Reports Methods 3, 100597, October 23, 2023 5 Article ll OPEN ACCESS Management of results: Final curation before biological analysis As a result of this full process, we obtained 353 segmented cysts (293 low-resolution from phase 3 + 60 tests). Next, with the purpose of conducting a detailed quantitative analysis of 3D epithelial packing, we performed a semi-automatic final curation to correct small defects in the segmentation (Figure S3 and STAR Methods). As an outcome, we obtained our final ground-truth dataset, with accurate feature values that we used in the biological analysis (Figure 2;Tables S1 and S4). On average, each imperfect cyst took 12 ±6 min to be curated. Thanks to our ground-truth dataset we could compare the values of the biological features extracted before and after the final curation step, thus measuring their impact on the results. Nevertheless, the user may opt for skipping this semi-manual step and keep a fully automatic processing pipeline. In our specific case, we automatically selected for analysis the cysts released after phase 5 whose epithelial monolayer was completely tiled by cells, i.e., without gaps produced by under-segmentation, namely, 307 ‘‘closed cysts’’ out of the 353 segmented cysts (Figures 1 and S3). The mean relative error of all the geometric features extracted from the closed cysts was 5.6% ±4.1%. In the case of connectivity features, the mean relative error was larger: 13.3% ±13.9%. (Figure S3,Table S4,STAR Methods, and discussion). Epithelial cysts adopt different shapes in 3D culture The pipeline that we have developed is able to realistically reconstruct the whole cyst and its lumen. In this way, it allowed us to identify the real shape of the complete 3D epithelial structure. We found that the full set of processed cysts (our ‘‘ground-truth dataset’’) can present a high heterogeneity in terms of shape (Figure S2). We analyzed the 3D structure of the total number of single-lumen cysts pointing to two geometrical features: axes lengths and solidity (a curvature index) (Figure S4 and STAR Methods). We considered two axes to be similar when the difference in their lengths was inferior to 10%. According to these considerations, we established five types of shapes (Figures 2A and S2;STAR Methods). When the three axes of symmetry were similar in length, cysts were classified as spheres (4.0%). When two of the three axes were similar, but one differed, they were called spheroids. Moreover, spheroids were divided into prolate (18.1%) when the different axis was the major one and oblate when it was the minor one (34.3%). Cysts with three axes of different lengths were classified as ellipsoids (39.9%). Finally, we categorized cysts as presenting negative curvature (3.7%) when solidity was less than 0.9, independently of axes values (Figures 2A and S4;STAR Methods). We also quantified the frequency of each type of shape in 4-, 7-, and 10-day cysts. We found that at 4and 7-day time points the most frequent shape was ellipsoid, with the percentage of oblate cysts increasing with the time of culture (Figure 2B). At 10 days there were 37.9% ellipsoids and 44.8% oblate cysts. Single-cell geometric analysis reveals cell morphology patterns in the MDCK cysts Reconstructing the 3D outlines of all the cells allowed the precise quantification of a large number of geometrical and connectivity characteristics (Table S1 and Figure S4). We designed a strategy to visualize these data in two ways: 3D maps of the surface of the reconstructed cysts and 2D plots that represent the position of all the cells analyzed within the cysts (Figure 2C). We named this approach single-cell cartography. In the case of the 3D maps, we obtained a readout of seven cell geometric features by plotting their normalized values using a color palette (Figures 2D, 2F, and S4;STAR Methods). Regarding the ‘‘cell height’’ there was a clear gradient ‘‘top to bottom,’’ with the cells with lower values on the top of the cyst and a progressive increase of their height toward the base (Figure 2D). Importantly, detailed examination of the entire surface of the cysts revealed Figure 2. Realistic high-content 3D segmentation reveals different morphologies and cell morphology patterns in MDCK cysts (A) Shape classification of single-layer, single-lumen cysts. (Top) The 3D rendering of representative segmented cysts. (Bottom) Schematic representation of the morphological cyst classification. (B) Frequency distribution of the different cyst shapes at 4, 7, and 10 days. (C) Schematic representation of the relation between 3D reconstructions of the cyst and the 2D plots. The drawing shows how the position and angle of the cell centroids is represented. (D and E) Single-cell cartography of ‘‘cell height’’ feature. (D) Computer rendering of three representative segmented cysts: front, back and bottom views. The cell color scale symbolizes the value of the cell height (normalized by cyst). (E) (Left) Polar scatter showing the normalized distance and angle of cell centroids regarding the cyst centroid (considering 0any vector contained in the xy plane of the cyst centroid; positive angles correlate with cell centroids placed over the cyst centroid, and cells located below the cyst centroid xy plane are represented with negative angles, as indicated in C), with a heatmap coloring the normalized value of cell height. 20,391 cells from 353 cysts were represented. (Right) Cell sorting by the normalized cell height values, from left to right: 0–0.25, 0.25–0.5, 0.5– 0.75, 0.75–1. On top is a scatter polar diagram showing the angle and distance of cell centroids regarding the cyst centroid, with a heatmap coloring based on the quantification of the normalized cell height; on the bottom is a polar histogram accounting for the frequency of cell positions. (F and G) Single-cell cartography of the ‘‘apical area’’ feature. The same color code and plotting properties as in (D) and (E) are used. (H–J) Cyst morphology affects the cartography of cell height at 10 days. (H) Cell z position versus cell height (normalized per cyst). The color of the dots is determined by the kind of cyst to which the cells belong. 8,528 cells from 116 cysts at 10 days were analyzed. (I) Scatterplots showing the comparison in cell z position versus cell height in 10-day ellipsoid, oblate, and prolate cysts. (J) Cell proportion in the highlighted quadrant of normalized cell z position <0.15 and normalized cell height <0.5 in 10-day ellipsoid, oblate, and prolate cysts. The percentage is computed with respect to the total number of cells for each morphology and experiment. Data are presented as mean ±SD. *p < 0.05 (Student’s t test). (K–M) Single-cell cartography of scutoids. (K) Computer rendering of three representative segmented cysts in three perspectives, showing scutoids in black and non-scutoid cells in light gray. (L) Histogram representing the cell proportion at different intervals of the normalized xy-position distance with respect to the xy coordinate of the cyst centroid. Percentage of all cells from the 353 segmented cysts. (M) Histogram representing the cell proportion of scutoid and non-scutoid cells at different intervals of the z-position distance with respect to the z coordinate of the cyst centroid. Percentage of all cells from the 353 segmented cysts. See also Figures S4–S6;Tables S1 and S5. 6Cell Reports Methods 3, 100597, October 23, 2023 Article ll OPEN ACCESS Figure 3. CartoCell segmentation on other epithelial tissue datasets (A–D) Hypoxia can induce cell morphology pattern changes in MDCK cysts. (A) Middle sections of top-to-bottom confocal microscopy z-stack images of representative cysts cultured under normoxic (left) or hypoxic (1% O 2 ) (right) conditions at 4 days. Cell contours were stained with Alexa Fluor 647 phalloidin (legend continued on next page) Cell Reports Methods 3, 100597, October 23, 2023 7 Article ll OPEN ACCESS a subpattern: just in the center of the bottom region, the cells were shorter (Figure 2D). These complex cell morphology patterns were consistent among a high number of cysts and appeared at different time points (Figure 2D). To confirm that the pattern was general, we leveraged high-content analysis and visualized the height of each cell of all processed cysts. We plotted the ‘‘cell height’’ data (with the color code used in the 3D cysts) on a scatter polar diagram considering the angle and radius of the cell centroids with respect to the centroid of the whole cyst (Figure 2E). In this way, we were able to visualize the distribution of the ‘‘cell height’’ feature in all cells, from all cysts, at the same time. The plots confirmed the patterns observed in the individual cysts: ‘‘shorter’’ cells (yellower colors) were located at the top and center of the bottom of the cysts, while ‘‘taller’’ cells (pink-purple colors) were found at the periphery of the lateral and bottom part of the cysts. A similar cell morphology pattern, although not as evident, was observed in the distribution of ‘‘cell basal area,’’ ‘‘cell volume,’’ and ‘‘cell surface area’’ values (Figure S5). Furthermore, we found a different pattern involving the distribution of the ‘‘cell apical area’’ values (Figures 2F and 2G). In this case, cells with a bigger apical area were enriched on the top and at the bottom of the cysts. Meanwhile, cells with the smaller apical area were located in the middle region of the cysts. We also found neither ‘‘cell solidity’’ nor ‘‘cell aspect ratio’’ characteristic showed any pattern (Figure S5). The cell morphology pattern of some features can correlate with the shape of the whole cysts Our high-content approach revealed that the MDCK cysts can present different shapes and also intrinsic cell patterns. To test whether the cell morphology patterns can be affected in some way by the global shape of the cysts, we plotted the values of the features against the position of the centroid of the cells in the z axis of the cyst. In this way, we can quantify differences in populations of cells between the three more abundant categories: ellipsoids, spheroids oblate, and spheroids prolate. In the case of ‘‘cell height,’’ there was a clear and robust gradient from ‘‘shorter’’ to ‘‘taller’’ cells, from the top toward the bottom, on the three types of shapes analyzed (Figures 2H and 2I). However, a more detailed analysis of the bottom-left side of the graphs (corresponding to the shorter cells in the base of the cysts) revealed significant differences on the 10-day cysts (Figures 2I and 2J) but not on 4-day and 7-day cysts (Figure S6, Table S5, and STAR Methods). We also obtained differences in the case of the ‘‘basal area’’ feature, but again only on 10-day cysts (Figure S6 and Table S5). Conversely, we did not find differences at any time point with the ‘‘cell apical area’’ feature (Figure S6 and Table S5). Our results suggest that the shape of the whole cyst could correlate with changes in cell morphology patterns (see discussion). The emergence of cell-packing patterns in the cysts Motivated by the finding of cell morphology patterns in the distribution of the values of cell geometric features, we examined the presence of particular arrangements linked to the connectivity of the cells. To this end, we obtained single-cell cartography representations of the distribution of scutoids 11 in the cysts (Figures 2K–2M). In this case, the high heterogeneity in the number of scutoids per cyst and their distribution did not enable the identification of any clear pattern using the 3D reconstructions of the cysts (Figure 2K). We then plotted the total number of cells and analyzed the distribution of their position along the x and y axes (Figure 2L) and the z axis (Figure 2M) of the cysts. We did not find differences in the distributions in the first case. However, our analysis detected a significant increase of the proportion of scutoids from top to bottom of the cysts (Figure 2M) that was not observed in non-scutoid cells (chi-squared test) (Table S5 and STAR Methods). Our results suggest that cells pack following self-organization patterns in the MDCK cysts. Environmental perturbations can alter cell morphology patterns After finding cell morphology patterns in MDCK cysts, we wondered whether they remained consistent when the cell-culture environment was modified. Specifically, we investigated the impact of hypoxia ([O 2 ] = 1%) on the architecture and organization of 3D cysts (Figure 3A and STAR Methods), since it has been shown that hypoxia can affect cystogenesis by altering polarity in MDCK cell cultures 65 and even inducing epithelial branching. 66 In our cultures, we found a significant reduction in cell number and the cyst and lumen sizes throughout the whole experimental observation period (Figure 3B and Table S6). We then employed the single-cell cartography approach to analyze 7,729 cells from 206 segmented and curated hypoxic cysts and compared the patterns in both normoxia and hypoxia. We observed similar patterns in the case of ‘‘cell height’’ (compare Figures 2D, 2E, and S6;Table S6;STAR Methods), as well as ‘‘cell basal,’’ ‘‘cell volume,’’ ‘‘cell solidity,’’ and ‘‘cell aspect ratio’’ (Table S6 and STAR Methods). However, we found a different pattern on the distribution of the ‘‘cell apical area’’ values (compare Figures 2F, 2G, 3C, and 3D). Under hypoxic conditions, cells with the bigger apical area were enriched at the top and bottom regions (as in normoxia), but also in the middle region of the cysts (Figure 3D, Table S6, and STAR Methods). Furthermore, in contrast to normoxia, where significant differences in the distribution along the z axis of the cyst were found (magenta) and anti-b-catenin antibody (green) (STAR Methods). Scale bars, 100 mm. (B) Quantification of cyst size in normoxic (NX, black) and hypoxic (1% O 2 ) (HX, gray) cysts at 4, 7, and 10 days. Mean and SD are shown for number of cells per cyst, cyst volume, and lumen volume. Data were obtained from over 60 segmented cysts at each time point, from at least three independent experiments. **p %0.05, ****p %0.001 (Mann-Whitney U test). (C and D) Same representations as Figures 2D–2G to map the feature ‘‘cell apical area’’ in hypoxic (1% O 2 ) MDCK cysts. 7,729 cells from 206 segmented cysts were represented. (E and F) CartoCell segmentation of mouse embryoids (E) and Drosophila egg chambers (F). (Top left) Middle section of confocal microscopy z-stack images, where in (E), cell contours were stained with Alexa Fluor 647 phalloidin (magenta) and anti-b-catenin antibody (green) and in (F) with Resille-GFP (STAR Methods). (Top right) 2D segmentation of the previous section. (Bottom left) Half projection of z-stack images with the same stains for cell contours as on top. (Bottom right) 3D computer rendering. Scale bars, 10 mm. See also Figure S6;Tables S6 and S7. 8Cell Reports Methods 3, 100597, October 23, 2023 Article ll OPEN ACCESS STAR+METHODS KEY RESOURCES TABLE RESOURCE AVAILABILITY Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Luis M. Escudero ([email protected]). Materials availability No new materials were generated in this study. Data and code availability dAll data used in our analysis has been deposited at Mendeley Data and are publicly available as of the date of publication. DOIs are listed in the key resources table. REAGENT or RESOURCE SOURCE IDENTIFIER Antibodies Donkey anti-rabbit Alexa Fluor 488 Thermo Fisher Scientific Thermo Fisher Scientific Cat# A-21206, RRID: AB_2535792 Goat anti-rabbit Alexa Fluor 488 Thermo Fisher Scientific Thermo Fisher Scientific Cat# A-11034, RRID:AB_2576217 Rabbit monoclonal anti-b-catenin Sigma-Aldrich Sigma-Aldrich Cat# C2206, RRID:AB_476831 Rabbit monoclonal anti-b-catenin Santacruz Biotechnology Santa Cruz Biotechnology Cat# sc-7199, RRID:AB_634603 Chemicals, peptides, and recombinant proteins Phalloidin Alexa Fluor 647 Thermo Fisher Scientific Thermo Fisher Scientific Cat# A-22287 Phalloidin Alexa Fluor 555 Thermo Fisher Scientific Thermo Fisher Scientific Cat# A34055 DAPI Sigma-Aldrich Sigma-Aldrich Cat# 268298 Deposited data Raw and segmented images Mendeley Data repository https://doi.org/10.17632/7gbkxgngpm https://data. mendeley.com/datasets/7gbkxgngpm Experimental models: Cell lines Dog: Madin-Darby canine kidney cell line type II ECACC (ECACC Cat# 00062107, RRID: CVCL_0424) Experimental models: Organisms/strains Mouse: C57BL6/J The Jackson Laboratory RRID: IMSR_JAX:000664 D.melanogaster: Resille-GFP: P{PTT-un}CG8668117-2 Bloomington Drosophila Stock Center RRID:BDSC_99503; FlyBase: FBti0141278. Software and algorithms Fiji Schindelin et al., 2012 62 https://imagej.net/software/fiji/ MATLAB 2021a MathWorks https://mathworks.com PlantSeg Wolny et al., 2020 23 https://github.com/hci-unihd/plant-seg StarDist Weigert et al., 2020 25 https://github.com/stardist/stardist Cellpose Stringer et al., 2021 26 https://github.com/MouseLand/cellpose Original code (CartoCell pipeline) This paper, GitHub, FrancoBarranco et al. (2023) 82 https://github.com/danifranco/BiaPy https://biapy.readthedocs.io/en/latest/ tutorials/cartocell.html Original code (Single Cell Cartography Tools) This paper, Mendeley Data, GitHub https://doi.org/10.17632/7gbkxgngpm https://github.com/ ComplexOrganizationOfLivingMatter/NaturalVariation/ releases/tag/CartoCell_SingleCellCartographyTools_1.0.0 GraphPad Prism version 8.4.2 Dotmatics www.graphpad.com Other 3D ResU-Net Franco-Barranco et al., 2021 63 https://github.com/danifranco/EM_Image_Segmentation Cell Reports Methods 3, 100597, October 23, 2023 e1 Article ll OPEN ACCESS dAll original code used in our analysis has been deposited at Mendeley Data and is publicly available as of the date of publication. DOIs are listed in the key resources table. dAny additional information required to reanalyze the data reported in this paper is available from the lead contact upon request. EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS MDCK cyst cell culture Type II MDCK (Madin-Darby canine kidney) cells were maintained in minimum essential medium (MEM) containing GlutaMAX (Gibco) and supplemented with 10% fetal bovine serum (FBS), 100 U/ml penicillin and 100 mg/mL streptomycin, in a 5% CO 2 humidified incubator at 37C. For cyst formation, MDCK cells (2500 cells/well) were suspended in a complete medium containing 2% Matrigel (Corning, Life Sciences). Cell suspension was plated in a 4-well culture slide (Corning, Life Science) on a thin layer coating of 100% Matrigel. The plates were kept at 37C in a humidified atmosphere of 5% CO 2 for 4, 7 or 10 days and the medium was changed every 2 days. The cysts under hypoxia conditions ([O 2 ] = 1%) were maintained in incubators that allowed a precise and stable control of temperature, as well as O 2 and CO 2 concentrations. The exposure times to hypoxia were 4, 7, or 10 days, and the medium was changed every 2 days. Drosophila egg chambers For Drosophila ovary in vivo model, we used Resille-GFP (II) transgenic flies expressing a membrane marker tagged with the green fluorescence protein (GFP) ubiquitously to visualize the follicle cells of the egg chambers that compose the ovaries. 83 Drosophila egg chambers progress through 14 morphologically distinct stages of development. For in vivo acquisition, 1-2-day-old females were kept in a new food vial with yeast for optimal ovary development for 48 h. The ovaries were dissected in M3 insect culture medium (S83981L, Sigma) supplemented with 10% FBS (10270106; Gibco) and 0.20 mg/mL insulin (I550050G; Sigma). Ovarioles were isolated by removing the muscle sheath covering the ovary to avoid contraction movements during image acquisition. The ovarioles from different ovaries were mounted in a glass-based dish (Thermo Fisher Scientific), covered with supplemented medium to prevent the evaporation during in vivo acquisition. We obtained and processed confocal stacks from stage 2 to stage 7 egg chambers (20 egg chambers analyzed, Table S7). It should be noted that we distinguished between the early stages 2–3 (with an average cell count of 107.6 ±31.5), middle stages 4–5 (315.5 ±103.6 cells/egg chamber) and late stages 6–7 (853.8 ±126.2 cells/egg chamber) (Table S7). Mouse embryoids mES cell derivation For WT mouse embryonic stem (mES) cell derivation, eight-cell stage mouse embryos were recovered from the oviducts of pregnant females and cultured in KSOM (MR-020P-5F, Millipore) containing the inhibitors 2i/LIF to preserve naive pluripotency: 1mM MEK inhibitor PD0325901 (72182, STEMCELL Technologies), 3mM GSK3 inhibitor CHIR99021 (72052, STEMCELL Technologies) and 10 ngml1 LIF (Qk019, Qkine). After 24h, the medium was changed to N2B27 medium with 2i/LIF for 48h. N2B27 medium was comprised of a 1:1 mix of DMEM F12 (21331-020, Thermo Fisher Scientific) and neurobasal A (10888-022, Thermo Fisher Scientific) supplemented with 1% v/v B27 (10889-038, Thermo Fisher Scientific), 0.5% v/v N2 (homemade), 100mMb-mercaptoethanol (31350010, Thermo Fisher Scientific), penicillin–streptomycin (15140122, Thermo Fisher Scientific) and GlutaMAX (35050061, Thermo Fisher Scientific). Hatched blastocysts were next plated on mitomycin C-treated MEFs in Fc medium containing DMEM (41966, Thermo Fisher Scientific), 15% FBS (Stem Cell Institute), penicillin–streptomycin (15140122, Thermo Fisher Scientific), GlutaMAX (35050061, Thermo Fisher Scientific), MEM non-essential amino acids (11140035, Thermo Fisher Scientific), sodium pyruvate (11360070, Thermo Fisher Scientific) and 100mMb-mercaptoethanol (31350-010, Thermo Fisher Scientific. 2i/LIF was also added to Fc medium. Two days later, blastocysts outgrowths were trypsinized and plated to obtain mES cell colonies. mES cell culture mES cells were routinely cultured in gelatin-coated plates in Fc medium supplemented with 2i/LIF at 37C, 5% CO 2 , 21% O 2 . Cells were routinely tested for mycoplasma contamination. For embryoid formation, mES cells (20,000 cells/well) were suspended in a complete medium containing 5% Matrigel (Corning, Life Sciences). Cell suspension was plated in an 8-well culture slide (Corning, Life Science) on a thin layer coating of 100% Matrigel. The plates were kept at 37C in a humidified atmosphere of 5% CO 2 for 72h and the medium was changed every 2 days. Immunostaining and confocal imaging Cysts grown on 4-well chamber slides (Corning, Life Science) were fixed with 4% paraformaldehyde in PBS and permeabilized with 0.5% Triton X-100 in Dubelcco’s Phosphate Buffered Saline (DPBS, Sigma-Aldrich) for 15 min at room temperature (RT). After blocking with a solution of 0.02% Saponin (Sigma) and 3% BSA (Applichem) in DPBS for 2h at RT, cysts were incubated overnight at 4C with anti-b-catenin antibody (1:1000 in DPBS-0.02% Saponin-3% BSA; rabbit, Sigma-Aldrich). The following day, the cysts were washed with DPBS-0.02% Saponin-3% BSA solution (3x, 5 min each) and incubated for 90 min at RT in this solution plus anti-rabbit conjugated to Alexa Fluor 488 (1:800, Thermo Fisher Scientific) and phalloidin-Alexa-Fluor 647 (0.08 mM, Thermo Fisher Scientific). e2 Cell Reports Methods 3, 100597, October 23, 2023 Article ll OPEN ACCESS After washing with DPBS, plastic chambers were removed from microscope slides and coverslips were mounted onto slides using Fluoromount-G mounting solution (SouthernBiotech). Cysts were imaged using a Nikon Eclipse Ti-E laser scanning confocal microscope. High-resolution z stack confocal images were captured using a dry 340 objective (0.95 NA) and variable zoom (3.5–5.5), with a step size of 0.5 mm per slice and a scan speed of 0.25 ms, from top to bottom. Then, they were exported as nd2 files with an XY resolution ranging between 0.06 and 0.08 mm per pixel and an image size of 1024 31024 pixels. Low resolution images were captured (from top to bottom) using 320 oil objective (0.75 NA), step size of 0.7 mm per slice, scan speed of 0.5 ms and exported as nd2 files with an XY resolution of 0.62 mm per pixel and an image size of 1024 31024 pixels. Drosophila egg chambers, images were acquired using a Leica Stellaris 8 FALCON Confocal microscope at 25C, with 325 water objective (0.95 NA). The whole ovarioles were captured with a resolution of (0.61 mm per pixel in XY and 1.02 mm between Z slices, using the optimal z-step size. Laser compensation was applied to maintain similar levels of GFP intensity along the z axis. Embryoids grown on 8-well chamber slides (Corning, Life Science) were fixed with 4% paraformaldehyde in PBS and permeabilized with PBS +0.2% Triton Tx-100 + 0.2% SDS for 10 min at RT. After blocking with a solution of PBS +3% BSA for 2h at RT, cysts were incubated overnight at 4C with anti-b-catenin antibody (1:1000 in 3% BSA; rabbit, Santacruz). The following day, the cysts were washed three times with PBS and incubated for 120 min at RT in PBS +3% BSA with Alexa Fluor 488 anti-Rabbit (1:1000, Thermo Fisher Scientific), Alexa Fluor 555 Phalloidin (1:1000, Thermo Fisher Scientific) and DAPI (1:2000, Sigma-Aldrich). Cysts were imaged using a Nikon Eclipse Ti-E laser scanning confocal microscope. Low resolution images were captured (from top to bottom) using 320 oil objective (0.75 NA), step size of 0.7 mm per slice, scan speed of 0.5 ms and exported as nd2 files with an XY resolution of 0.62 mm per pixel and an image size of 1024 31024 pixels. METHOD DETAILS Custom DNN architecture: 3D ResU-Net Building upon the state of the art, we have designed 3D ResU-Net, a stable 3D residual U-Net 63 to segment epithelial cysts at the cell level. The architecture is presented in Figure S1. More specifically, 3D ResU-Net is formed by full pre-activation residual blocks (two 333 convolutional layers with a shortcut as shown in Figure S1), with 52 filters in the first level and adding 16 more at each level, dropout of 0.1 at each block. Down-sampling (max-pooling) operators are performed only in 2D, since the input volumes are anisotropic. The total number of trainable parameters is 1.3M. The network received raw cyst images as input and outputs three different channels: i) cell masks, with the probability that a voxel belongs to an individual cell, ii) contour, containing the probabilities of cell outlines, and iii) cell region, representing the foreground probability of the complete cyst. Network optimization To find the best solutions with our custom 3D ResU-Net, we made an exhaustive search of hyperparameters (Table S2) and training configurations, exploring different loss functions, optimizers, learning rates, batch sizes, and data augmentation techniques. In particular, we minimized the binary cross-entropy (BCE) loss using the Adam optimizer, with a learning rate of 0.0001, a batch size value of 2 and using a patch size of 80 380380 voxels. We used a Tesla P40 GPU card to train the network until convergence, i.e., for 1300 epochs with a patience established at 50 epochs monitoring the validation loss and picking up the model that performs best in the validation set (2 samples of ‘‘training high-resolution dataset’’ were used for model M1 and model M2 validation). Moreover, we applied on-the-fly data augmentation with random rotations, vertical, horizontal and z axis flips and brightness distortions. Image preprocessing Before accessing the network, all raw images were preprocessed for contrast homogenization using Fiji 62 macros. In this preprocessing, a contrast adjustment was performed using the ’enhance contrast’ function with 0.3% of saturated pixels. Additionally, an 8-bit transformation is applied to them. In the case of high-resolution images, a down-sampling was applied using Fiji macros, transforming the variable images resolutions (0.06–0.08 mm per pixel in XY and 0.5 mm between Z slices) to have the pixel size of the low-resolution images (0.62 mm per pixel in XY and 0.7 mm between Z slices). This was done by calculating a correction factor that multiplies the size of the original image: ½newSizeX;newSizeY;newSizeZ/SizeX$lowResXY highResXY ;SizeY$lowResXY highResXY ;SizeZ$lowResZ highResZ Both preprocessing Fiji macros are available at the public repository of the laboratory (see data and code availability section). Note that for label down-sampling the resize command ‘‘Size’’ must have the interpolation method set as ‘‘None’’ and the option of ‘‘average when down-sampling’’ disabled. For raw images, however, the interpolation method should be set as ‘‘Bilinear’’ and the option of ‘‘average when down-sampling’’ ticked. In the case of the low-resolution images, before applying the aforementioned automatic contrast enhancement process, a manual preprocessing was performed in which each cyst was cropped. This procedure was carried out with Fiji by drawing the region of interest (ROI) using the "rectangle" tool and cropping using the "Crop" or "Duplicate" command. Cell Reports Methods 3, 100597, October 23, 2023 e3 Article ll OPEN ACCESS Inference and postprocessing Each patch of 80 380380 voxels was processed by the network to reconstruct the output to the original cyst image size applying a padding of 16 316 316. In this way, we avoided border effects on every patch, as reported in. 63 The DNN produced three different outputs representing the probabilities of the individual cell masks, contours and the whole cell region. We binarized the first two outputs based on a fixed set of threshold values (0.2 was experimentally found to work best) and created instance seeds to be fed to a marker-controlled watershed algorithm. As a final post-processing step, we ran the Voronoi algorithm along the three dimensions, using the resulting instances of the watershed algorithm as Voronoi seeds, and the cell region output (binarized using Otsu thresholding, 84 as the bounded region to be occupied by Voronoi cells (Figure S3). Thus, the unoccupied intercellular space of the cell region was filled by the nearest individual cells (Voronoi seeds). Training and test datasets acquisition ‘‘High-resolution label images’’ were obtained after segmentation of the ‘‘high-resolution raw images’’ (21 cysts) using LimeSeg, 61 a plugin of Fiji 62 for 3D segmentation, based on surface elements (‘‘Surfels’’). This software was sourced from a set of seeds, manually placed over the volumetric image to localize every single cell. These seeds grow until the cell outlines are identified by detection of intensity gradient changes. The output of LimeSeg was processed using an in-house MATLAB program (2021a MathWorks) to detect and curate imperfections during cysts segmentation (see proofreading of segmented cysts section). A down-sampled version of these 21 segmented cysts (see images preprocessing section) was used for training the model M1 (Phase 2, Figure 1): 19 cysts composing the training dataset, and the remaining 2, making up the validation dataset. After running by default our high-content pipeline, we used the trained model M2 (Phase 5, Figure 1) to infer and subsequently segment the ‘‘test raw images’’ (60 cysts acquired at low resolution). These segmented images were manually curated using our in-house MATLAB proofreading program, obtaining the ‘‘test label images’’. Proofreading of segmented cysts A custom program developed in MATLAB and available in the public repository of the laboratory (see data and code availability section) was designed for the proofreading of segmented cysts. The software includes a user-friendly graphical user interface (GUI) that allows to remove, modify, merge or create labels by drawing on the two-dimensional slices of the image stack, also allowing the interpolation between labels on different slices for faster curations. Both cell and lumen labels can be modified using the GUI, which also has specific tools for each of them to ensure a proper visualization. In view that our biological study was developed on single-layer and single-lumen cysts, the proofreading software relies on a segmentation error detection tool specific to our purpose. The GUI displays the cell IDs of cells that do not contact the apical and/or basal surface of the cyst. The software was designed to work quickly on batches of cysts. Once the cyst stops displaying errors in the GUI and is marked by the user as fixed, the next cyst will be displayed in the GUI to be corrected. The procedure we carried out for the cyst curation started with the creation of 3 folders: One of them containing the batch of labels predicted by our pipeline, another one, the batch of raw images and a third one reserved for the curated labels. The software merges the raw images and the labels, and displays the result in a GUI along with information on possible segmentation errors. An expert reviewed the displayed image by carefully comparing each label with the staining of the cell membrane of the raw image, and adjusting the labels until a perfect segmentation was achieved. QUANTIFICATION AND STATISTICAL ANALYSIS Comparison with the state of the art Three state-of-the-art methods were tested against our protocol, being these methods PlantSeg, 23 Cellpose (Stringer et al., 2021) and StarDist 3D. 25 For a more robust comparison, each one of the methods were trained 10 times using default configuration values, and the 21 low-resolution cysts used to train our model M1 (Phase 2, Figure 1) were inputted as training dataset. Each of the trained models was evaluated on the same test set, composed by 60 perfectly annotated cysts, to obtain measures of the error in the results yielding the table (Table S3). The training of Cellpose was conducted locally by using a GPU (Graphics Processing Unit) and not using a pretrained model, as per the instructions provided in the training documentation (https://cellpose.readthedocs.io/en/latest/train.html). Inference was performed by following the instructions given in the command line documentation (https://cellpose.readthedocs.io/en/latest/ command.html) and using the diameter suggested by the Cellpose GUI. StarDist 3D training was performed in Google Colab using the official ZeroCostDL4Mic 72 implementation (https://github.com/ HenriquesLab/ZeroCostDL4Mic/wiki) using default values except for the following parameters: patch size, which was changed to 48 and patch height, which was changed to 32 for convenience given the size of the images to be used. PlantSeg training was performed locally following the training documentation instructions (https://github.com/hci-unihd/plant-seg) using as default configuration the 3D U-Net example for confocal imaging (https://github.com/wolny/pytorch-3dunet/blob/master/ resources/3DUnet_confocal_boundary/train_config.yml) replacing patch size to [32, 64, 64] as [z, x, y] for convenience given the size of the images used and the minimum values allowed by PlantSeg. Further to the training of the network, the PlantSeg GUI has e4 Cell Reports Methods 3, 100597, October 23, 2023 Article ll OPEN ACCESS modifiable parameters for the postprocessing. This part was performed in a custom way trying to optimize the watershed (done with Simple ITK) output obtained from the probability maps predicted by the network. The parameters used were: Under-/OverSegmentation Factor = 0.75, Run Watershed 2D = False, CNN Prediction Threshold = 0.113, Watershed Seeds Sigma = 1.0, Watershed Boundary Sigma = 0.4, Superpixels Minimum Size = 1, Cell Minimum Size = 5. These three methods with identical training configurations as previously described, have also been evaluated after training them with the same 314 cysts used for training M2 model (Phase 4, Figure 1). To ensure a robust comparison, ten models were trained for each method. Cyst features and shape classification Using an in-house MATLAB code, we quantified a set of geometrical and topological parameters of the segmented epithelial cysts (Table S1) as is graphically described in Figure S4. We carried out a classification of cysts depending on the morphology and differences between axes lengths (Figure 2A). We considered that two axes lengths were different if they differed more than 10%. We classified all cysts into 5 groups: 1. Sphere, when the lengths of the three axes of symmetry were similar. 2. Oblate, when two axes lengths were similar and the different one was the shortest axis length. 3. Prolate, when two axes lengths were similar and the different one was the longest axis length. 4. Ellipsoid, when the three axes lengths were different. 5. Negative curvature, when the solidity (volume/ convex volume) of the cyst was inferior to 0.9. Error evaluation of biological features We extracted the features of both manually curated cysts (ground-truth) and the output of our high-content segmentation pipeline (without proofreading). Some of the segmented cysts without curation presented under-segmentation that promoted gaps in the segmented tissue. This defective segmentation was called "cyst opening". These gaps prevented the identification of the lumen of the ‘‘open cysts’’ automatically, and thus some biological features could not be extracted. The 13% (46 cysts) of the automatically segmented cysts presented this defect, and they were not used in the comparison of the biological features values (Figure S3). For the remaining 87% of cysts (307 cysts), features were automatically extracted and compared with the features extracted from manually curated cysts. For each cyst, measurements of every feature were compared by computing the relative error calculated as relative error = jpredicted groundtruthj groundtruth . In the particular case of percentage of scutoids, we could not calculate the relative error because in some cases this feature represented a 0%, resulting in indetermination. Therefore, for the calculation of its relative error, we defined the complementary of this feature (100% - percentage of scutoids) such that we did not find any cyst with the 100% of cells being scutoids. Finally, we calculated the mean and standard deviation of the errors for every feature (Table S4). Single-cell cartography representation We performed an analysis of the spatial distribution of features from more than 20,000 cells from 353 segmented cysts using our Single-cell Cartography tools available in the public repository of the laboratory (see data and code availability section). Different types of representations arose from the use of these tools. Computer rendering of 3D cysts We displayed a 3D visualization of the segmented cysts and, using a gradient of color over the cell surfaces (Figures 2D, 2F, 2K, 3C, and S6). We can plot the normalized value of individual cell features, using our custom MATLAB function Paint3D. A batch processing of cysts allowed the creation of large sheets with the previously described three-dimensional representations of cysts on which to perform a visual pattern analysis. Polar plots We used two types of two-dimensional polar plots. Polar scatterplots and polar histograms were used to represent the relative spatial position (Figure 2C) and frequencies of a normalized cell feature for all cells of all cysts simultaneously (Figures 2E, 2G, 3D, S5,andS6). For the creation of these two-dimensional polar plots (both polar scatterplot and polar histograms) we proceeded as follows: The polar coordinate center for each cyst was set at the centroid of the cyst. The radius was normalized from 0 to 1, being 1 the distance to the farthest cell centroid from the centroid of the cyst. The colatitude angle (representing the height on the vertical axis) was calculated with respect to the horizontal plane passing through the cyst at the centroid, thus having positive angles for cells above the cyst centroid and negative angles for cells below the cyst (Figure 2C). The azimuthal angle (which rotates around the vertical axis) was ignored since the scope of the study was to search for patterns along the vertical axis. Disregarding this angle led to a two-dimensional representation. This approach consisted of 5 polar scatterplots and 4 polar histogram plots. First, a general polar scatterplot was shown in which all cells were represented (Figures 2E, 2G, 3D, S5,andS6). The value of the features was represented by a color gradient as in the previous case. The rest of the plots were dedicated to different ranges of the normalized feature to be studied: 0–0.25, 0.25–0.50, 0.50–0.75, 0.75–1. Each of the ranges was analyzed with a polar scatterplot and a polar histogram plot showing, normalized, the distribution of cells along the colatitude. In this way, we were able to visualize all the values of a particular cell feature distributed along the cysts vertical axis. Cell Reports Methods 3, 100597, October 23, 2023 e5 Article ll OPEN ACCESS Normalized cell spatial data For each cell, the following data were represented: the Z-position of the cell centroid with respect to the centroid of the lowest cell in the cyst (Figures 2H, 2I, 2M, and S6); the distance from the cell centroid to the vertical (Z) axis passing through the centroid of the cyst (Figures 2L and S6) and the value of the cellular characteristic to be studied regarding its spatial position (Figures 2H, 2I, and S6). The cell features were normalized regarding the maximum and minimum value of the feature in the whole cyst. Evaluation metrics To evaluate our results, we used common metrics to measure instance segmentation performance in 2D and 3D images, which are calculated by matching the ground-truth and prediction segmentation masks with an Intersection Over Union (IoU) value over a certain threshold. In particular, we show values that require at least 30%, 50% and 75% IoU with the ground-truth for a detection to be a true positive (TP) (Table S3). More specifically, we used the following metrics: Precision, defined as precision =TP TP+FP where TP and FP are the number of true and false positives, respectively. Recall, defined as recall =TP TP+FN where FN is the number of false negatives. Accuracy, defined as accuracy =TP TP+FP+FN F1 or F-score, defined as F1=23precision 3recall precision+recall Panoptic quality, a unified metric to express both segmentation and recognition quality, defined as in Equation 1 of 85 panoptic = Sðp;gÞeTPIoUðp;gÞ jTPj+1 2jFPj+1 2jFNj where p and g are the predicted and ground-truth segments, respectively. Therefore, 1 jTPjP ðp;gÞeTP IoUðp;gÞis the average of matched segments, and 1 2jFPj+1 2jFNjin the denominator penalizes segments without matches. Statistical analysis At each time-point sampled (4 days, 7 days and 10 days cysts), at least seven independent cultures were carried out for normoxic cysts, and three independent cultures under for hypoxic cysts. For comparisons of the features values on certain cell populations between different categories of cyst shapes, as well as to compare some features between hypoxic and normoxic cysts, we used a univariate statistical protocol. First, samples were evaluated for normal distribution and similar variance by using the ShapiroWilk test and the two-sample F-test, respectively. If samples followed a normal distribution and similar variance, we employed the two unpaired Student’s ttest; whereas data had a normal distribution but not equal variance, we used the two-tailed Welch test. Finally, when data not adjusting to a normal distribution, we employed the non-parametric Mann-Whitney Utest. Data were represented in a bar graph as mean ±SD (standard deviation) and p %0.05 was considered statistically significant (Figures 2J, 3B, and S6). ‘‘*’’, ‘‘**’’, ‘‘***’’ and ‘‘****’’ indicating p %0.05, p %0.01, p %0.001 and p %0.0001 respectively (Tables S5 and S6). In a different statistical analysis, we tested cell spatial distribution similarity from bottom to top (in z axis) or from cyst centroid to outside (in XY axis) of in the proportion of scutoids and non-scutoidal cells (Figures 2L, 2M, and S6;Tables S5 and S6). Similarly, we also compared the similarity in the spatial distribution of cells, ranging from a colatitude angle of 90 to 90, between normoxic and hypoxic cysts regarding the proportion of cells within each range of values of the normalized feature (Figures 2,3,S5, and S6). Following the guidelines from, 12,86 we used the chi-square test for the trend across all samples to determine if there is a linear trend for the proportional data, considering statistically significant p %0.05 (Tables S5 and S6). ‘‘*’’, ‘‘**’’, ‘‘***’’ and ‘‘****’’ indicating p %0.05, p %0.01, p % 0.001 and p %0.0001, respectively. Statistical analyzes and graphs were performed using GraphPad Prism version 8.4.2. (GraphPad Software, La Jolla California, USA, www.graphpad.com). e6 Cell Reports Methods 3, 100597, October 23, 2023 Article ll OPEN ACCESS