3D cell neighbour dynamics in growing pseudostratified epithelia

During morphogenesis, epithelial sheets remodel into complex geometries. How cells dynamically organise their contact with neighbouring cells in these tightly packed tissues is poorly understood. We have used light-sheet microscopy of growing mouse embryonic lung explants, three-dimensional cell segmentation, and physical theory to unravel the principles behind 3D cell organisation in growing pseudostratified epithelia. We find that cells have highly irregular 3D shapes and exhibit numerous neighbour intercalations along the apical-basal axis as well as over time. Despite the fluidic nature, the cell packing configurations follow fundamental relationships previously described for apical epithelial layers, that is, Euler's polyhedron formula, Lewis’ law, and Aboav-Weaire's law, at all times and across the entire tissue thickness. This arrangement minimises the lateral cell-cell surface energy for a given cross-sectional area variability, generated primarily by the distribution and movement of nuclei. We conclude that the complex 3D cell organisation in growing epithelia emerges from simple physical principles.


INTRODUCTION
Common to all animals and plants, epithelia are a fundamental tissue type whose expansion, budding, branching, and folding is key to the morphogenesis of organs and body cavities. Characterized by apicalbasal polarity (Figure 1a), epithelial cells adhere tightly to their apical neighbours in a virtually impermeable adhesion belt, form lateral cell-cell junction complexes along the apico-basal axis to provide mechanical stabilization, and bind tightly to the basal lamina and extracellular matrix (ECM) on the basal side (Drubin and Nelson, 1996;Rodriguez-Boulan and Macara, 2014;Shin and Margolis, 2006). How cell neighbour relationships are organised in these tightly adherent layers, and how these change during tissue and concomitant cell shape changes is poorly understood, despite their importance for cell-cell signalling and the fluidity of the tissue.
Cell neighbour relationships can be most easily studied on epithelial surfaces, and the polygonal arrangements of apical surfaces (Figure 1b) have been meticulously analysed (Classen et al., 2005;Escudero et al., 2011;Etournay et al., 2015;Farhadifar et al., 2007;Gibson et al., 2006;Gómez-Gálvez et al., 2018;Heller et al., 2016;Kokic et al., 2019;Sanchez-Gutierrez et al., 2016). Widely considered to be a reliable proxy for three-dimensional (3D) cell shape, 3D epithelial cell shapes are often depicted as prisms with polygonal faces that retain the same neighbour relationships along the entire apico-basal axis ( Figure 1c). Cells in curved epithelia are pictured as frustra, which have the same number of sides, but different apical and basal areas. If the curvature differs substantially along the principal axes, as is the case in epithelial tubes, neighbour relationships must change along the apical-basal axis. Prismatoids accommodate the neighbour change at the surface, while scutoids undergo the neighbour change somewhere along the apicalbasal axis (Gómez-Gálvez et al., 2018). However, even though the curvature is the same in both principal directions of spherically shaped epithelia, the neighbour relationships still differ between the apical and basal sides (Gómez-Gálvez et al., 2018), suggesting that effects other than curvature must determine the 3D neighbour arrangements of cells in epithelia.
Given the challenges in visualising 3D neighbour arrangements, most studies to date have focused on apical cell arrangements, and have revealed striking regularities. First, even though the frequencies of neighbour numbers differ widely between epithelial tissues (Figure 1e), cells have on average (close to) six neighbours (Classen et al., 2005;Escudero et al., 2011;Etournay et al., 2015;Farhadifar et al., 2007;Gibson et al., 2006;Heller et al., 2016;Kokic et al., 2019;Sanchez-Gutierrez et al., 2016) (Figure 1f). This can be explained with topological constraints in contiguous polygonal lattices, as expressed by Euler's Formula (Gibson et al., 2006;Rivier and Lissowski, 1982). Thus, if three cells meet at each vertex, the average number of neighbours in infinitely large contiguous polygonal lattices is exactly ̅ = 6. (1) While the average number of neighbours in the entire lattice is (close to) to six, the local averages deviate from six, and instead rather closely follow a phenomenological relationship, termed Aboav-Weaire's law (Aboav, 1970). According to Aboav-Weaire's law (Figure 1g), the average number of neighbours of all n cells that border a cell with n neighbours follows as Finally, the average apical area, ̅̅̅̅ , of cells with n neighbours is linearly related to the number of cell neighbours, n ( Figure 1h, black line), a relation termed Lewis' law (Lewis, 1928), Here, ̅ refers to the average apical cell area in the tissue.
We have recently shown that Aboav-Weaire's law and Lewis' law are a direct consequence of a minimisation of the lateral cell-cell contact surface energy Vetter et al., 2019). The lowest lateral cellcell contact surface energy is obtained in a regular polygonal lattice because regular polygons have the smallest perimeter per polygonal area. The distribution of apical cell sizes that emerges from cell growth and division is, however, such that epithelial tissues cannot organise into perfectly regular polygonal lattices. By adhering to Aboav-Weaire's law and Lewis' law, cells assume the most regular lattice. In particular, by following, Aboav-Weaire's law, the internal angles are closest to that of a regular polygon, while adding up to 360° at each tricellular junction (Figure 1i) (Vetter et al., 2019). And by following the relationship between polygon area and polygon type as stipulated by Lewis' law (Eq. 2), the difference in side lengths, ̅ ̅ ⁄ , is minimized between cells (Figure 1j) . The side lengths would be equal (Figure 1j, yellow line), if cells followed a quadratic relation of the form This quadratic relation (Figure 1h, yellow line), however, requires a larger area variability than is observed in most epithelia imaged to date. Accordingly, the predicted quadratic relation had not been previously reported, but could be confirmed experimentally by us by increasing the apical area variability .
Given the relationship between apical area and neighbour numbers as stipulated by Eqs. 1,3,4, the apical area variability emerges as the key determinant of apical epithelial organisation, and the theory correctly predicts how the fraction of hexagons in the tissue depends on the apical area variability, as can be quantified by the coefficient of variation (CV = std/mean) ( Figure 1k) . As such, growth and cell division determine the variability of the apical areas and thus determine apical organisation indirectly. Taken together, the apical organisation of epithelia can be understood based on the principles of lateral cell-cell contact surface energy minimisation.
In this work, we leverage these theoretical insights along with light-sheet fluorescence microscopy to study 3D epithelial organisation, both in cleared and growing pseudostratified epithelia. We find that cells have complex 3D shapes with numerous neighbour transitions along their apical-basal axis as well as over time.
We show that much as on the apical side, the variation of the cross-sectional areas along the apical-basal axis define the epithelial organisation at all times and across the entire tissue thickness. The observed neighbour arrangement minimizes the lateral cell-cell surface energy for a given cross-sectional area variability. The cross-sectional areas vary as a result of cell growth, division and interkinetic nuclear migration (IKNM). We conclude that the complex 3D cell organization in growing epithelia emerges from simple physical principles.

