3D cell neighbour dynamics in growing pseudostratified epithelia
Abstract
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 lightsheet microscopy of growing mouse embryonic lung explants, threedimensional 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 apicalbasal 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 AboavWeaire's law, at all times and across the entire tissue thickness. This arrangement minimises the lateral cellcell surface energy for a given crosssectional 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. Characterised by apicalbasal polarity (Figure 1a), epithelial cells adhere tightly to their apical neighbours in a virtually impermeable adhesion belt, form lateral cellcell junction complexes along the apicobasal axis to provide mechanical stabilisation, and bind tightly to the basal lamina and extracellular matrix (ECM) on the basal side (Drubin and Nelson, 1996; RodriguezBoulan 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 cellcell 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ómezGálvez et al., 2018; Heller et al., 2016; Kokic et al., 2019; Ramanathan et al., 2019; SánchezGutiérrez et al., 2016). Widely considered to be a reliable proxy for threedimensional (3D) cell shape, 3D epithelial cell shapes are often depicted as prisms with polygonal faces that retain the same neighbour relationships along the entire apicobasal 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 apicalbasal axis. Prismatoids accommodate the neighbour change at the surface, while scutoids undergo the neighbour change somewhere along the apicalbasal axis (Figure 1c,d; GómezGá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ómezGá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; Ramanathan et al., 2019; SánchezGutiérrez et al., 2016; Figure 1f). This can be explained with topological constraints in contiguous polygonal lattices, as expressed by Euler's polyhedron 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
While the average number of neighbours in the entire lattice is (close to) six, the local averages deviate from six, and instead rather closely follow a phenomenological relationship, termed AboavWeaire's law (Aboav, 1970). According to AboavWeaire'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, $\overline{{A}_{n}}$, 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, $\overline{A}$ refers to the average apical cell area in the tissue.
We have recently shown that AboavWeaire's law and Lewis’ law are a direct consequence of a minimisation of the lateral cellcell contact surface energy (Kokic et al., 2019; 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 AboavWeaire's law and Lewis’ law, cells assume the most regular lattice. In particular, by following AboavWeaire'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 (Equation 2), the difference in side lengths, $\overline{S}}_{n}/\overline{S$, is minimised between cells (Figure 1j; Kokic et al., 2019). 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 (Kokic et al., 2019).
Given the relationship between apical area and neighbour numbers as stipulated by Equations 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; Kokic et al., 2019). As such, active processes such as growth, cell division, cell death and extrusion, cell intercalation and apical constriction 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 cellcell contact surface energy minimisation.
In this work, we leverage these theoretical insights along with lightsheet 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 apicalbasal axis as well as over time. We show that much as on the apical side, the variation of the crosssectional areas along the apicalbasal axis defines the epithelial organisation at all times and across the entire tissue thickness. The observed neighbour arrangement minimises the lateral cellcell surface energy for a given crosssectional area variability. The crosssectional areas vary as a result of active cell processes, most prominently including interkinetic nuclear migration (IKNM). We conclude that the complex 3D cell organisation in growing epithelia emerges from simple physical principles.
Results
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 CUBICcleared mouse embryonic (E12.5) lung rudiments from a Shh^{GC/+}; ROSA^{mT/mG} background using lightsheet 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, that is, the total segmented apical area is about 5fold smaller than the basal area (Figure 2b). We detected less than half as many cells on the apical side, and the mean crosssectional cell area of apical cells is therefore on average only twofold 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 apicalbasal axis, both in the tube and tip datasets. This observation is consistent with previous reports (GómezGálvez et al., 2018; Ramanathan et al., 2019; Rupprecht et al., 2017). 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 curvature alone, since prismatic cells fully accommodate the cell dilation from equal curvature changes in both directions (GómezGálvez et al., 2018; Rupprecht et al., 2017), 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 cellcell contact surface energy (Kokic et al., 2019; Vetter et al., 2019). AboavWeaire's law (Equation 2, Figure 1g) and the linear Lewis’ law or the quadratic relationship (Equations 3,4, Figure 1h) emerge as global organisation 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 these hold not only for the apical, but also for the basal datasets (Figure 2e,f). 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 cellcell contact surface energy. Accordingly, the observed difference in overall neighbour relationships (Figure 2e,f) is a consequence of the difference in the crosssectional 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 apicalbasal 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 CUBICcleared, lightsheet imaged embryonic lung explants (Figure 3a, Figure 2—figure supplement 1, Figure 2—video 1, Figure 3—videos 1–3). 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 apicalbasal axis, accurate volumetric reconstructions of cell morphology were obtained that allowed for the extraction of morphometric quantifications along the apicalbasal axis. In both datasets, the 3D organisation of epithelial cells is highly complex, and cell neighbour relationships change continuously along the apicalbasal 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 eight in the tip, between adjacent crosssections along the apicalbasal cell axis (Figure 3d). We will refer to these neighbour changes as lateral T1 transitions, or T1L. The mean relative apicalbasal 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, that is, 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 chisquared test also confirms that the number of apicalbasal T1 transitions per cell is Poissondistributed (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 organisation of epithelial cells.
The large number of observed T1L transitions and their distribution along the apicalbasal axis challenges the recently popularised notion of curvaturedriven scutoids as cell building blocks for epithelia (GómezGálvez et al., 2018). To further examine the potential influence of tissue curvature on T1L transitions, we measured the apicalbasal 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 ce ll portions from the apical end to the first transition and from the last transition to the basal end to reduce boundary effects, that is, only interior segments between transitions were considered. The mean apicalbasal 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. (Figure 3f). Epithelial tubes in the developing lung are often collapsed (Conrad et al., 2021), making their apical and basal surfaces nearly elliptic in shape. The semiaxes of the fitted ellipses were then averaged over all sections to obtain the semiaxes a_{apical}, a_{basal}, b_{apical}, b_{basal}. 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 apicalbasal position of the transition. The minor curvature of an ellipse with major and minor semiaxes 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 apicalbasal 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 R_{1} and R_{2} the radii of curvature between two adjacent T1 transitions along the apicalbasal axis of a cell, we find that the distribution of curvature fold change R_{2}/R_{1} 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 apicalbasal T1L transition per cell is <R_{2}/R_{1}> = 1.142 ± 0.010 (95% CI). By extending the theory of scutoids (GómezGá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 apicalbasal axis are driven by changes in crosssectional area variation
Other than curvature effects, what else could drive the observed changes in neighbour relationships along the apicalbasal axis? We notice that much as the apical and basal layers, each layer along the apicalbasal axis behaves according to the three relationships previously described for the apical side, that is, Euler's polyhedron formula (Equation 1, Figure 4a), AboavWeaire's law (Equation 2, Figure 4b), and Lewis’ law (Equations 3, 4, Figure 4c). As predicted by the theory based on the minimisation of the lateral cellcell energy (Kokic et al., 2019), 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 crosssectional area variability (Figure 4d). We conclude that the neighbour relationships of epithelial cells along the entire apicalbasal axis can be explained with a minimisation of the lateral cellcell contact surface energy, as previously revealed for the apical layer.

Figure 4—source data 1
 https://cdn.elifesciences.org/articles/68135/elife68135fig4data1v2.xls

Figure 4—source data 2
 https://cdn.elifesciences.org/articles/68135/elife68135fig4data2v2.xls

Figure 4—source data 3
 https://cdn.elifesciences.org/articles/68135/elife68135fig4data3v2.xls
If epithelial cell neighbour relationships are indeed driven by a minimisation of the total lateral cellcell contact surface energy then the T1L transitions along the apicalbasal axes should be driven by changes in the crosssectional area along the apicalbasal axis. If we analyse four 3D segmented cells (Figure 4e) in detail, we indeed see how an increase in the crosssectional 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 crosssections (Figure 3e), then we find that the frequency of T1L transitions along the apicalbasal axis is indeed higher, the larger the increase in crosssectional area, and vice versa (Figure 4h,i).
Changes in crosssectional area as a result of interkinetic nuclear migration (IKNM)
So, what determines the crosssectional cell areas in each layer? In pseudostratified epithelia, mitosis is restricted to the apical surface (Gundersen and Worman, 2013). Depending on the average diameter of nuclei and the average apical crosssectional 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 apicalbasal axis, giving the tissue a pseudostratified configuration (Norden, 2017). 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 5a, Figure 5—figure supplement 1), and 3D segmented all nuclei within epithelial cells in a tube segment (Figure 5b, Figure 3—figure supplement 1, Figure 3—video 2). The nuclei were distributed along the entire apicalbasal axis (Figure 5c), and consistent with the pseudostratified appearance of the epithelium, nuclei in neighbouring cells had different positions along the apicalbasal axis (Figure 5b). The nuclear shapes, volumes, and crosssectional areas (Figure 5d–f) all varied along the apicalbasal axis. As expected, nuclei are largest and most spherical at the apical side, where they undergo mitosis (Figure 5d). Thus, a onesided, twosample Welch ttest revealed a significantly reduced ellipticity of nuclei located in the first 25% of the apicalbasal axis compared to those in the middle 50% (p=0.0002). The nuclear volumes are on average about 50% smaller than the cell volumes and largely correlate (r = 0.79, Figure 5f), likely reflecting parallel expansion during the cell cycle. Where present, the nuclear crosssectional areas are only slightly smaller than those of the entire cell, and the crosssectional areas of the cell and the nucleus are strongly correlated (r = 0.94, Figure 5g). The strong correlation can be accounted for by the opposing actions of cells and nuclei in the columnar epithelium. The nuclear volumes are too large to allow for a spherical nucleus to fit into a cylindrical cell of the measured height (Figure 5h). Accordingly, to fit into the cell, the nucleus necessarily has to deform. Nuclei respond to external forces with anisotropic shape changes (Haase et al., 2016; Neelam et al., 2016), which is consistent with the elliptical nuclear shapes that we observe (Figure 5d). However, there is a limit to how much the stiff nucleus can deform (Lammerding, 2011; Shah et al., 2021), resulting in a local widening of the cell where the nucleus is present. Cell sections without nucleus typically have smaller crosssectional areas, thereby leading to a higher frequency of small crosssections in cells compared to nuclei. Accordingly, as previously seen for the cell crosssectional areas, the observed changes in cell neighbour numbers correlates with the observed changes in nuclear crosssectional areas (Figure 5i) such that most T1L transitions occur where the nucleus starts and ends (Figure 5j).

