Many epithelial tissues pack cells into a honeycomb pattern to support their structural and functional integrity. Developmental changes in cell packing geometry have been shown to be regulated by both mechanical and biochemical interactions between cells; however, it is largely unknown how molecular and cellular dynamics and tissue mechanics are orchestrated to realize the correct and robust development of hexagonal cell packing. Here, by combining mechanical and genetic perturbations along with live imaging and Bayesian force inference, we investigate how mechanical forces regulate cellular dynamics to attain a hexagonal cell configuration in the Drosophila pupal wing. We show that tissue stress is oriented towards the proximal-distal axis by extrinsic forces acting on the wing. Cells respond to tissue stretching and orient cell contact surfaces with the stretching direction of the tissue, thereby stabilizing the balance between the intrinsic cell junction tension and the extrinsic force at the cell-population level. Consequently, under topological constraints of the two-dimensional epithelial sheet, mismatches in the orientation of hexagonal arrays are suppressed, allowing more rapid relaxation to the hexagonal cell pattern. Thus, our results identify the mechanism through which the mechanical anisotropy in a tissue promotes ordering in cell packing geometry.

Epithelial and epithelial-derived tissues often pack cells into a honeycomb pattern to realize proper functioning of tissues by supporting structural strength, promoting the correct formation of cell polarity, etc. (Thompson, 1917; Wootton, 1992; Raphael and Altschuler, 2003; Classen et al., 2005; Cooper et al., 2008; Ma et al., 2008). Cell geometry (e.g. the eventual fraction of hexagonal cells) is highly reproducible in each epithelial tissue and is regulated by forces acting along the plane of the adherens junction, such as tension, which shortens a cell contact surface, and pressure, which counteracts tension to maintain the size of a cell (Fig. 1A) (Graner and Sawada, 1993; Mofrad and Kamm, 2006; Farhadifar et al., 2007; Käfer et al., 2007; Lecuit and Lenne, 2007; Honda et al., 2008; Rauzi et al., 2008; Aigouy et al., 2010; Dahmann et al., 2011; Eaton and Jülicher, 2011; Lecuit et al., 2011; Kasza and Zallen, 2011; Bosveld et al., 2012; Sano et al., 2012). Theoretically, developmental changes in cell geometry have been described as the minimization of a potential energy, by which the balance of cell pressures and cell junction tensions is considered (Honda, 1983; Graner and Sawada, 1993; Farhadifar et al., 2007; Käfer et al., 2007; Aegerter-Wilmsen et al., 2010). The relaxation process towards a hexagonal cell configuration is often slowed down by the existence of local energy-minimum states that satisfy the force balance, which suppresses reconnection of cell contact surfaces to change nonhexagonal cells into hexagonal ones. At present, little is understood about how such disorder in cell geometry is resolved and/or escaped to realize the correct and robust development of hexagonal cell packing. Identification of such a mechanism of cell packing could shed light on aspects of the physical basis of multicellular pattern formation.

The Drosophila pupal wing epithelium, which prominently increases the fraction of hexagonal cells during the course of its development, provides an excellent model system in which to study cell packing (Fig. 1B). In phase I [from ∼15-24 hours after puparium formation (APF)] of wing development, cells undergo one or two rounds of cell division, and the initial, nearly isotropic morphology of the cells becomes elongated along the proximal-distal (PD) axis of the wing (Fig. 1C) (Fristrom and Fristrom, 1993; Classen et al., 2005; Aigouy et al., 2010). In phase II (from ∼24-32 hours APF), the bias in the lengths of cell contact surfaces exhibits a moderate decrease. In addition, extensive cell rearrangement occurs, and the formation of new cell contact surfaces is biased towards the PD direction (Fig. 1D). Studies based on sophisticated image analyses have clarified the regulation of cell rearrangement by planar cell polarity (PCP) proteins (Strutt and Strutt, 2009), Rap1 GTPase (Knox and Brown, 2002) and tissue stretching along the PD axis by hinge constriction (Classen et al., 2005; Aigouy et al., 2010; O’Keefe et al., 2012). However, the mechanism by which tissue stretching directs cell rearrangement is yet to be clarified. Moreover, the link between directional cell rearrangement and hexagonal cell packing is missing. Answering these questions require space-time maps of force/stress with molecular and cellular quantifiers.

Various in vivo mechanical measurement methods have been developed. In particular, laser ablation of individual cell junctions is used as a tool to evaluate the tension acting on a cell contact surface. Such measurements have shed light on how cell shape and rearrangement are regulated by the activity and/or localization of force-generating molecular machinery within a cell (Hutson et al., 2003; Rauzi et al., 2008; Paluch and Heisenberg, 2009; Bosveld et al., 2012). A new complementary approach to such subcellular and invasive measurements is force inference from cell shape and connectivity, which offers cellular resolution and is both global and noninvasive (Chiou et al., 2012; Ishihara and Sugimura, 2012). In our previous study, we formulated a Bayesian framework of force inference, in which all the cell junction tensions, differences in pressures among cells, and tissue stress are simultaneously inferred from the observed geometry of cells, up to a scaling factor (supplementary material Appendix S1). We have shown that inferred force and stress values are consistent with those obtained using other methods, such as laser ablation of cortical actin cables, quantification of myosin concentration and photo-elasticity (Nienhaus et al., 2009), and large-scale tissue ablation (Bonnet et al., 2012; Ishihara and Sugimura, 2012; Ishihara et al., 2013). The global and noninvasive nature of the Bayesian force-inference method uniquely enables us to quantify space-time maps of force/stress in tissues and to relate the maps to hexagonal cell packing processes.

In this study, we address how forces and stress promote ordering of cell packing geometry in the Drosophila pupal wing. Force-inference analysis showed that tissue stress remained highly anisotropic until early phase II of wing development, when hexagonal cell packing occurred. We determined that anisotropic tissue stress promoted hexagonal cell packing by suppressing mismatches in the orientation of hexagonal arrays. Furthermore, our results indicated that a PCP protein, Flamingo (Fmi; Starry night, Stan - FlyBase), was indispensable for hexagonal cell packing after the balance between the extrinsic and intrinsic forces was nearly stabilized.

Drosophila strains