Apical and basal epithelial organisation
We started by exploring the apical and basal cellular organisation in epithelial tubes and buds ( Figure 2a). To this end, we imaged CUBIC-cleared mouse embryonic (E12.5) lung rudiments from a Shh GC/+ ; ROSA mT/mG background using light-sheet microscopy, and segmented the fluorescent membrane boundaries of over 400 cells per dataset in 2.5D (Figure 2b, Figure 2 -figure supplement 2). The apical and basal surfaces are both curved and thus differ in their total areas, i.e., the total segmented apical area is about 5-fold smaller than the basal area ( Figure 2b). We detected less than half as many cells on the apical side, and the mean cross-sectional cell area of apical cells is therefore on average only 2-fold smaller than that of basal cells, while the area variability, measured as area CV, is higher ( Figure 2c). Notably, the frequencies of the different neighbour numbers are not identical on the apical and basal side (Figure 2d), suggesting that the neighbour relationships change along the apical-basal axis, both in the tube and tip datasets. This observation is consistent with previous reports (Gómez-Gálvez et al., 2018). The change in neighbour relationships has previously been attributed to a curvature effect in tubes, but the neighbour changes in spherical geometries cannot be explained with such an effect (Gómez-Gálvez et al., 2018), suggesting that mainly other effects determine epithelial organisation.
So how can we explain the difference in apical and basal epithelial organisation in both datasets? We have previously shown that the apical organisation emerges from the minimisation of the overall lateral cell-cell contact surface energy Vetter et al., 2019). Aboav-Weaire's law (Eq. 2, Figure 1g) and Lewis' law (Eqs. 3,4, Figure 1h) emerge as global organization laws from this physical constraint, and ensure that the angles are closest to that of a regular polygon (Figure 1i, yellow lines), and that the side lengths are the most equal ( Figure 1j, yellow line). We now find that Aboav-Weaire's law ( Figure 2e) and Lewis' law ( Figure 2f) hold not only for the apical, but also for the basal datasets. Consistent with our theory, the apical layers, which have a larger area variability than the basal layers (Figure 2c), follow the quadratic law (yellow line) rather than the linear Lewis' law (black line).
We conclude that basal layers follow the same organisational principles as apical layers, such that their organisation can also be explained with a minimisation of the lateral cell-cell contact surface energy. Accordingly, the observed difference in overall neighbour relationships (Figure 2e,f) is a consequence of the difference in the cross-sectional area distributions ( Figure 2c). So, why do the normalised area distributions differ between the apical and basal sides in both the tube segment and the bud, and how do they change along the apical-basal axis?