Figure 5—source data 1
 https://cdn.elifesciences.org/articles/68135/elife68135fig5data1v2.xls

Figure 5—source data 2
 https://cdn.elifesciences.org/articles/68135/elife68135fig5data2v2.xls

Figure 5—source data 3
 https://cdn.elifesciences.org/articles/68135/elife68135fig5data3v2.xls
We conclude that the positions of nuclei can explain much of the observed variability in the crosssectional 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 crosssectional areas of the cells along the apicalbasal axis, one would expect continuous spatialtemporal 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 min for a total of 10 hr using lightsheet microscopy (Figure 6—video 1). We used a subset of this dataset (11 time points, >3 hr) (Figure 6—video 2) 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, Figure 6—video 3).

Figure 6—source data 1
 https://cdn.elifesciences.org/articles/68135/elife68135fig6data1v2.xls

Figure 6—source data 2
 https://cdn.elifesciences.org/articles/68135/elife68135fig6data2v2.xls

Figure 6—source data 3
 https://cdn.elifesciences.org/articles/68135/elife68135fig6data3v2.xls
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 crosssectional 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 polyhedron formula (Figure 6c), AboavWeaire's law (Figure 6d), and Lewis’ law (Figure 6e). Furthermore, the fraction of hexagons also followed from the variability of the crosssectional 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 apicalbasal axis and interpolated to reconstruct 3D morphology for each time point (Figure 7—figure supplement 1, Figure 7—video 1). All planar segmentations were then pooled into five groups along the apicalbasal axis to enable morphometric analysis in different tissue regions. Over the time course, the volume of individual cells varied between roughly 400 and 800 µm^{3} (Figure 7b), and the apicalbasal length varied between roughly 20 and 30 µm (Figure 7c). We note that all layers conform to AboavWeaire's law (Figure 7d) and Lewis’ law (Figure 7e) in all time points. Moreover, consistent with our theory, the fraction of hexagons follows from the variability of the crosssectional area, though more deviations are observed, given the small number of cells analysed (Figure 7f).