The flies used in this study were sqhp-sqhGFP (Royou et al., 2004), DEcadherin-GFP knock-in (Huang et al., 2009), UAS-Dα-catenin-TagRFP (Ishihara and Sugimura, 2012), tubP-gal80ts (Bloomington Stock Center #7108), UAS-sqh dsRNA (Bloomington Stock Center #33892) and UAS-flamingo dsRNA (Shimada et al., 2006). Flies were raised at 25°C unless otherwise noted.

Immunohistochemistry

Anti-DCAD2 (Shg - FlyBase) (Oda et al., 1994) and anti-Zipper (kindly provided by Fumio Matsuzaki, RIKEN CDB, Japan) antibodies were used in this study. Pupae at appropriate ages were dissected, and wings were fixed at room temperature for 30 minutes in PBS containing 4% paraformaldehyde. After washing with PBS containing 0.1% Triton X-100, these preparations were incubated with the antibodies indicated above.

Image collection and analysis

Preparation of samples of the Drosophila pupal wing and scutum for image collection was conducted as previously described (Shimada et al., 2006; Koto et al., 2009; Ishihara and Sugimura, 2012). Images were acquired using an inverted confocal microscope (FV1000D; Olympus) equipped with an Olympus 60×/NA1.2 SPlanApo water-immersion objective at 25°C unless otherwise noted. An inverted confocal microscope (A1R; Nikon) equipped with a 60×/NA1.2 Plan Apochromat water-immersion objective was also used (Fig. 2; Fig. 3E; Fig. 4; Fig. 7C,D; supplementary material Fig. S9). After imaging, we confirmed that the pupae survived to at least the pharate stage. Images were processed and analyzed as previously described (Ishihara and Sugimura, 2012).

Laser ablation experiment

To measure the relative value of cell junction tension, laser ablation of the single contact surface was performed as described previously (Ishihara and Sugimura, 2012). The displacement of the vertices at 16 seconds after laser irradiation was used to calculate Vmax of the vertices (the vertex displacement at 11 seconds was used in the previous study). To evaluate the anisotropy of global stress, a two-photon laser tuned to 720 nm wavelength (Chameleon; Coherent) was used to ablate a group of cells. Images were acquired using an inverted confocal microscope (SP5; Leica) equipped with a Leica 60×/NA1.2 HCX water-immersion objective. The displacement of cells at 14 seconds after laser irradiation was measured to calculate the difference between the velocity of cell groups along the PD and anterior-posterior (AP) axes (Vx - Vy).

Mechanical manipulation to relax tissue stretch

To test the possibility that anisotropic tissue stress plays an important role in hexagonal pattern formation, the wing was detached with the hinge by forceps at 23.5-24 hours APF, as previously described (Aigouy et al., 2010). We also employed a different method of relaxing tissue stretch: a two-photon laser tuned to 800 nm wavelength (Chamereon; Coherent) cut the wing along the anterior cross vein at 23.5-24 hours APF. A decrease in tension anisotropy was confirmed by laser cutting (data not shown). We adjusted the strength of laser irradiation to avoid damage to the wing. Only pupae that normally developed prehairs of wing cells were included in further analyses (100% of the samples cut by forceps and >60% of the samples cut by laser).

Statistics

Data are presented as mean ± s.d. P-values were calculated based on the Wilcoxon rank sum test. We characterized anisotropy in the inferred values for tension, the lengths of the cell contact surfaces, and the signal intensity of MRLC-GFP towards the PD axis using circular statistics (Fisher, 1993).

Force inference and numerical simulations

Detailed information on these methods is given in supplementary material Appendix S1.

The alignment of hexagonal cells proceeded over the same time course as the increase in hexagonal cells in the Drosophila wing

In the Drosophila pupal wing, the fraction of hexagons gradually decreased after cell division began (phase I; Fig. 1E). Then, the increase in hexagonal cells became prominent after the cells stopped dividing (phase II). We characterized developmental changes in hexagonal cell shape using the following quantities: (1) the orientation of hexagonal cells, which was defined by the average of cos6θj, where θ1-6 are the edge angles (⟨cos6θj⟩; supplementary material Fig. S1A) and (2) the orientation of cell shape anisotropy, which was defined by the angle of the longest axis of a fitted ellipse (cosφ; supplementary material Fig. S1B). We found that hexagonal cells increasingly pointed to the PD axis during phase II (more blue cells at later time points in Fig. 1F; Fig. 1G), whereas the orientation of cell shape anisotropy continued pointing to the PD axis from mid-phase I onwards (Fig. 1H). Thus, the PD alignment, but not the PD elongation of hexagonal cells, proceeded over the same time course during which the fraction of hexagonal cells increased in Drosophila wings. The mechanical process, which underlies the packing and alignment of hexagonal cells during phase II of Drosophila pupal development, is addressed below.

Fig. 1.

Packing and alignment of hexagonal cells in the Drosophila pupal wing. (A) The structure and forces of an epithelial sheet. Cortical actin-myosin cables (red) run along the plane of the adherens junctions. Tension shortens a cell contact surface (an edge). Pressure maintains the size of a cell. (B) Left: Schematic of an adult Drosophila. Right: A wide image of the pupal wing at 32 hours APF. As spatial landmarks, we used three sensory organs on the L3 vein (L3-1 to L3-3) and analyzed the intervein region between L3-1 and L3-3 (rectangle). In all images of the pupal wing, the vertical and horizontal directions are aligned with the anterior-posterior (AP) and proximal-distal (PD) axes, respectively. We set a line connecting L3-1 and L3-3 as the PD axis. (C) Cell-level dynamics in the Drosophila pupal wing (Aigouy et al., 2010). Wing cells at each stage are colored according to their number of edges (red, square; green, pentagon; gray, hexagon; blue, heptagon; magenta, octagon). (D) The average angular distribution of the newly formed edges during time-lapse imaging from 23.5 to 30.5 hours APF was quantified at 30.5 hours APF (n=3). (E) Developmental changes in the fraction of hexagonal cells in the wing (red) and scutum (blue). Mean ± s.d. are plotted. The number of flies examined is indicated. (F) Hexagonal cells are color-coded by ⟨cos6θj⟩ (the average of cosj, where θ1-6 represents the angle of the edges belonging to the cell; supplementary material Fig. S1A). PD- and AP-oriented hexagons are colored blue and yellow, respectively. Nonhexagonal cells are shown in black. Time-lapse imaging was conducted from 23.5 to 30.5 hours APF. (G,H) The average values of ⟨cos6θj⟩ (G) (supplementary material Fig. S1A) and cosφ (H) (supplementary material Fig. S1B) in each wing or scutum were calculated, and the mean ± s.d. among samples at each developmental stage is plotted. The number of control flies examined is the same as in E. The fly genotype was sqhp-sqhGFP, apterous-gal4/sqhp-sqhGFP, UAS-Dα-catenin-TagRFP.

Fig. 1.

Packing and alignment of hexagonal cells in the Drosophila pupal wing. (A) The structure and forces of an epithelial sheet. Cortical actin-myosin cables (red) run along the plane of the adherens junctions. Tension shortens a cell contact surface (an edge). Pressure maintains the size of a cell. (B) Left: Schematic of an adult Drosophila. Right: A wide image of the pupal wing at 32 hours APF. As spatial landmarks, we used three sensory organs on the L3 vein (L3-1 to L3-3) and analyzed the intervein region between L3-1 and L3-3 (rectangle). In all images of the pupal wing, the vertical and horizontal directions are aligned with the anterior-posterior (AP) and proximal-distal (PD) axes, respectively. We set a line connecting L3-1 and L3-3 as the PD axis. (C) Cell-level dynamics in the Drosophila pupal wing (Aigouy et al., 2010). Wing cells at each stage are colored according to their number of edges (red, square; green, pentagon; gray, hexagon; blue, heptagon; magenta, octagon). (D) The average angular distribution of the newly formed edges during time-lapse imaging from 23.5 to 30.5 hours APF was quantified at 30.5 hours APF (n=3). (E) Developmental changes in the fraction of hexagonal cells in the wing (red) and scutum (blue). Mean ± s.d. are plotted. The number of flies examined is indicated. (F) Hexagonal cells are color-coded by ⟨cos6θj⟩ (the average of cosj, where θ1-6 represents the angle of the edges belonging to the cell; supplementary material Fig. S1A). PD- and AP-oriented hexagons are colored blue and yellow, respectively. Nonhexagonal cells are shown in black. Time-lapse imaging was conducted from 23.5 to 30.5 hours APF. (G,H) The average values of ⟨cos6θj⟩ (G) (supplementary material Fig. S1A) and cosφ (H) (supplementary material Fig. S1B) in each wing or scutum were calculated, and the mean ± s.d. among samples at each developmental stage is plotted. The number of control flies examined is the same as in E. The fly genotype was sqhp-sqhGFP, apterous-gal4/sqhp-sqhGFP, UAS-Dα-catenin-TagRFP.

Tissue stress remained highly anisotropic until early phase II of wing development, when hexagonal cell packing occurred

We started by evaluating maps of cell junction tension and tissue stress in the wing at several developmental stages (supplementary material Figs S2-S6) and by analyzing temporal changes in tension/stress maps from time-lapse data in phase II (Fig. 2). Force-inference analysis showed that distribution of the inferred tension was biased at the beginning of phase II (24 hours APF). Indeed, cell contact surfaces with higher and lower tension values were primarily located along the PD and AP directions, respectively (in Fig. 2A, arrows and arrowheads, respectively; see supplementary material Fig. S5 for results of laser cutting of a single cell contact surface). The myosin regulatory light chain-green fluorescent protein (MRLC-GFP) signal was more enriched on the PD edges than on the AP edges [in Fig. 2A′, yellow and blue arrowheads, respectively; see supplementary material Fig. S2F for the enrichment of myosin heavy chain protein (Kiehart et al., 1989) on the PD edges]. Thus, all of our data indicated that cell junction tension was stronger on the PD edges at 24 hours APF (see supplementary material Fig. S2A-E for patterns of inferred tension and myosin at each stage of wing development). We characterized the anisotropies of inferred tension, edge length and myosin towards the PD axis in each wing, where larger R-values indicated larger anisotropy (see supplementary material Fig. S3A for the definition of R-values). The quantification of R-values from time-lapse data indicated that the anisotropies of the inferred tension, cell shape and myosin gradually decreased during phase II (Fig. 2D; see supplementary material Fig. S3B and Fig. S4 for the average of R-values at each stage of wing development and their statistical tests).

Fig. 2.

Stress in the wing was stronger along the PD axis, when hexagonal cell packing occurred. (A-C′) Time-lapse analysis of the control wing from 24 to 30 hours APF (n=5). The inferred tension (A-C), MRLC-GFP (A′-C′), inferred local cell stress (A′-C′) and inferred tissue stress (insets in A′-C′; the longest axis of the stress ellipse represents the maximum stress direction) are shown. The PD and AP edges are indicated with arrows and arrowheads (A) and with yellow and blue arrowheads (A′), respectively. (D-D′) The anisotropy of the inferred tension (RT; D), the length of the edge (RL; D′) and the MRLC-GFP signal intensity (RS; D′) are plotted (different colors represent different samples). (E) Schematics of global (tissue) stress and local cell stress (supplementary material Appendix S1) (Ishihara and Sugimura, 2012). The anisotropy of tissue stress was determined by geometrical and tension anisotropies (i.e. directional bias in the distribution of cell contact surfaces and that in tension values). This means that tissue stress can be anisotropic when tension is constant for all cell contact surfaces. (F) Tissue stress anisotropy was measured from global ablation data in the wing and scutum (supplementary material Movies 1, 2). The difference between the velocity of cell groups along the horizontal and vertical axes (Vx - Vy) was calculated from the displacement of cells at 14 seconds after laser irradiation, and the mean ± s.d. among samples at each developmental stage is plotted. Fly genotypes are sqhp-sqhGFP, apterous-gal4/sqhp-sqhGFP, UAS-Dα-catenin-TagRFP (A-D) and DECad-GFP knock-in (F). Scale bar: 20 μm.

Fig. 2.

Stress in the wing was stronger along the PD axis, when hexagonal cell packing occurred. (A-C′) Time-lapse analysis of the control wing from 24 to 30 hours APF (n=5). The inferred tension (A-C), MRLC-GFP (A′-C′), inferred local cell stress (A′-C′) and inferred tissue stress (insets in A′-C′; the longest axis of the stress ellipse represents the maximum stress direction) are shown. The PD and AP edges are indicated with arrows and arrowheads (A) and with yellow and blue arrowheads (A′), respectively. (D-D′) The anisotropy of the inferred tension (RT; D), the length of the edge (RL; D′) and the MRLC-GFP signal intensity (RS; D′) are plotted (different colors represent different samples). (E) Schematics of global (tissue) stress and local cell stress (supplementary material Appendix S1) (Ishihara and Sugimura, 2012). The anisotropy of tissue stress was determined by geometrical and tension anisotropies (i.e. directional bias in the distribution of cell contact surfaces and that in tension values). This means that tissue stress can be anisotropic when tension is constant for all cell contact surfaces. (F) Tissue stress anisotropy was measured from global ablation data in the wing and scutum (supplementary material Movies 1, 2). The difference between the velocity of cell groups along the horizontal and vertical axes (Vx - Vy) was calculated from the displacement of cells at 14 seconds after laser irradiation, and the mean ± s.d. among samples at each developmental stage is plotted. Fly genotypes are sqhp-sqhGFP, apterous-gal4/sqhp-sqhGFP, UAS-Dα-catenin-TagRFP (A-D) and DECad-GFP knock-in (F). Scale bar: 20 μm.

We next evaluated tissue stress and local cell stress tensors by integrating inferred tensions and pressures of individual cells (Fig. 2E; supplementary material Appendix S1) (Batchelor, 1970; Ishihara and Sugimura, 2012), where the maximum stress direction along which the stress was strongest in a cell/tissue was indicated by the longest axis of stress ellipses. Our force-inference analysis showed that the maximum stress direction of a cell/tissue pointed to the PD axis throughout phase II (Fig. 2A′-C′, magenta crosses, and blue crosses in insets, respectively). To examine experimentally tissue stress anisotropy, the group of cells shown in supplementary material Movie 1 was ablated, and the ablation contour was tracked; the anisotropic extension of the cell group contour confirmed the stronger tissue stress along the PD axis of the wing.

We next investigated whether tissue stress anisotropy decayed during phase II, as cell junction tension anisotropy did. The magnitude of the tissue stress anisotropy was quantified by measuring the difference between the extension of the cell group contour along the PD and AP axes from large-scale tissue ablation data (e.g. supplementary material Movie 1) and by inferring the normal stress difference σA ≡ (σxx - σyy)/2, which was used previously (Ishihara et al., 2013) to cross-validate the Bayesian force inference with large-scale tissue ablation (Bonnet et al., 2012). Both measurements indicated that tissue stress remained highly anisotropic until early phase II (∼26-27 hours APF; Fig. 2F; supplementary material Fig. S6) when tension and myosin anisotropies had already started decreasing (Fig. 2A-D; supplementary material Fig. S2 and Fig. S3B).

Note that tissue stress anisotropy was determined by the directional bias both in the distribution of cell contact surfaces and in cell junction tension values (i.e. geometrical anisotropy and tension anisotropy, respectively; Fig. 2E). Accordingly, the discrepancy in the timing of the decrease of anisotropies of tissue stress and cell junction tension can be explained as follows: the decrease in tension anisotropy was compensated for by the increase in PD cell contact surfaces by cell rearrangement, thereby maintaining the anisotropy of tissue stress at high levels until early phase II, when hexagonal cell packing occurred.

Anisotropic tissue stress was required and sufficient for promotion of hexagonal cell packing

Next, we examined the role of anisotropic tissue stress in promoting the formation of a hexagonal cell array. To highlight the importance of anisotropic tissue stress in controlling the polygonal distribution of cells, we compared cell geometry and force/stress maps of the wing with those of the scutum (anterior part of the dorsal thorax; Fig. 1B). We previously reported that tissue stress in the scutum is greater along the AP axis (supplementary material Fig. S3B-B′) (Ishihara and Sugimura, 2012); however, compared with the wing, the amplitude of the tissue stress anisotropy is weak (compare supplementary material Movies 1 and 2; Fig. 2F). The final fraction of hexagonal cells in the scutum was close to that in the wing at 13 hours APF when the tissue stress was isotropic (∼60%; Fig. 1E), whereas the fraction of hexagons increased to >70% in the wing (Fig. 1E). In addition, the alignment of hexagonal cells to the tissue stretch axis was suppressed in the scutum (Fig. 1G). These observations suggested the possibility that strong anisotropic stress in tissue supported efficient hexagonal packing. Consistent with this hypothesis, a decrease in the number of hexagonal cells was observed by severing the wing from the hinge at an earlier developmental stage (Aigouy et al., 2010); however, further studies are necessary because pleiotropic effects, such as alternation of cell division patterns, are also induced by this technique (Aigouy et al., 2010; data not shown). We thus cut the wing at 23.5-24 hours APF to exclude effects on cell division. This mechanical manipulation significantly decreased the anisotropy of tissue stress (Fig. 3A) and resulted in a smaller fraction of hexagonal cells (Fig. 3B-D; supplementary material Fig. S7A-C). Moreover, the orientation of hexagonal cells was more disperse in the cut wings (Fig. 3C, circles label some of the AP-oriented cells), indicating that the orientation of hexagonal cells was improperly assigned (Fig. 3E,F; supplementary material Fig. S7D). Strongly anisotropic tissue stress was therefore required for the packing and alignment of hexagonal cells.

Fig. 3.

Anisotropic tissue stress promoted hexagonal packing. (A) The inferred normal stress difference σA ≡ (σxx - σyy)/2 is plotted for control wings and cut wings severed by forceps at 23.5-24 hours APF (different colors represent different samples). (B,C) Images of control (B) and cut (C) wings. Cells are colored according to the number of edges. The percentage of hexagonal cells was: 72.8% (B) and 59.6% (C). Circles in C label some of the AP-oriented cells. (D) Quantification of the fraction of hexagonal cells in control wings (as in B; left bar) and wings cut by forceps (as in C; right bar). (E) Quantification of cell rearrangement in the control and cut wings (n=4 and n=3, respectively). The angles of newly formed cell contact surfaces are measured from time-lapse data taken at 24-25 hours APF. The direction of each edge is classified as illustrated by the semicircle (e.g. red class I for the PD edges and blue class III for the AP edges). (F) The average value of ⟨cos6θj⟩ in each wing or wing cut at 23.5-24 hours APF was calculated, and the mean ± s.d. among samples at each developmental stage is plotted. (G) Diagram of the relationship between potential energy and cell configuration. (H) Patterns of cells obtained by the numerical simulation of cell rearrangement without stretch, with isotropic stretch and with uniaxial (horizontal) stretch. t, time in numerical simulation steps. (I) Time evolution of the fraction of hexagonal cells in numerical simulation. The time scale and magnitude of noise are controlled by τ and z, respectively (1/τ=0.1 and z=30%; see supplementary material Appendix S1, equations S2-S5). (J) Fractions of hexagonal cells at t=2000 in the numerical simulations are plotted for several values of and (1/τ=0.1 and z=30%; n=8 for each set of parameters). and control the strength of the tension. The results of the simulations under isotropic stretch and horizontal stretch are indicated with crosses and circles, respectively. In D-F, the numbers of flies examined are indicated at the top. The fly genotype is as described in the legend of Fig. 1. Scale bar: 20 μm.

Fig. 3.

Anisotropic tissue stress promoted hexagonal packing. (A) The inferred normal stress difference σA ≡ (σxx - σyy)/2 is plotted for control wings and cut wings severed by forceps at 23.5-24 hours APF (different colors represent different samples). (B,C) Images of control (B) and cut (C) wings. Cells are colored according to the number of edges. The percentage of hexagonal cells was: 72.8% (B) and 59.6% (C). Circles in C label some of the AP-oriented cells. (D) Quantification of the fraction of hexagonal cells in control wings (as in B; left bar) and wings cut by forceps (as in C; right bar). (E) Quantification of cell rearrangement in the control and cut wings (n=4 and n=3, respectively). The angles of newly formed cell contact surfaces are measured from time-lapse data taken at 24-25 hours APF. The direction of each edge is classified as illustrated by the semicircle (e.g. red class I for the PD edges and blue class III for the AP edges). (F) The average value of ⟨cos6θj⟩ in each wing or wing cut at 23.5-24 hours APF was calculated, and the mean ± s.d. among samples at each developmental stage is plotted. (G) Diagram of the relationship between potential energy and cell configuration. (H) Patterns of cells obtained by the numerical simulation of cell rearrangement without stretch, with isotropic stretch and with uniaxial (horizontal) stretch. t, time in numerical simulation steps. (I) Time evolution of the fraction of hexagonal cells in numerical simulation. The time scale and magnitude of noise are controlled by τ and z, respectively (1/τ=0.1 and z=30%; see supplementary material Appendix S1, equations S2-S5). (J) Fractions of hexagonal cells at t=2000 in the numerical simulations are plotted for several values of and (1/τ=0.1 and z=30%; n=8 for each set of parameters). and control the strength of the tension. The results of the simulations under isotropic stretch and horizontal stretch are indicated with crosses and circles, respectively. In D-F, the numbers of flies examined are indicated at the top. The fly genotype is as described in the legend of Fig. 1. Scale bar: 20 μm.

To determine whether anisotropic tissue stress was sufficient for promoting hexagonal packing, we performed numerical simulation of cell rearrangement under three conditions: no stretch, isotropic stretch and horizontal stretch (Fig. 3G-J; supplementary material Fig. S8; note that experimentally applying stretch to the pupal wing is technically very difficult). For numerical simulations, we modified the cell vertex model to control the anisotropic stress environment (supplementary material Appendix S1) (Honda, 1983; Graner and Sawada, 1993; Ouchi et al., 2003; Farhadifar et al., 2007; Käfer et al., 2007; Rauzi et al., 2008). Note that the hexagonal configuration was probably the energy-minimum state in each of the three conditions (Hales, 2001) (Fig. 3G), but local energy-minimum states hindered relaxation to the hexagonal state. Fluctuation could help the relaxation process by preventing a system from being trapped in a local minimum. We therefore introduced noise to the system by varying parameters in the line tension term [in the wing, estimated tensions fluctuate about ±10% in 10 minutes (Sugimura et al., 2013)]. We found that unless the simulations were conducted under large noise with a slow time correlation, the eventual fraction of hexagonal cells under horizontal stretch was significantly larger than that under isotropic stretch (Fig. 3H-J; supplementary material Fig. S8B,C and Appendix S1 for results about dependence on noise parameters). These data indicated that horizontal stretch could drive an increase in hexagons, whereas the fluctuation in forces led to packing of cells into a hexagonal pattern in other conditions.

Cells rearranged cell contact surfaces to better align local cell stress and cell junction tension with global tissue stress

To understand the mechanism by which anisotropic tissue stress realized a rapid relaxation towards a hexagonal configuration, we addressed how cells responded to tissue stretch and changed cell packing geometry. For this, we manipulated extrinsic and intrinsic forces acting on the wing and monitored cellular responses to changes in the stress field in the tissue. Upon relaxation of tissue stress by detaching the wing from the hinge at 24 hours APF (Fig. 3A), cell shape anisotropy diminished, and cells decreased the degree of anisotropy of inferred tension, local cell stress, and myosin (compare Fig. 4A-D with Fig. 2A-D), indicating that the mechanical anisotropy of a cell was dependent on the tissue stretch. By contrast, when the cell junction tension was weakened by RNA interference (RNAi) against the spaghetti squash (sqh) gene, which encodes MRLC, excess shear deformation of cells was observed (supplementary material Fig. S9A,B). These results suggested that cells redistributed myosin and generated stronger intrinsic cell junction tension on the PD edges in order to resist tissue stretch that was oriented towards the PD axis by extrinsic force (supplementary material Fig. S9C). This response of wing cells to balance the extrinsic force began in phase I, when the hinge constriction begins, because cutting the wing in phase I also diminished tension and myosin anisotropies (data not shown).

Fig. 4.

The pattern of local cell stresses/cell junction tensions was dependent on anisotropic tissue stress. (A-C′) A pupal wing was detached from the hinge region by forceps at 23.5-24 hours APF and was observed at 24, 27 and 30 hours APF (n=3; the same set of samples was analyzed as in Fig. 3A). The inferred tension (A-C), MRLC-GFP (A′-C′), inferred local cell stress (A′-C′) and inferred tissue stress (insets in A′-C′; the longest axis of the stress ellipse represents the maximum stress direction) are shown. (D-D′) The anisotropy of the inferred tension (RT; D), the length of the edge (RL; D′) and the MRLC-GFP signal intensity (RS; D′) are plotted (different colors represent different samples). Gray lines are the corresponding results for control wings (Fig. 2D-D′). Similar results were obtained by cutting a wing at 15 hours APF (n=11). The fly genotype was sqhp-sqhGFP, apterous-gal4/sqhp-sqhGFP, UAS-Dα-catenin-TagRFP. Scale bar: 20 μm.

Fig. 4.

The pattern of local cell stresses/cell junction tensions was dependent on anisotropic tissue stress. (A-C′) A pupal wing was detached from the hinge region by forceps at 23.5-24 hours APF and was observed at 24, 27 and 30 hours APF (n=3; the same set of samples was analyzed as in Fig. 3A). The inferred tension (A-C), MRLC-GFP (A′-C′), inferred local cell stress (A′-C′) and inferred tissue stress (insets in A′-C′; the longest axis of the stress ellipse represents the maximum stress direction) are shown. (D-D′) The anisotropy of the inferred tension (RT; D), the length of the edge (RL; D′) and the MRLC-GFP signal intensity (RS; D′) are plotted (different colors represent different samples). Gray lines are the corresponding results for control wings (Fig. 2D-D′). Similar results were obtained by cutting a wing at 15 hours APF (n=11). The fly genotype was sqhp-sqhGFP, apterous-gal4/sqhp-sqhGFP, UAS-Dα-catenin-TagRFP. Scale bar: 20 μm.

Next, we tracked force/stress and cell geometry from time-lapse data in phase II, when hexagonal cell packing proceeds, and examined how cells alter their force/stress and geometry to achieve balance against tissue stretch. Specifically, we evaluated the alignment and magnitude of cell junction tension/local cell stress with respect to tissue stress before and after the cell rearrangement/alignment process, keeping in mind that tissue stress was oriented towards the PD axis by extrinsic force (Fig. 5A). Force-inference analysis showed that the stress of individual cells deviated significantly at ∼24 hours APF; some cells developed strongly anisotropic stress, whereas stress of others was almost isotropic (Fig. 5B). As hexagonal cells are pointing to the PD axis (Fig. 5C, horizontal axis), local cell stress exhibited similar anisotropy and magnitude along the direction of maximum global stress (Fig. 5B-E; supplementary material Fig. S10A,B). Over the same time course, the magnitude of inferred tension became uniform (Fig. 5F), which possibly resulted from the fact that the increase of PD cell contact surfaces through cell rearrangement led to the decrease of mechanical load on each PD cell contact surface. These results showed that cells responded to anisotropic tissue stretch by orienting cell contact surfaces to align local cell stress/cell junction tension with the stretching direction of the tissue, thereby stabilizing the balance of intrinsic and extrinsic forces at the cell-population level (supplementary material Fig. S10C). Thus, global force balance assigned the orientation of wing cells, which concomitantly proceeded with hexagonal cell packing in the wing.

Fig. 5.

The alignment of local cell stresses/cell junction tensions with global tissue stress during hexagonal cell packing. (A) The anisotropy of local cell stress and its magnitude along the maximum stress direction of the tissue were calculated by a - b and c for each cell, respectively. (B-F) Developmental changes in inferred patterns of local cell stress and cell junction tension. Time-lapse data obtained from 23.5 to 30.5 hours APF were extracted at 12-minute intervals. (B) Images showing local stress of each cell, in which the crossed lines are the longest and shortest axes of the local stress ellipse of each cell. Initially, some cells exerted strong local stress along the PD axis, whereas others generated relatively weak and/or isotropic local stress. Then, the variance of local cell stress decreased, and individual cells eventually generated forces with similar anisotropy and magnitude along the maximum stress direction of the cell population (i.e. the PD direction). (C) The magnitude of local cell stress along the maximum stress direction of the tissue (c in A) for each hexagonal cell is plotted against its ⟨cos6θj⟩. (D-F) Temporal changes in the standard deviation of the magnitude of local cell stress along the maximum stress direction of the tissue (D); in the standard deviation of the anisotropy of the local stress (E) (a - b in A); and in the standard deviation of inferred tension are plotted (different colors represent different samples). The fly genotype was sqhp-sqhGFP, apterous-gal4/sqhp-sqhGFP, UAS-Dα-catenin-TagRFP. Scale bar: 20 μm.

Fig. 5.

The alignment of local cell stresses/cell junction tensions with global tissue stress during hexagonal cell packing. (A) The anisotropy of local cell stress and its magnitude along the maximum stress direction of the tissue were calculated by a - b and c for each cell, respectively. (B-F) Developmental changes in inferred patterns of local cell stress and cell junction tension. Time-lapse data obtained from 23.5 to 30.5 hours APF were extracted at 12-minute intervals. (B) Images showing local stress of each cell, in which the crossed lines are the longest and shortest axes of the local stress ellipse of each cell. Initially, some cells exerted strong local stress along the PD axis, whereas others generated relatively weak and/or isotropic local stress. Then, the variance of local cell stress decreased, and individual cells eventually generated forces with similar anisotropy and magnitude along the maximum stress direction of the cell population (i.e. the PD direction). (C) The magnitude of local cell stress along the maximum stress direction of the tissue (c in A) for each hexagonal cell is plotted against its ⟨cos6θj⟩. (D-F) Temporal changes in the standard deviation of the magnitude of local cell stress along the maximum stress direction of the tissue (D); in the standard deviation of the anisotropy of the local stress (E) (a - b in A); and in the standard deviation of inferred tension are plotted (different colors represent different samples). The fly genotype was sqhp-sqhGFP, apterous-gal4/sqhp-sqhGFP, UAS-Dα-catenin-TagRFP. Scale bar: 20 μm.

Comparison of extrinsic tissue stretch-driven mechanisms of cell rearrangement and intrinsic cell junction tension-driven mechanisms for promoting hexagonal cell packing

To support further the role of anisotropic tissue stress in controlling cell packing geometry, we numerically compared the mechanisms of cell rearrangement driven by extrinsic tissue stretch with those driven by intrinsic cell junction tension (Kasza and Zallen, 2011; Lecuit et al., 2011) for promoting hexagonal cell packing. For this, we performed two numerical simulations: (1) assuming that applied stretch was anisotropic, with no angular bias in parameters of cell junction tension (Fig. 6A; as was performed in Fig. 3H), and (2) vice versa (Fig. 6A′; supplementary material Appendix S1). Other parameters were the same in both numerical simulations. As mentioned before, the fraction of hexagonal cells increased in response to horizontal stretch (Fig. 6B; Fig. 3I). By contrast, fewer hexagonal cells were found in the biased tension parameter model (Fig. 6B′). The average of ⟨cos6θj⟩ increased by horizontal stretch (Fig. 6C), as was observed in the wing (Fig. 1G), but not in the biased tension parameter model (Fig. 6C′). Cells were relatively elongated along the horizontal direction in both the anisotropic tissue stretch model and the biased tension parameter model (Fig. 6D,D′). However, the angular dependence of tension was opposite between them; horizontal cell contact surfaces had larger tension than did vertical cell contact surfaces owing to the extrinsic tissue stretch in the former (intrinsic tension parameters were the same for all edges). By contrast, by definition of the model, horizontal cell contact surfaces had lower tension in the latter. These results indicated that numerical simulation with anisotropic tissue stretch agreed better with the experimental observations in terms of cell geometry and force.

Fig. 6.

Comparison between two mechanisms of cell rearrangement in hexagonal cell packing by numerical simulation. (A,A′) Numerical simulation of cell rearrangements. (A) Tissue stretch model. Isotropic or anisotropic tissue stretch by external forces was applied (as was performed in Fig. 3H). There was no angular bias in tension parameters. (A′) No external force was applied. σ, a parameter to control the strength of tension, depended on the angle of cell contact surface, where the magnitude of the angular bias was controlled by μ. Here, vertical cell contact surfaces were set to generate larger tension than horizontal ones. Details of these numerical simulations are described in supplementary material Appendix S1. (B,B′) Fractions of hexagonal cells at t=2000 according to numerical simulations of the tissue stretch model (B; see also Fig. 3I) and the biased tension parameter model (B′). The eventual cell configuration is shown in the right panel for μ=0.8 and =0.04. (C,C′) The mean ± s.d. of ⟨cos6θj⟩ among samples at each condition is plotted for the tissue stretch model (C) and the biased tension parameter model (C′). (D,D′) The anisotropy of the length of the edge (RL) is plotted for the tissue stretch model (D) and the biased tension parameter model (D′). n=10 for each condition (B-D).

Fig. 6.

Comparison between two mechanisms of cell rearrangement in hexagonal cell packing by numerical simulation. (A,A′) Numerical simulation of cell rearrangements. (A) Tissue stretch model. Isotropic or anisotropic tissue stretch by external forces was applied (as was performed in Fig. 3H). There was no angular bias in tension parameters. (A′) No external force was applied. σ, a parameter to control the strength of tension, depended on the angle of cell contact surface, where the magnitude of the angular bias was controlled by μ. Here, vertical cell contact surfaces were set to generate larger tension than horizontal ones. Details of these numerical simulations are described in supplementary material Appendix S1. (B,B′) Fractions of hexagonal cells at t=2000 according to numerical simulations of the tissue stretch model (B; see also Fig. 3I) and the biased tension parameter model (B′). The eventual cell configuration is shown in the right panel for μ=0.8 and =0.04. (C,C′) The mean ± s.d. of ⟨cos6θj⟩ among samples at each condition is plotted for the tissue stretch model (C) and the biased tension parameter model (C′). (D,D′) The anisotropy of the length of the edge (RL) is plotted for the tissue stretch model (D) and the biased tension parameter model (D′). n=10 for each condition (B-D).

In the wing, time-lapse data of cell rearrangement indicated that the inferred tension of shrinking AP edges was not as large as that of newly formed PD edges (Fig. 7A,B; supplementary material Movie 3). The level of MRLC-GFP signal intensity remained almost constant in newly formed PD edges [Fig. 7C,D; the ratio of MRLC-GFP signal intensity between shrinking AP edges and newly formed PD edges was 0.99±0.15 (n=22)]. Moreover, sqh RNAi in wing cells did not affect the PD bias in cell rearrangement [86% of the newly formed edges belonged to class I (see Fig. 3E for the classification of edges) (n=3); supplementary material Fig. S9 legend]. These results implied that intrinsic cell junction tension of shrinking AP edges did not primarily determine the direction of cell rearrangement in the wing. Taken together, both our computational and our experimental data further supported the role of anisotropic tissue stress in promoting hexagonal cell packing through rearrangement and alignment of cells.

Fig. 7.

Dynamics of inferred tension and myosin during the PD intercalation of cells. (A-B) Tension of the cells was inferred from the time-lapse data from 23.5 to 30.5 hours APF. (A,A′) Patterns of the inferred tension. Asterisks indicate the corresponding cells in each frame. Arrowheads and arrows point to the AP and PD edges, respectively. (B) Temporal changes in the inferred tension of the shrinking AP edge and the newly formed PD edges in the wing shown in supplementary material Movie 3. The PD intercalation of cells occurred between -3 and 0 minutes. When an edge underwent multiple rounds of remodeling, data were divided into single series of the PD intercalation. (C,D) MRLC-GFP during junctional remodeling was tracked from the time-lapse data from 24 to 25 hours APF. (C) Time-lapse images of MRLC-GFP (left) and Dα-catenin-TagRFP (right). Arrowheads and arrows point to the AP and PD edges, respectively. (D) The ratio of the MRLC-GFP signal intensity between shrinking AP edges and newly formed PD edges. The MRLC-GFP signal intensity at 1.5, 1 and 0.5 minutes before PD cell rearrangement was measured and its average at the three time-points was calculated. The average of MRLC-GFP signal intensity at 0, 0.5 and 1 minute after PD cell rearrangement was calculated, and was divided by the value before junctional exchange. The fly genotype was sqhp-sqhGFP, apterous-gal4/sqhp-sqhGFP, UAS-Dα-catenin-TagRFP. Scale bars: 5 μm.

Fig. 7.

Dynamics of inferred tension and myosin during the PD intercalation of cells. (A-B) Tension of the cells was inferred from the time-lapse data from 23.5 to 30.5 hours APF. (A,A′) Patterns of the inferred tension. Asterisks indicate the corresponding cells in each frame. Arrowheads and arrows point to the AP and PD edges, respectively. (B) Temporal changes in the inferred tension of the shrinking AP edge and the newly formed PD edges in the wing shown in supplementary material Movie 3. The PD intercalation of cells occurred between -3 and 0 minutes. When an edge underwent multiple rounds of remodeling, data were divided into single series of the PD intercalation. (C,D) MRLC-GFP during junctional remodeling was tracked from the time-lapse data from 24 to 25 hours APF. (C) Time-lapse images of MRLC-GFP (left) and Dα-catenin-TagRFP (right). Arrowheads and arrows point to the AP and PD edges, respectively. (D) The ratio of the MRLC-GFP signal intensity between shrinking AP edges and newly formed PD edges. The MRLC-GFP signal intensity at 1.5, 1 and 0.5 minutes before PD cell rearrangement was measured and its average at the three time-points was calculated. The average of MRLC-GFP signal intensity at 0, 0.5 and 1 minute after PD cell rearrangement was calculated, and was divided by the value before junctional exchange. The fly genotype was sqhp-sqhGFP, apterous-gal4/sqhp-sqhGFP, UAS-Dα-catenin-TagRFP. Scale bars: 5 μm.

A physical mechanism of hexagonal cell packing

To connect the mechanical force balance with cell packing geometry, we sought to determine how the orientation of hexagonal cells, which we showed was assigned by the mechanical force balance in the wing, led to efficient hexagonal cell packing. Our analysis of cell configurations clarified two characteristics of wing cells: (1) a spatial map of ⟨cos6θj⟩ indicated that hexagonal cells re-oriented to match their neighbors’ orientation better (hexagonal array; Fig. 8A,B), and (2) nonhexagonal cells appeared at the boundaries of hexagonal arrays that point to different orientations (i.e. segregation of the hexagonal array; Fig. 8A). We characterized segregation of the hexagonal array by quantifying the difference between the largest and smallest ⟨cos6θj⟩ among adjacent cells (denoted as Di; cumulative distributions are plotted in Fig. 8C). In Fig. 8C, the rightward shift of the blue lines at later developmental stages indicated that nonhexagonal cells present between the mismatched hexagonal cell arrays were more predominant after the cells were rearranged during hexagonal cell packing.

Fig. 8.

Changes in the configuration of cells during hexagonal cell packing. (A) Images of control and cut (as in supplementary material Fig. S7) wings. Each hexagonal cell is colored according to ⟨cos6θj⟩. Nonhexagonal cells are shown in black. (B,C) Quantitative analysis of the segregation of hexagonal cell arrays. For clarity, hi denotes ⟨cos6θj⟩ of the ith cell. (B) Developmental changes in correlation coefficients of hi (i.e. ⟨cos6θj⟩) between neighboring hexagonal cells. The correlation coefficients were computed by C = Σ[ij] (hihj - 2)/NHH2, where Σ[ij] indicates taking the sum over adjacent pairs of hexagonal cells, and NHH is the number of the pairs. The correlation coefficient increased up to ∼0.7, indicating local alignment of hexagonal cell orientation (i.e. formation of a hexagonal cell array). The number of control flies examined was the same as in Fig. 1E. (C) The difference between the largest and smallest hi (i.e. ⟨cos6θj⟩) among adjacent cells is denoted as Di. Note that Di is defined for both hexagonal and nonhexagonal cells. Cumulative frequencies of Di for hexagonal (red) and nonhexagonal (blue) cells in numerical simulation or at the stages indicated are shown. The values of parameters in numerical simulation were the same as those employed in Fig. 3H. From 21 to 23 hours APF, the distributions of Di were almost identical between hexagonal and nonhexagonal cells (P=0.21, Wilcoxon rank sum test). During the later stages, Di for nonhexagonal cells was increased whereas that for hexagonal cells was not (P<10-8). The fly genotype was sqhp-sqhGFP, apterous-gal4/sqhp-sqhGFP, UAS-Dα-catenin-TagRFP. Scale bar: 20 μm.

Fig. 8.

Changes in the configuration of cells during hexagonal cell packing. (A) Images of control and cut (as in supplementary material Fig. S7) wings. Each hexagonal cell is colored according to ⟨cos6θj⟩. Nonhexagonal cells are shown in black. (B,C) Quantitative analysis of the segregation of hexagonal cell arrays. For clarity, hi denotes ⟨cos6θj⟩ of the ith cell. (B) Developmental changes in correlation coefficients of hi (i.e. ⟨cos6θj⟩) between neighboring hexagonal cells. The correlation coefficients were computed by C = Σ[ij] (hihj - 2)/NHH2, where Σ[ij] indicates taking the sum over adjacent pairs of hexagonal cells, and NHH is the number of the pairs. The correlation coefficient increased up to ∼0.7, indicating local alignment of hexagonal cell orientation (i.e. formation of a hexagonal cell array). The number of control flies examined was the same as in Fig. 1E. (C) The difference between the largest and smallest hi (i.e. ⟨cos6θj⟩) among adjacent cells is denoted as Di. Note that Di is defined for both hexagonal and nonhexagonal cells. Cumulative frequencies of Di for hexagonal (red) and nonhexagonal (blue) cells in numerical simulation or at the stages indicated are shown. The values of parameters in numerical simulation were the same as those employed in Fig. 3H. From 21 to 23 hours APF, the distributions of Di were almost identical between hexagonal and nonhexagonal cells (P=0.21, Wilcoxon rank sum test). During the later stages, Di for nonhexagonal cells was increased whereas that for hexagonal cells was not (P<10-8). The fly genotype was sqhp-sqhGFP, apterous-gal4/sqhp-sqhGFP, UAS-Dα-catenin-TagRFP. Scale bar: 20 μm.

Based on these results, we propose a physical mechanism of hexagonal cell packing. In an isotropic stress field, local relaxation to hexagonal cells often causes a mismatch in the orientation of hexagonal arrays (e.g. white arrows in Fig. 9A point to different orientations in top and bottom hexagonal arrays); thus, nonhexagonal cells appear at the boundaries of the arrays because of topological constraints. By contrast, anisotropic tissue stress suppresses mismatches in the orientation of hexagonal arrays by endowing the cells with an inherent orientation (Fig. 9B). This accounts for a mechanism that allows for more rapid relaxation to the hexagonal cell pattern.

Fig. 9.

A physical mechanism of hexagonal cell packing. (A,B) Schematics of cell configurations under isotropic (A) and anisotropic (B) stress conditions. In the isotropic condition, nonhexagonal cells appeared between the mismatched top and bottom hexagonal arrays. A large force fluctuation was required to prevent the systems from being trapped in local energy minimum states.

Fig. 9.

A physical mechanism of hexagonal cell packing. (A,B) Schematics of cell configurations under isotropic (A) and anisotropic (B) stress conditions. In the isotropic condition, nonhexagonal cells appeared between the mismatched top and bottom hexagonal arrays. A large force fluctuation was required to prevent the systems from being trapped in local energy minimum states.

Flamingo (Fmi) was required for hexagonal cell packing when the global force balance was nearly stabilized

Finally, we studied the respective timings of actions of genetic and mechanical regulation in hexagonal cell packing. Among the PCP mutants examined, the flamingo (fmi) (Usui et al., 1999) mutant clone exhibited a relatively strong defect in hexagonalization (Classen et al., 2005) (see supplementary material Fig. S11A and its legend for evidence that that fmi RNAi did not significantly change mechanical parameters of the wing). Supplementary material Fig. S11B shows that the number of hexagonal cells increased in fmi RNAi wings, as was also observed in control wings, until 27.5 hours APF. After that, the increase and alignment of hexagonal cells was disrupted in fmi RNAi wings (supplementary material Fig. S11B,C, asterisks). Interestingly, anisotropic tissue stress began to weaken after 27 hours APF (Fig. 3A), and the alignment of local cellular forces to tissue stress was mostly completed around 27 hours APF (Fig. 5D-F). These results suggested that both of these forces and Fmi regulated hexagonal packing through the PD alignment of hexagonal cells, but that Fmi was indispensable for the process when the balance between the extrinsic and intrinsic forces was nearly stabilized (summarized in supplementary material Fig. S12A).

The present study aimed at unveiling the physical basis of tissue integrity via hexagonal cell packing. For this purpose, we utilized the Bayesian force-inference method to take a data-driven approach, in which the interplay between global stress in a tissue and local forces and shape changes of cells was analyzed. We determined that tissue stress remained highly anisotropic until early phase II of wing development when hexagonal cell packing occurred and that the mechanical anisotropy in a tissue promoted ordering of cell packing geometry (Fig. 9).

In the Drosophila wing, hexagonal cell packing is primarily achieved by cell rearrangement (Classen et al., 2005). Previous studies in other model systems, using subcellular and invasive force measurements, such as laser ablation of individual cell junctions, showed that the nonuniform localization of force-generating molecular machinery generates the angular bias in intrinsic cell junction tension, which triggers directional cellular rearrangements underlying tissue morphogenesis (supplementary material Fig. S12B, arrows 1 and 2) (e.g. Rauzi et al., 2008; Bosveld et al., 2012). By mapping forces inside the tissue, the present study highlighted a pivotal role of another physical ingredient, namely tissue stress, in multicellular pattern formation; the maximum stress direction at the tissue level provided directional information for assignment of the orientation of hexagonal cells, leading to the organization of the hexagonally packed, structurally stable wing (supplementary material Fig. S12B, arrows 2 and 3).

As anisotropic stress of the wing emerges from the hinge constriction, mechanical interactions between different tissues or body parts can serve as input for the mechanical processes of cells and molecules. Such passive mechanical reactions affect the active mechanical processes of the force-generating molecular machinery that triggers cell and tissue deformation (Butler et al., 2009; Fernandez-Gonzalez et al., 2009; Pouille et al., 2009; Blanchard and Adams, 2011; Hoffman et al., 2011). Our finding that myosin redistributed in order to resist the extrinsic force in the wing (supplementary material Fig. S12B, arrows 3 and 4) represented one of these feedback processes. Future studies will examine whether similar mechanisms function in cell packing in other epithelial tissues and will further clarify the physical and biological principles behind the active and passive mechanical processes that coordinate multicellular pattern formation and tissue development.

We thank Tadashi Uemura (Kyoto University) for his continuous support and encouragement on this work; Kunihiko Kaneko (University of Tokyo) and Atsushi Miyawaki (RIKEN BSI) for giving us the opportunity to conduct this study; and Yohanns Bellaïche, Philippe Marcq (Institut Curie), François Graner (University Paris VII), James Alan Hejna (Kyoto University) and Shigeo Hayashi (RIKEN CDB) for providing critical comments regarding the manuscript. We are also grateful to Yang Hong (University of Pittsburgh), Roger Karess (Institut Jacques Monod), Fumio Matsuzaki (RIKEN CDB), the Bloomington Stock Center and the Drosophila Genetic Resource Center for reagents; Reiko Takahashi (RIKEN BSI) and Yuri Tsukahara (Kyoto University) for technical assistance; Minako Izutsu (Kyoto University) for fly food; and the RIKEN BSI-Olympus Collaboration Center and the iCeMS Imaging Center for imaging equipment.

Funding

This work was supported by grants from the RIKEN Incentive Research Grant program [K.S.]; The Ministry of Education, Culture, Sports, Science and Technology (MEXT) [K.S. and S.I.]; and PRESTO Japan Science and Technology Agency (JST) [S.I.]. K.S. was a RIKEN special postdoctoral researcher.

Author contributions

K.S. and S.I. designed the project, performed data analysis and wrote the paper. K.S. conducted experiments and S.I. conducted numerical simulations.

Aegerter-Wilmsen
T.
,
Smith
A. C.
,
Christen
A. J.
,
Aegerter
C. M.
,
Hafen
E.
,
Basler
K.
(
2010
).
Exploring the effects of mechanical feedback on epithelial topology
.
Development
137
,
499
506
.
Aigouy
B.
,
Farhadifar
R.
,
Staple
D. B.
,
Sagner
A.
,
Röper
J. C.
,
Jülicher
F.
,
Eaton
S.
(
2010
).
Cell flow reorients the axis of planar polarity in the wing epithelium of Drosophila
.
Cell
142
,
773
786
.
Batchelor
G. K.
(
1970
).
The stress system in a suspension of force-free particles
.
J. Fluid Mech.
41
,
545
570
.
Blanchard
G. B.
,
Adams
R. J.
(
2011
).
Measuring the multi-scale integration of mechanical forces during morphogenesis
.
Curr. Opin. Genet. Dev.
21
,
653
663
.
Bonnet
I.
,
Marcq
P.
,
Bosveld
F.
,
Fetler
L.
,
Bellaïche
Y.
,
Graner
F.
(
2012
).
Mechanical state, material properties and continuous description of an epithelial tissue
.
J. R. Soc. Interface
9
,
2614
2623
.
Bosveld
F.
,
Bonnet
I.
,
Guirao
B.
,
Tlili
S.
,
Wang
Z.
,
Petitalot
A.
,
Marchand
R.
,
Bardet
P. L.
,
Marcq
P.
,
Graner
F.
, et al. 
. (
2012
).
Mechanical control of morphogenesis by Fat/Dachsous/Four-jointed planar cell polarity pathway
.
Science
336
,
724
727
.
Butler
L. C.
,
Blanchard
G. B.
,
Kabla
A. J.
,
Lawrence
N. J.
,
Welchman
D. P.
,
Mahadevan
L.
,
Adams
R. J.
,
Sanson
B.
(
2009
).
Cell shape changes indicate a role for extrinsic tensile forces in Drosophila germ-band extension
.
Nat. Cell Biol.
11
,
859
864
.
Chiou
K. K.
,
Hufnagel
L.
,
Shraiman
B. I.
(
2012
).
Mechanical stress inference for two dimensional cell arrays
.
PLOS Comput. Biol.
8
,
e1002512
.
Classen
A. K.
,
Anderson
K. I.
,
Marois
E.
,
Eaton
S.
(
2005
).
Hexagonal packing of Drosophila wing epithelial cells by the planar cell polarity pathway
.
Dev. Cell
9
,
805
817
.
Cooper
M. A.
,
Son
A. I.
,
Komlos
D.
,
Sun
Y.
,
Kleiman
N. J.
,
Zhou
R.
(
2008
).
Loss of ephrin-A5 function disrupts lens fiber cell packing and leads to cataract
.
Proc. Natl. Acad. Sci. USA
105
,
16620
16625
.
Dahmann
C.
,
Oates
A. C.
,
Brand
M.
(
2011
).
Boundary formation and maintenance in tissue development
.
Nat. Rev. Genet.
12
,
43
55
.
Eaton
S.
,
Jülicher
F.
(
2011
).
Cell flow and tissue polarity patterns
.
Curr. Opin. Genet. Dev.
21
,
747
752
. (
Ï
).
Farhadifar
R.
,
Röper
J. C.
,
Aigouy
B.
,
Eaton
S.
,
Jülicher
F.
(
2007
).
The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing
.
Curr. Biol.
17
,
2095
2104
.
Fernandez-Gonzalez
R.
,
Simoes
S. M.
,
Röper
J. C.
,
Eaton
S.
,
Zallen
J. A.
(
2009
).
Myosin II dynamics are regulated by tension in intercalating cells
.
Dev. Cell
17
,
736
743
.
Fisher
N. I.
(
1993
).
Statistical Analysis of Circular Data
.
Cambridge
:
Cambridge University Press
.
Fristrom
D.
,
Fristrom
J. W.
(
1993
).
The metamorphic development of the adult epidermis
. In
The Development of Drosophila Melanogaster
(ed.
Bate
M.
,
Arias
A. Martinez
), pp.
843
897
.
Cold Spring Harbor, NY
:
Cold Spring Harbor Laboratory Press
.
Graner
F.
,
Sawada
Y.
(
1993
).
Can surface adhesion drive cell rearrangement? Part II: a geometrical model
.
J. Theor. Biol.
164
,
477
506
.
Hales
T. C.
(
2001
).
The honeycomb conjecture
.
Discrete Comput. Geom.
25
,
1
22
.
Hoffman
B. D.
,
Grashoff
C.
,
Schwartz
M. A.
(
2011
).
Dynamic molecular processes mediate cellular mechanotransduction
.
Nature
475
,
316
323
.
Honda
H.
(
1983
).
Geometrical models for cells in tissues
.
Int. Rev. Cytol.
81
,
191
248
.
Honda
H.
,
Motosugi
N.
,
Nagai
T.
,
Tanemura
M.
,
Hiiragi
T.
(
2008
).
Computer simulation of emerging asymmetry in the mouse blastocyst
.
Development
135
,
1407
1414
.
Huang
J.
,
Zhou
W.
,
Dong
W.
,
Watson
A. M.
,
Hong
Y.
(
2009
).
From the Cover: Directed, efficient, and versatile modifications of the Drosophila genome by genomic engineering
.
Proc. Natl. Acad. Sci. USA
106
,
8284
8289
.
Hutson
M. S.
,
Tokutake
Y.
,
Chang
M. S.
,
Bloor
J. W.
,
Venakides
S.
,
Kiehart
D. P.
,
Edwards
G. S.
(
2003
).
Forces for morphogenesis investigated with laser microsurgery and quantitative modeling
.
Science
300
,
145
149
.
Ishihara
S.
,
Sugimura
K.
(
2012
).
Bayesian inference of force dynamics during morphogenesis
.
J. Theor. Biol.
313
,
201
211
.
Ishihara
S.
,
Sugimura
K.
,
Cox
S. J.
,
Bonnet
I.
,
Bellaïche
Y.
,
Graner
F.
(
2013
).
Comparative study of non-invasive force and stress inference methods in tissue
.
Eur. Phys. J. E
36
,
45
.
Käfer
J.
,
Hayashi
T.
,
Marée
A. F.
,
Carthew
R. W.
,
Graner
F.
(
2007
).
Cell adhesion and cortex contractility determine cell patterning in the Drosophila retina
.
Proc. Natl. Acad. Sci. USA
104
,
18549
18554
.
Kasza
K. E.
,
Zallen
J. A.
(
2011
).
Dynamics and regulation of contractile actin-myosin networks in morphogenesis
.
Curr. Opin. Cell Biol.
23
,
30
38
.
Kiehart
D. P.
,
Lutz
M. S.
,
Chan
D.
,
Ketchum
A. S.
,
Laymon
R. A.
,
Nguyen
B.
,
Goldstein
L. S.
(
1989
).
Identification of the gene for fly non-muscle myosin heavy chain: Drosophila myosin heavy chains are encoded by a gene family
.
EMBO J.
8
,
913
922
.
Knox
A. L.
,
Brown
N. H.
(
2002
).
Rap1 GTPase regulation of adherens junction positioning and cell adhesion
.
Science
295
,
1285
1288
.
Koto
A.
,
Kuranaga
E.
,
Miura
M.
(
2009
).
Temporal regulation of Drosophila IAP1 determines caspase functions in sensory organ development
.
J. Cell Biol.
187
,
219
231
.
Lecuit
T.
,
Lenne
P. F.
(
2007
).
Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis
.
Nat. Rev. Mol. Cell Biol.
8
,
633
644
.
Lecuit
T.
,
Lenne
P. F.
,
Munro
E.
(
2011
).
Force generation, transmission, and integration during cell and tissue morphogenesis
.
Annu. Rev. Cell Dev. Biol.
27
,
157
184
.
Ma
D.
,
Amonlirdviman
K.
,
Raffard
R. L.
,
Abate
A.
,
Tomlin
C. J.
,
Axelrod
J. D.
(
2008
).
Cell packing influences planar cell polarity signaling
.
Proc. Natl. Acad. Sci. USA
105
,
18800
18805
.
Mofrad
M. R. K.
,
Kamm
R. D.
(
2006
).
Cytoskeletal Mechanics
.
Cambridge
:
Cambridge University Press
.
Nienhaus
U.
,
Aegerter-Wilmsen
T.
,
Aegerter
C. M.
(
2009
).
Determination of mechanical stress distribution in Drosophila wing discs using photoelasticity
.
Mech. Dev.
126
,
942
949
.
O’Keefe
D. D.
,
Gonzalez-Niño
E.
,
Edgar
B. A.
,
Curtiss
J.
(
2012
).
Discontinuities in Rap1 activity determine epithelial cell morphology within the developing wing of Drosophila
.
Dev. Biol.
369
,
223
234
.
Oda
H.
,
Uemura
T.
,
Harada
Y.
,
Iwai
Y.
,
Takeichi
M.
(
1994
).
A Drosophila homolog of cadherin associated with armadillo and essential for embryonic cell-cell adhesion
.
Dev. Biol.
165
,
716
726
.
Ouchi
N. B.
,
Glazier
J. A.
,
Rieu
J. P.
,
Upadhyaya
A.
,
Sawada
Y.
(
2003
).
Improving the realism of the cellular Potts model in simulations of biological cells
.
Physica A
329
,
451
458
.
Paluch
E.
,
Heisenberg
C. P.
(
2009
).
Biology and physics of cell shape changes in development
.
Curr. Biol.
19
,
R790
R799
.
Pouille
P. A.
,
Ahmadi
P.
,
Brunet
A. C.
,
Farge
E.
(
2009
).
Mechanical signals trigger Myosin II redistribution and mesoderm invagination in Drosophila embryos
.
Sci. Signal.
2
,
ra16
.
Raphael
Y.
,
Altschuler
R. A.
(
2003
).
Structure and innervation of the cochlea
.
Brain Res. Bull.
60
,
397
422
.
Rauzi
M.
,
Verant
P.
,
Lecuit
T.
,
Lenne
P. F.
(
2008
).
Nature and anisotropy of cortical forces orienting Drosophila tissue morphogenesis
.
Nat. Cell Biol.
10
,
1401
1410
.
Royou
A.
,
Field
C.
,
Sisson
J. C.
,
Sullivan
W.
,
Karess
R.
(
2004
).
Reassessing the role and dynamics of nonmuscle myosin II during furrow formation in early Drosophila embryos
.
Mol. Biol. Cell
15
,
838
850
.
Sano
H.
,
Kunwar
P. S.
,
Renault
A. D.
,
Barbosa
V.
,
Clark
I. B.
,
Ishihara
S.
,
Sugimura
K.
,
Lehmann
R.
(
2012
).
The Drosophila actin regulator ENABLED regulates cell shape and orientation during gonad morphogenesis
.
PLoS ONE
7
,
e52649
.
Shimada
Y.
,
Yonemura
S.
,
Ohkura
H.
,
Strutt
D.
,
Uemura
T.
(
2006
).
Polarized transport of Frizzled along the planar microtubule arrays in Drosophila wing epithelium
.
Dev. Cell
10
,
209
222
.
Strutt
H.
,
Strutt
D.
(
2009
).
Asymmetric localisation of planar polarity proteins: Mechanisms and consequences
.
Semin. Cell Dev. Biol.
20
,
957
963
.
Sugimura
K.
,
Bellaïche
Y.
,
Graner
F.
,
Marcq
P.
,
Ishihara
S.
(
2013
).
Robustness of force and stress inference in an epithelial tissue
.
Conf. Proc. IEEE Eng. Med. Biol. Soc.
(
in press
)
Thompson
D. W.
(
1917
).
On Growth and Form
.
Cambridge
:
Cambridge University Press
.
Usui
T.
,
Shima
Y.
,
Shimada
Y.
,
Hirano
S.
,
Burgess
R. W.
,
Schwarz
T. L.
,
Takeichi
M.
,
Uemura
T.
(
1999
).
Flamingo, a seven-pass transmembrane cadherin, regulates planar cell polarity under the control of Frizzled
.
Cell
98
,
585
595
.
Wootton
R. J.
(
1992
).
Functional morphology of insect wings
.
Annu. Rev. Entomol.
37
,
113
140
.

Competing interests statement

The authors declare no competing financial interests.

Supplementary information