3D organisation of epithelia
To explore the physical principles behind 3D epithelial cell organisation, we 3D segmented 140 cells from a tube segment and 59 cells from a bud segment in CUBIC-cleared, light-sheet imaged embryonic lung explants ( Figure 3a, Figure 2 -figure supplement 1, Video 1-4). By interpolating between equally spaced sequential contour surfaces (every 1.66 µm in the tube and every 1.72 µm in the bud dataset) along the apical-basal axis, accurate volumetric reconstructions of cell morphology were obtained that allowed for the extraction of morphometric quantifications along the apical-basal axis. In both datasets, the 3D organisation of epithelial cells is highly complex, and cell neighbour relationships change continuously along the apical-basal axis ( Figure 3b). As a result, cells are in direct physical contact not only with the cells that are neighbours on the apical side, but also with cells that appear two or even three cell diameters apart ( Figure 3c).
Remarkably, we record up to 14 cell neighbour changes per cell in the tube and up to 8 in the tip, between adjacent cross-sections along the apical-basal cell axis ( Figure 3d). We will refer to these neighbour changes as lateral T1 transitions, or T1L. The mean relative apical-basal position for the lateral T1 transitions is 0.4890.020 (95% CI), and there is no clear apical or basal tendency, though fewer transitions are observed close to the basal surface ( Figure 3e). The dispersion index, i.e. the ratio of the variance  2 and the mean number µ of transitions per cell, which equals unity for a Poisson distribution, is close to unity for both samples (Figure 3d). The chi-squared test also confirms that the number of apical-basal T1 transitions per cell is Poisson-distributed ( Figure 3d). A Poisson distribution models the probability of a number of independent random events occurring in a given interval at a constant average rate. The consistency with a Poisson distribution, therefore, suggests a stochastic basis to the 3D organization of epithelial cells.
The large number of observed T1L transitions and their distribution along the apical-basal axis challenges the recently popularized notion of curvature-driven scutoids as cell building blocks for epithelia (Gómez-Gálvez et al., 2018). To further examine the potential influence of tissue curvature on T1L transitions, we measured the apical-basal distance between two consecutive neighbour number changes for each cell in the tube dataset and recorded at which local tissue curvature they occur. For this analysis, we excluded cell portions from the apical end to the first transition and from the last transition to the basal end to reduce boundary effects, i.e., only interior segments between transitions were considered. The mean apical-basal distance between two transitions is 17.89  0.66 µm (95% CI). Local tissue curvature was approximated by fitting ellipses to the apical and basal surface boundaries of the tubular epithelium in 624 equidistant sections perpendicular to the main tube axis. The semi-axes of these ellipses were then averaged over all sections to obtain the semi-axes aapical, abasal, bapical, bbasal. Since our sample of 140 cells was segmented from a region close to the cusp of the nearly elliptical tube, a reasonably close estimate of the local tissue curvature where a T1L transition occurs is given by a linear interpolation between the curvature at the minor vertices of the apical and basal ellipses, according to the relative apical-basal position of the transition. The minor curvature of an ellipse with major and minor semi-axes a and b is given by b/a 2 . Therefore, we estimate the local radius of curvature R by where x[0,1] is the relative apical-basal location of the T1L transition. The examined tissue exhibits an average curvature fold change of R(1)/R(0)=2.21 from the basal to the apical side. Denoting by R1 and R2 the radii of curvature between two adjacent T1 transitions along the apical-basal axis of a cell (Figure 3f), we find that the distribution of curvature fold change R2/R1 shows no significant dependency on the number of neighbours n the cell has along that portion of the cell ( Figure 3g). The mean curvature fold change per apical-basal T1L transition per cell is <R2/R1> = 1.142  0.010 (95% CI). By extending the theory of scutoids (Gómez-Gálvez et al., 2018) to multiple T1L transitions per cell, we have derived a quantitative estimate of how tissue curvature would translate into the number of neighbour exchanges within that framework (Supplementary Material). If curvature changes were a main driver of T1L transitions, cells with smaller neighbour numbers n would be expected to change n over a much larger curvature fold change than cells with many neighbours (Figure 3g, blue line). However, we observe no systematic dependency of the curvature fold change on the number of neighbours the cell has along that portion of the cell in the developing mouse lung epithelium ( Figure 3g). From this, we conclude that tissue curvature affects cell neighbourhood rearrangements through the tissue thickness at most mildly.