Figure 7—source data 1
 https://cdn.elifesciences.org/articles/68135/elife68135fig7data1v2.xls

Figure 7—source data 2
 https://cdn.elifesciences.org/articles/68135/elife68135fig7data2v2.xls

Figure 7—source data 3
 https://cdn.elifesciences.org/articles/68135/elife68135fig7data3v2.xls
Much as in the static dataset (Figures 3 and 4), we observe up to 14 neighbour number changes (T1L transitions) along the apicalbasal axis (Figure 8a,b). The average number of T1L transitions is relatively constant over time (Figure 8b). The mean relative apicalbasal 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 crosssectional area correlates with a change in neighbour number (Figure 8d). The neighbour relationships are, of course, not determined by the local cell crosssection, but by the overall crosssectional 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 lightsheet 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 apicalbasal axis, with some cells having up to 14 changes in their neighbour contacts along their apicalbasal 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 apicalbasal bias (Figure 3e), and the prevalence of contact remodelling 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, that is, Euler's polyhedron formula, Lewis' law, and AboavWeaire's law across the entire tissue and at all times (Figures 2, 4, 6 and 7). This arrangement minimises the lateral cellcell surface energy in each plane along the apicalbasal axis, given the variability in the cell crosssectional areas (Figure 1; Kokic et al., 2019; Vetter et al., 2019). Where present, the stiff nucleus determines the cell crosssectional area, as is apparent from the strong correlation between the cell crosssectional and the nuclear crosssectional areas (Figure 5g). Accordingly, most changes in neighbour relationships occur at the apical and basal limits of the nucleus where crosssectional areas change sharply (Figure 5i). As the nucleus moves along the apicalbasal 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 apicalbasal axis, which gives rise to the cell crosssectional 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 apicalbasal 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 apicalbasal 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, that is, the distance between the apical and basal surfaces are still unknown, but signalling downstream of Fibroblastic Growth Factor (FGF), Sonic Hedgehog (SHH), Bone Morphogenetic Protein (BMP)/ transforming growth factorbeta (TGFβ), and WNT has been observed to affect cell height, presumably via an effect on cell tension and/or cellcell adhesion (GritliLinde et al., 2002; Hirashima and Matsuda, 2021; Kadzik et al., 2014; Kondo and Hayashi, 2015; Widmann and Dahmann, 2009).
Cellbased 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 cellbased 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 neighbour 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 neighbour relationships are dynamic and change as result of cell growth and cell cyclelinked IKNM. The complex 3D cell neighbour 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 cellcell signalling, epithelial cohesion, and energetic modelling of developing epithelial layers in both healthy and disease contexts.
Materials and methods
Ethical statement
Request a detailed protocolPermission to use animals was obtained from the veterinary office of the Canton BaselStadt (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 DBSSE/UniBasel facility under standard water, chow, enrichment, and 12 hr light/dark cycles.
Animals
To investigate 3D cellular dynamics during mouse embryonic lung development, we used mouse lung rudiments from animals homozygous for the Rosa26^{mTmG} and heterozygous for the Shh^{GFPCre} allele (Shh^{GFPCre/+}; Rosa26^{mTmG}). The doublefluorescent Shhcontrolled Cre reporter mouse expresses membranetargeted tandem dimer Tomato (mT) before CREmediated excision and membranetargeted 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.
Immunofluorescence
Request a detailed protocolE12.5 mouse lungs were fixated for 1 hr in 4% paraformaldehyde in PBS, and subsequently incubated with Lamin B1 (Thermo; Material No. 702972; 1:200) at 4°C for 3 days. As a structural component of the nuclear lamina, LaminB1 immunostaining makes crowded nuclei clearly distinguishable and easily segmentable. After washing in DPBS, lungs were incubated with conjugated fluorescent secondary Alexa Fluor 555 donkey antimouse IgG (H+L) (Abcam; Material No. ab150106; 1:250) for 2 days at 4°C.
Optical clearing and lightsheet imaging
Request a detailed protocolOptical clearing of embryonic lung rudiments enabled the 3D segmentation of numerous epithelial cells from single image stacks. To this extent, the wholemount 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: CUBIC1 [25% (w/w) urea, 25% ethylenediamine, 15% (w/w) Triton X100 in distilled water], and CUBIC2 [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 CUBIC1 (CUBIC1:H2O=1:1) for four days, and in 1X CUBIC1 until they became transparent. All explants were subsequently washed several times in PBS and treated with 1/2 CUBIC2 (CUBIC2:PBS=1:1) for around four days. Lastly, incubation in 1X CUBIC2 was done until the desired transparency was achieved. All solutions were changed daily, and CUBIC1 steps were performed on a shaker at 37°C while CUBIC2 steps at room temperature. Cleared samples were then embedded in 2% low melting point solid agarose cylinders and immersed in CUBIC2 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 .
Timelapse lightsheet acquisitions
Request a detailed protocolLightsheet acquisitions of live epithelial cell morphology enabled the study of 3D organisation dynamics. Following dissection in DPBS at room temperature, E12.5 lung explants were cultured in sterile Dulbecco's modified Eagle's medium w/o phenol red (DMEM) (Life Technologies Europe BV; 11039021) containing 10% Fetal Bovine Serum (FBS) (SigmaAldrich Chemie GmbH; F9665500ML), 1% Glutamax (Life Technologies Europe BV; A1286001), and 1% penicillin/streptomycin (Life Technologies Europe BV; 10378–016). All specimens were equilibrated at 37°C with 5% CO2 in a humidified incubator for 1 hr.
Following a 1 hr 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 minimise 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 nearnative 3D environment and supported cell growth and survival. All cylinders were kept at 37°C with 5% CO_{2} in culture media for 1 hr 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% CO_{2} for at least 2 hr 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 lightsheet to enable fine optical sectioning with optimal lateral resolution. For this study, all live imaging was done with a 20x/1.0 PlanAPO water immersion objective.
Image processing
Request a detailed protocolTo efficiently process the resulting volumetric CZI datasets (10s100s of GBs), all image stacks were transferred to a storage server and subsequently processed in remote workstations (Intel Xeon CPU E52650 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
Request a detailed protocolCell morphology on the apical and basal membranes of embryonic lung epithelia was investigated using the opensource 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 distortionfree, 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 AboavWeaire plots (Gómez et al., 2021).
To render timelapse 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 isosurfaces, 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 apicalbasal axis. Imaris was also used to generate highresolution videos, which, despite being strongly downsampled to accommodate vast timelapse datasets, presented little noticeable loss in image quality. Furthermore, to extract cell areas and the number of neighbours along the apicalbasal 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 (Gómez et al., 2021).
Appendix 1
A simple theory for lateral T1 transitions in curved tissues
Drawing inspiration from the notion of ‘scutoids’ (GómezGálvez et al., 2018), we develop here a simple geometrical theory to estimate the effect of tissue curvature, if any, on the occurrence and number of lateral T1 transitions of cells along the apicalbasal axis in a pseudostratified epithelium. The purpose of this effort is to verify whether an idea brought forward in the aforementioned article, which is that tissue curvature could be responsible for the occurrence of up to one lateral T1 transition per cell between its apical and basal sides, is consistent with our observations in the developing, pseudostratified mouse lung epithelium.
We start by considering, like the authors of GómezGálvez et al., 2018, the vicinity of an edge shared by two adjacent cells in a twodimensional cross section of the tissue perpendicular to the apicalbasal axis. This vicinity is defined as the quadrilateral (Appendix 1—figure 1A, orange) spanned by the four cell junctions (green) that are directly linked to the edge of interest (blue), covering a fraction of four cells in such a twodimensional projection. As the projection plane is moved along the apicalbasal axis, the edge length can shrink to zero and reappear with different orientation, leaving the cell neighbourhood of all four cells modified by one. It is this process which we refer to as lateral neighbourhood transition of type T1. Each such edge vicinity is characterised by an aspect ratio $\u03f5=w/h$, where $w=({w}_{1}+{w}_{2})/2$ is the average width and $h=({h}_{1}+{h}_{2})/2$ the average height. We now proceed to quantitatively estimate with a simple geometrical model how a change in tissue curvature translates into a change of aspect ratio of such motifs, and how that in turn would be expected to change the cell neighbourhood if all T1 transitions were effectively governed by the geometric effect of changing tissue curvature.
We consider two different tissue topologies: (i) a cylindrical tube and (ii) a spherical vesicle or hemispherical tube cap. The crucial difference between these two cases is that the former only has a mean curvature but no Gaussian curvature, whereas the latter has both. As the radius $r$ of these surfaces changes, the change of aspect ratio follows by the quotient rule of calculus as
We now consider the general case where the edge motif is arbitrarily oriented on the surface, and will later average over all possible orientations. In general, the edge of interest spans an angle $\theta $ with the first principal curvilinear axis of the surface, which we take as the main cylinder axis for the tubular epithelium (Appendix 1—figure 1), and as any great circle on the spherical epithelium. The change of width and height then read
where w_{x}, w_{y} and h_{x}, h_{y} are the components of the width and height vectors, respectively (Appendix 1—figure 1B, gray):
Their response to a change in tissue radius is given by
where ${\kappa}_{1}$ and ${\kappa}_{2}$ are the principal curvatures of the surface in the respective directions. Substituting Equation 8 into Equations 7, 6 and 5 yields
A cylinder with radius $r$ has the principal curvatures ${\kappa}_{1}=0$ and ${\kappa}_{2}=1/r$, so
For a sphere with radius $r$, on the other hand, the two principal curvatures are both equal, ${\kappa}_{1}={\kappa}_{2}=1/r$, and thus the aspect ratio does not change with varying radius,
independent of the orientation $\theta $.
From here onward, we consider the region along the apicalbasal axis between two T1 transitions on a given cell, irrespective of whether these T1 transitions increase or decrease the neighbour number. In our experimental quantifications of lateral T1 transition in the mouse lung epithelium, they were counted in the same way. The sign of the T1 transition is encoded in the sign of $\partial \u03f5/\partial r$, hence we consider only the absolute value $\mathrm{\partial}\u03f5/\mathrm{\partial}r$ hereafter, to quantify the change of edge aspect ratio between T1 transitions along the tissue radius.
Carrying on with the tubular geometry, we now average over all orientations $\theta $ to calculate the mean change of aspect ratio for the entire cell ensemble in the tissue. The set of possible orientations is $\theta \in [0,\pi /2)$, because all angles $\theta \in [\pi /2,\pi )$ just switch the roles of $w$ and $h$. Therefore,
This is an ordinary differential equation which can easily be solved by separation of variables. The result is
and therefore
where ${\u03f5}_{1}$ is the aspect ratio at a given cylinder radius R_{1}. For the sphere, on the other hand, we trivially have a constant aspect ratio $\u03f5(r)={\u03f5}_{1}$.
To link this change of aspect ratio with lateral T1 transitions, a relationship between $\u03f5$ and the number of neighbours a cell locally has, $n$, needs to be established. A simple rough estimate can be found by recalling that cells typically tend toward a regular polygonal shape in a crosssectional projection (Kokic et al., 2019; Vetter et al., 2019). In a honeycomb lattice, the aspect ratio is
For general regular polygons with $n\ge 4$, the relationship reads
as can be recognised from Appendix 1—figure 1B. Equation 10 holds locally in the region between two consecutive T1 transitions, where $n$ is constant. From this, we can estimate the change of aspect ratio between two lateral T1 transition as
On average, for hexagonal cells ($n=6$), this evaluates to
So how many T1 transitions can be expected along the apicalbasal axis? We can apply the chain rule of calculus to write
Substituting Equations 9 and 11 into Equation 12, we find
for a cylindrical epithelium, and $dn/dr=0$ for a spherical one. A spherical topology is thus entirely compatible with cells being frusta that do not change their neighbour number from their apical to the basal side, whereas a cylindrical topology is not. For tubes, a change of radius ${R}_{1}\to {R}_{2}$ lets the neighbour number change by
per cell edge. How far must R_{1} and R_{2} be apart to expect one T1 transition per cell? Over that distance, each of the cell’s $n$ edges will have undergone a change of $1/n$ on average. The answer can therefore be found by setting $\mathrm{\Delta}n=1/n$ with ${\u03f5}_{1}=\u03f5(n)$ from Equation 10, and solving for the radius fold change:
This is the expected fold change of mean curvature inside a tubular epithelium between two consecutive lateral T1 transitions per cell, if the T1 transitions were mainly driven by a change of curvature. Note that the two cases with different signs in the exponent reflect the ${R}_{1}\leftrightarrow {R}_{2}$ symmetry. As shown in Appendix 1—figure 2, Equation 13 predicts a strong dependency of the required change of curvature between adjacent lateral T1 transitions on the number of neighbours the cell has, $n$. Cells with a large neighbour number are expected to transition at much smaller changes of curvature on average than cells with few neighbours.
Equation 13 relates the fold change of tissue curvature on a tubular sheet between consecutive lateral T1 transitions along the apicalbasal axis of a single cell, to the cell’s neighbour number $n$ in any cross section perpendicular to the radius between the two transitions. $n$ here is to be understood locally within the tissue, not as the cell connectivity on either the apical or basal surfaces of the cell. Equation 13 does not directly lend itself to estimating the number of lateral T1 transitions to be expected for cells in a curved tissue, although our theory might be adaptable to this case with some modifications. The sign of T1 transitions is ignored here, leaving the possibility of frequent backandforth transitions with zero net effect on $n$—a phenomenon that we indeed do observe in the developing lung epithelium. Equation 13 quantifies only the net contribution of curvature on the number of T1 transitions a cell undergoes. This net effect strongly depends on the the local neighbour number $n$, as shown in Appendix 1—figure 2.
Data availability
The source code and plotted data files are available as a git repository at https://git.bsse.ethz.ch/iber/Publications/2021_gomez_3d_cell_neighbour_dynamics.git. The raw data is publicly available as openBIS repository at https://openbisdatarepo.ethz.ch/openbis/webapp/elnlims/?user=observer&pass=openbis under the name 3D Epithelium.

openBISID 3D Epithelium. 3D Epithelium.
References

The arrangement of grains in a polycrystalMetallography 3:383–390.https://doi.org/10.1016/00260800(70)900388

Segmentation and quantitative analysis of epithelial tissuesMethods in Molecular Biology 1478:227–239.https://doi.org/10.1007/9781493963713_13

Epithelial organisation revealed by a network of cellular contactsNature Communications 2:526.https://doi.org/10.1038/ncomms1536

Vertex models of epithelial morphogenesisBiophysical Journal 106:2291–2304.https://doi.org/10.1016/j.bpj.2013.11.4498

Mechanisms of cell height changes that mediate epithelial invaginationDevelopment, Growth & Differentiation 57:313–323.https://doi.org/10.1111/dgd.12224

Vertical uniformity of cells and nuclei in epithelial monolayersScientific Reports 6:19689.https://doi.org/10.1038/srep19689

Pseudostratified epithelia  cell biology, diversity and roles in organ formation at a glanceJournal of Cell Science 130:1859–1863.https://doi.org/10.1242/jcs.192997

On the correlation between sizes and shapes of cells in epithelial mosaicsJournal of Physics A: Mathematical and General 15:L143–L148.https://doi.org/10.1088/03054470/15/3/012

Organization and execution of the epithelial polarity programmeNature Reviews Molecular Cell Biology 15:225–242.https://doi.org/10.1038/nrm3775

Geometric constraints alter cell arrangements within curved epithelial tissuesMolecular Biology of the Cell 28:3582–3594.https://doi.org/10.1091/mbc.e17010060

Fiji: an opensource platform for biologicalimage analysisNature Methods 9:676–682.https://doi.org/10.1038/nmeth.2019

Advanced CUBIC protocols for wholebrain and wholebody clearing and imagingNature Protocols 10:1709–1727.https://doi.org/10.1038/nprot.2015.085
Article and author information
Author details
Funding
Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (CRSII5_170930)
 Dagmar Iber
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work has been supported through an SNF Sinergia grant to DI. We acknowledge the many discussions and varied input of Marco Kokic, Anđela Markovic, Odyssé Michos, and Erand Smakaj, and we thank Richard Smith for help with MorphographX.
Ethics
Animal experimentation: Permission to use animals was obtained from the veterinary office of the Canton BaselStadt (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.
Version history
 Received: March 5, 2021
 Preprint posted: March 7, 2021 (view preprint)
 Accepted: September 27, 2021
 Accepted Manuscript published: October 5, 2021 (version 1)
 Version of Record published: November 5, 2021 (version 2)
Copyright
© 2021, Gómez et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,427
 views

 367
 downloads

 39
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Cell Biology
The spindle assembly checkpoint (SAC) temporally regulates mitosis by preventing progression from metaphase to anaphase until all chromosomes are correctly attached to the mitotic spindle. Centrosomes refine the spatial organization of the mitotic spindle at the spindle poles. However, centrosome loss leads to elongated mitosis, suggesting that centrosomes also inform the temporal organization of mitosis in mammalian cells. Here, we find that the mitotic delay in acentrosomal cells is enforced by the SAC in a MPS1dependent manner, and that a SACdependent mitotic delay is required for bipolar cell division to occur in acentrosomal cells. Although acentrosomal cells become polyploid, polyploidy is not sufficient to cause dependency on a SACmediated delay to complete cell division. Rather, the division failure in absence of MPS1 activity results from mitotic exit occurring before acentrosomal spindles can become bipolar. Furthermore, prevention of centrosome separation suffices to make cell division reliant on a SACdependent mitotic delay. Thus, centrosomes and their definition of two spindle poles early in mitosis provide a ‘timely twoness’ that allows cell division to occur in absence of a SACdependent mitotic delay.

 Cell Biology
Fertilization occurs before the completion of oocyte meiosis in the majority of animal species and sperm contents move long distances within the zygotes of mouse and C. elegans. If incorporated into the meiotic spindle, paternal chromosomes could be expelled into a polar body resulting in lethal monosomy. Through live imaging of fertilization in C. elegans, we found that the microtubule disassembling enzymes, katanin and kinesin13 limit longrange movement of sperm contents and that maternal ataxin2 maintains paternal DNA and paternal mitochondria as a cohesive unit that moves together. Depletion of katanin or double depletion of kinesin13 and ataxin2 resulted in the capture of the sperm contents by the meiotic spindle. Thus limiting movement of sperm contents and maintaining cohesion of sperm contents within the zygote both contribute to preventing premature interaction between maternal and paternal genomes.