Neighbour changes along the apical-basal axis are driven by changes in cross-sectional area variation
Other than curvature effects, what else could drive the observed changes in neighbour relationships along the apical-basal axis? We notice that much as the apical and basal layers, each layer along the apical-basal axis behaves according to the three relationships previously described for the apical side, i.e., Euler's formula (Eq. 1, Figure 4a), Aboav-Weaire's law (Eq. 2, Figure 4b), and Lewis' law (Eqs. 3, 4, Figure 4c). As predicted by the theory based on the minimisation of the lateral cell-cell energy , the layers with a large area variability ( Figure 4a) follow the quadratic law (yellow line) and those with a lower area variability the linear Lewis' law (black line). The fraction of hexagons also follows the predicted relationship with the cross-sectional area variability ( Figure 4d). We conclude that the entire 3D organisation of epithelia can be explained with a minimisation of the lateral cell-cell contact surface energy, as previously revealed for the apical layer.
If epithelial cell neighbour relationships are indeed driven by a minimisation of the total lateral cell-cell contact surface energy, then the T1L transitions along the apical-basal axes should be driven by changes in the cross-sectional area along the apical-basal axis. If we analyse four 3D segmented cells ( Figure 4e) in detail, we indeed see how an increase in the cross-sectional area results in an increase in the neighbour number, and vice versa ( Figure 4f) via lateral T1 transitions ( Figure 4g). As the cell neighbour arrangements represent global minima, the local analysis does, of course, not provide a perfect correlation. When we consider all 140 segmented cells in the tube segment with their 746 cell neighbour exchanges between adjacent cross-sections (Figure 3e), then we find that the frequency of T1L transitions along the apicalbasal axis is indeed higher, the larger the increase in cross-sectional area, and vice versa (Figure 4h,i).

Changes in cross-sectional area as a result of interkinetic nuclear migration (IKNM)
So, what determines the cross-sectional cell areas in each layer? In epithelia, mitosis is restricted to the apical surface. Depending on the average diameter of nuclei and the average apical cross-sectional area, there is insufficient space for all nuclei to be accommodated apically. Therefore, as a cell exits mitosis, the nucleus moves from the apical towards the basal side (G1 and S phase) and then back to the apical side (G2 phase) to undergo another round of mitosis, a process referred to as interkinetic nuclear migration (IKNM) (Meyer et al., 2011). Consequently, nuclei are distributed along the entire apical-basal axis, giving the tissue a pseudostratified configuration. We wondered to what extent the nuclear distribution, and its effect on the 3D cell shape, explains the observed area distributions and lateral T1 transitions.
To this end, we stained the nuclear envelope with fluorescently tagged antibodies against lamin B1 ( Figure 5d). Thus, a one-sided, two-sample Welch t-test revealed a significantly reduced ellipticity of nuclei located in the first 25% of the apical-basal axis compared to those in the middle 50% (p=0.0002). While the nuclear volumes ( Figure 5f) and cross-sectional areas ( Figure 5g) are slightly smaller than for the entire cell, the cross-sectional areas of the cell and the nucleus are strongly correlated (r = 0.94, Figure 5g). This is consistent with the nucleus determining the cross-sectional cell area, where present. Cell sections without nucleus typically have smaller cross-sectional areas, thereby leading to a higher frequency of small cross-sections in cells compared to nuclei. Accordingly, as previously seen for the cell cross-sectional areas, the observed changes in cell neighbour numbers correlates with the observed changes in nuclear cross-sectional areas ( Figure 5h) such that most T1L transitions occur where the nucleus starts and ends ( Figure 5i).
We conclude that the positions of nuclei can explain much of the observed variability in the cross-sectional cell areas. During the cell cycle, nuclei migrate, and the cell volumes first increase, and subsequently halve due to cell division. As all these processes affect the cross-sectional areas of the cells along the apical-basal axis, one would expect continuous spatial-temporal T1L transitions in growing pseudostratified epithelia.

3D cell organisation in growing epithelia
To follow 3D cellular dynamics during epithelial growth and deformation, we cultured embryonic lungs from a Shh GC/+ ; ROSA mT/mG background and imaged every 20 minutes for a total of 10 hours using lightsheet microscopy (Video 6). We used a subset of this dataset (11 time points, >3 hours) (Video 7) to 2.5D segment the apical and basal surfaces, and to explore 3D cell shape dynamics and neighbour relationships in a growing lung bud ( Figure 6a, Video 8).
As the explant was growing, we readjusted the 2.5D segmented region such that the segmented surface area and cell numbers remained roughly constant over time ( Figure 6 -figure supplement 1a). Nonetheless, the segmented bud increased in volume as the thickness of the layer increased with time ( Figure 6 -figure supplement 1b). Much as in the static dataset, the neighbour number distributions (Figure 6b), and variability of cross-sectional areas (Figure 6c) differ between the apical and basal cell layers in all time points. However, for all time points, both the apical and basal layers conformed to Euler's formula (Figure 6c), Aboav-Weaire's law (Figure 6d), and Lewis' law ( Figure 6e). Furthermore, the fraction of hexagons also followed from the variability of the cross-sectional areas, as predicted by the theory (Figure 6f).
We next sought to analyse the 3D dynamics of segmented epithelial cells. As the tracking of packed cells in growing pseudostratified epithelia is challenging, we focused on a small patch with 15 cells in total ( Figure  7a). Sequential contour surfaces were drawn to follow cell membrane outlines on several planes along the apical-basal axis and interpolated to reconstruct 3D morphology for each time point ( Much as in the static dataset ( Figure 3,4), we observe up to 14 neighbour number changes (T1L transitions) along the apical-basal axis (Figure 8a,b). The average number of T1L transitions is relatively constant over time ( Figure 8b). The mean relative apical-basal position for T1L transitions is again roughly in the middle, but in this small dataset, we now observe more T1L transitions in the center of the cell than at the apical or basal boundaries (Figure 8c). By following a single cell over time, we can appreciate the dynamic cell shape changes, and how a change in the cross-sectional area correlates with a change in neighbour number (Figure 8d). The neighbour relationships are, of course, not determined by the local cell cross-section, but by the overall cross-sectional area distribution in that layer. Accordingly, the correlation between the crosssectional area and the neighbour number is not perfect for a single cell. By considering a patch of cells, we can, however, see how those T1L transitions occur dynamically in developing tissues (Figure 8e).

DISCUSSION
Epithelial tissues remodel into complex geometries during morphogenesis. We used light-sheet microscopy and 3D cell segmentation to unravel the physical principles that define the 3D cell neighbour relationships in pseudostratified epithelial tissues. Our analysis reveals that pseudostratified epithelial layers adopt a far more complex packing solution than previously anticipated: the 3D epithelial cell shapes are highly irregular, and cell neighbour relationships change multiple times along the apical-basal axis, with some cells having up to 14 changes in their neighbor contacts along their apical-basal axis ( Figure  3a). Curvature effects can result in neighbour changes, but the data does not show the dependency on cell neighbour numbers that would be expected if curvature effects played a dominating role (Figure 3g). There is also no apical-basal bias (Figure 3e), and the prevalence of contact remodeling is randomly distributed (Figure 4i).
Even though the neighbour relationships are uncorrelated between the apical and basal sides and appear random at first sight, they follow the same fundamental relationships that have previously been described for apical epithelial layers, i.e. Euler's formula, Lewis' law, and Aboav-Weaire's law across the entire tissue and at all times (Figure 2, 4, 6, 7). This arrangement minimizes the lateral cell-cell surface energy in each plane along the apical-basal axis, given the variability in the cell cross-sectional areas ( Figure 1) Vetter et al., 2019). Where present, the stiff nucleus determines the cell cross-sectional area, as is apparent from the strong correlation between the cell cross-sectional and the nuclear cross-sectional areas ( Figure 5g). Accordingly, most changes in neighbour relationships occur at the apical and basal limits of the nucleus where cross-sectional areas change sharply ( Figure 5i). As the nucleus moves along the apical-basal axis during the cell cycle, a process referred to as interkinetic nuclear migration (IKNM) (Meyer et al., 2011), cell neighbour relationships change continuously (Figure 8). We conclude that neighbour relationships in epithelia are fluidic, and the complex, dynamic 3D organisation of cells in growing epithelia follows simple physical principles.
Defining the physical principles behind cell neighbour relationships is only the first step in unravelling the determinants of epithelial 3D cell shapes. The second key aspect is the cell volume distribution along the apical-basal axis, which gives rise to the cell cross-sectional area distribution, which then determines the cell neighbour relationships (Figure 5e). The overall cell volume is determined by cell growth and division, but its distribution along the apical-basal axis depends on the nuclear dynamics (Figure 5g), and the epithelial cell heights. We find that the nucleus occupies, on average, 55% of the cell volume in the embryonic lung epithelia. As the cell nuclei move along the apical-basal axis during the cell cycle (Meyer et al., 2011), the cytoplasm fills the remaining space between the apical and basal surfaces, likely in a way that minimises the total surface area of all cells. The determinants of the epithelial thickness, i.e. the distance between the apical and basal surfaces are still unknown, but signalling factors that control cellular tension are known to affect cell height (Widmann and Dahmann, 2009).
Cell-based modelling frameworks are heavily used to investigate epithelial processes and how they result in morphological changes such as tissue bending, folding, fusion, and anisotropic growth during morphogenesis (Fletcher et al., 2014;Tanaka, 2015). Our data confirms many underlying assumptions of cell-based modelling frameworks and provides quantitative data to calibrate parameters. Once calibrated to reproduce the here identified 3D cell shape distributions, such simulation frameworks will help to reveal the determinants of 3D cell shapes, and will be invaluable in providing insight into how local changes in cell growth, adhesion, tension, or in the basal lamina affect cell shapes locally and within the remaining epithelial layer.
In summary, this study offers a detailed view of 3D cell neighbor relationship dynamics and packing in growing epithelial tissues, and demonstrates that the 3D cell shapes are much more complex than previously anticipated, and that cell neighbor relationships are dynamic and change as result of cell growth and cell cycle-linked IKNM. The complex 3D cell neighbor relationships can nonetheless be understood based on simple physical principles. Although we recognize that tissue architecture is a multifactorial process, our work carries vast implications for the study of cell-cell signaling, epithelial cohesion, and energetic modeling of developing epithelial layers in both healthy and disease contexts.

Ethical Statement
Permission to use animals was obtained from the veterinary office of the Canton Basel-Stadt (license number 2777/26711). Experimental procedures were performed in accordance with the Guide for the Care and Use of Laboratory Animals and approved by the Ethics Committee for Animal Care of ETH Zurich. All animals were housed at the D-BSSE/UniBasel facility under standard water, chow, enrichment, and 12hrs light/dark cycles.

Animals
To investigate 3D cellular dynamics during mouse embryonic lung development, we used mouse lung rudiments from animals homozygous for the ROSA mT/mG and heterozygous for the Shh-cre allele (Shh cre/+ ; ROSA mT/mG ). The double-fluorescent Shh-controlled Cre reporter mouse expresses membrane-targeted tandem dimer Tomato (mT) before Cre-mediated excision and membrane-targeted green fluorescent protein (mG) after excision (Muzumdar et al., 2007). As a result, only epithelial cell membranes are labelled by GFP, while all adjacent mesenchymal tissue is labelled by tdTomato.

Optical clearing and Lightsheet imaging
Optical clearing of embryonic lung rudiments enabled the 3D segmentation of numerous epithelial cells from single image stacks. To this extent, the whole-mount clearing of dissected E12.5 lung explants was performed with the Clear Unobstructed Brain/Body Imaging Cocktails and Computational Analysis (CUBIC) protocol (Susaki et al., 2015) (Figure 2 -figure supplement 3). Reagents for delipidation and refractive index (RI) matching were made as follows: CUBIC-1 [25% (w/w) urea, 25% ethylenediamine, 15% (w/w) Triton X-100 in distilled water], and CUBIC-2 [25% (w/w) urea, 50% (w/w) sucrose, 10% (w/w) nitrilotriethanol in distilled water], respectively. Following fixation and immunostaining, samples were incubated in 1/2 CUBIC-1 (CUBIC-1:H2O=1:1) for four days, and in 1X CUBIC-1 until they became transparent. All explants were subsequently washed several times in PBS and treated with 1/2 CUBIC-2 (CUBIC-2:PBS=1:1) for around four days. Lastly, incubation in 1X CUBIC-2 was done until the desired transparency was achieved. All solutions were changed daily, and CUBIC-1 steps were performed on a shaker at 37 °C while CUBIC-2 steps at room temperature. Cleared samples were then embedded in 2% low melting point solid agarose cylinders and immersed in CUBIC-2 for two more days to increase the agarose refractive index. 3D image stacks were acquired on a Zeiss Lightsheet Z.1 microscope using a Zeiss 20x/1.0 clearing objective (Supplementary Figure 1).
Following a 1hr equilibration period, 1.5% LMP hollow agarose cylinders were prepared (Udan et al., 2014). Hollow cylinders, in contrast to solid ones, accommodate unencumbered 3D embryonic growth, provide boundaries to minimize tissue drift, enable imaging from multiple orientations, and allow for better perfusion of gasses and nutrients. All specimens were suspended within each hollow cylinder in undiluted Matrigel (VWR International GmbH; 734-1101), an ECM-based optically clear hydrogel that provided a near-native 3D environment and supported cell growth and survival. All cylinders were kept at 37ºC with 5% CO2 in culture media for 1hr before mounting.
For each overnight culture, the imaging chamber was prepared by sonication at 80ºC, followed by ethanol and sterile PBS washes. After assembly, the chamber was filled with culture medium and allowed to equilibrate at 37ºC with 5% CO2 for at least 2hrs before a cylinder containing an explant was mounted for imaging. Furthermore, to compensate for evaporation and to maintain a fresh culture media environment, two peristaltic pumps were installed to supply 0.4 mL and extract 0.2 mL of culture medium per hour. Each lung explant was then aligned with the focal plane within the center of a thin light-sheet to enable fine optical sectioning with optimal lateral resolution. For this study, all live imaging was done with a 20x/1.0 Plan-APO water immersion objective.

Image processing
To efficiently process the resulting volumetric CZI datasets (10s-100s of GBs), all image stacks were transferred to a storage server and subsequently processed in remote workstations (Intel Xeon CPU E5-2650 with 512 GB memory). Deconvolution via Huygens Professional v19.04 (Scientific Volume Imaging, The Netherlands, http://svi.nl) improved overall contrast and resolution while Fiji (ImageJ v1.52t) (Schindelin et al., 2012) aided in accentuating cell membranes, enhancing local contrast, removing background fluorescence, and TIFF conversion.

Cell morphometric quantifications
Cell morphology on the apical and basal membranes of embryonic lung epithelia was investigated using the open-source software platform MorphoGraphX (MGX) (Barbier de Reuille et al., 2015). By meshing the curved boundaries of input 3D image stacks and projecting nearby signal onto it, MGX builds a curved 2.5D image projection that is distortion-free, unlike planar 2D projections that ignore curvature. We then proceeded to use a suitable implementation of the Watershed transform to extract individual cell geometries, with minimal manual curation, and quantify properties such as surface area and the number of cell neighbours. All border cells were excluded. Apical and basal cell meshes were exported as text files and traversed using the R Programming Environment to extract the neighbour relationships between cells as needed to generate Aboav-Weaire plots.
To render time-lapse datasets and extract 3D volumetric surface reconstructions of entire epithelial cells, we employed Imaris v9.1.2 (Bitplane, South Windsor, CT, USA). By computationally interpolating between cell membrane contour surfaces from successive transverse frames into iso-surfaces, faithful cell and nuclear 3D volumes were obtained. Quantified volumetric features included cell and nuclear volume, total surface area, sphericity, and nuclear position along the apical-basal axis. Imaris was also used to generate high-resolution videos, which, despite being strongly downsampled to accommodate vast time-lapse datasets, presented little noticeable loss in image quality. Furthermore, to extract cell areas and the number of neighbours along the apical-basal axis, transverse image frames were imported into ImageJ and processed using the interactive plugin TissueAnalyzer (Aigouy et al., 2016). Like this, cell segmentation masks across layers could be generated, and cell geometry and neighbour topology quantified.

Table of contents -Figures
• Figure 1 -Principles of epithelial organization. The legend provides the measured average number of cell neighbours for each tissue and the references to the primary data (Classen et al., 2005;Escudero et al., 2011;Etournay et al., 2015;Farhadifar et al., 2007;Gibson et al., 2006;Heller et al., 2016;Sanchez-Gutierrez et al., 2016). Data points for n < 3 were removed as they must present segmentation artefacts. (f) The measured average number of cell neighbours is close to the topological requirement ( ̅ = 6) in all tissues; see panel e for the colour code. (g) Epithelial tissues follow the AW law (black line). The AW law formulates a relationship between the average number of neighbours, n, that a cell has and that its direct neighbours have, . The product • can be determined by summing over all n i . (h) The relative average apical cell area, ̅̅̅̅ ̅ ⁄ , increases with the number of neighbours, n, and mostly follows the linear Lewis' law (Eq. 2, black line), or the quadratic relationship (Eq. 3, yellow line) in case of higher apical area variability. (i) The average internal angle by polygon type is close to that of a regular polygon, = ( − 2)⁄ • 180° (yellow lines). To form a contiguous lattice, the angles at each tricellular junction must add to 360, and the resulting observed deviation in the angles follows the prediction (red line).