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

How tissue shape emerges from the collective mechanical properties and behavior of individual cells is not understood. We combine experiment and theory to study this problem in the developing wing epithelium of Drosophila. At pupal stages, the wing-hinge contraction contributes to anisotropic tissue flows that reshape the wing blade. Here, we quantitatively account for this wing-blade shape change on the basis of cell divisions, cell rearrangements and cell shape changes. We show that cells both generate and respond to epithelial stresses during this process, and that the nature of this interplay specifies the pattern of junctional network remodeling that changes wing shape. We show that patterned constraints exerted on the tissue by the extracellular matrix are key to force the tissue into the right shape. We present a continuum mechanical model that quantitatively describes the relationship between epithelial stresses and cell dynamics, and how their interplay reshapes the wing. DOI: http://dx.doi.org/10.7554/eLife.07090.001


Introduction
The dynamic choreography of tissue shape changes that occur during development dramatically illustrates the fact that morphogenesis depends on organized cellular force generation. The mechanisms that control the orientation and patterning of these active processes and the corresponding tissue stresses are beginning to be explored in a variety of developmental systems, for review (Lecuit and Lenne, 2007;Keller, 2012;Heisenberg and Bellaiche, 2013). However, a complete understanding of the mechanical basis of morphogenesis will require not only a description of cell autonomously generated forces, but also quantitative insights into how cells respond to tissue stresses. Cells exert forces on extracellular matrices, but also on each other-this is particularly true of epithelial cells, which are tightly connected by specialized adhesive junctions. Thus, stresses generated by one epithelial cell can be transmitted to others throughout the tissue. In vitro experiments have shown that tissues respond to stress elastically over short time scales but that they can plastically remodel when subjected to stress over longer times (Beloussov et al., 2000;Harris et al., 2012). This can occur as a result of cell shape changes, cell rearrangements or both, and appears to vary with the cell types examined. Furthermore, experiments with cultured epithelial cells suggest that tissue compression can limit cell proliferation in vitro (Puliafito et al., 2012;Streichan et al., 2014). How these cellular responses might influence tissue size and shape in vivo is not clear. Nevertheless, these in vitro observations suggest that a complete and quantitative understanding of tissue morphogenesis will require new insights into tissue viscoelasticity in vivo and the cellular mechanisms that give rise to it.
Drosophila pupal wing morphogenesis is an ideal system in which to study the interplay of cellular force generation and tissue material properties in vivo. During pupal stages, anisotropic stresses along the proximal-distal (PD) axis of the wing blade epithelium help guide anisotropic tissue flows that reshape the blade-elongating it in the PD axis and narrowing it in the anterior-posterior (AP) axis, for review (Eaton and Julicher, 2011). The mechanisms that produce PD-oriented stresses in the wing blade are not fully understood. They are generated in part by contraction of cells in the wing hinge, which connects to the wing blade on its proximal side. However, we do not understand the origin of counterforces that restrain movement of the wing blade at the margin.
Analyzing cells in a subregion of the wing blade showed that tissue flows are associated with cell shape changes, cell divisions and cell rearrangements that are oriented along the PD axis (Aigouy et al., 2010). To quantitatively understand the cellular basis of this tissue shape change, we must determine the global patterns of these cellular events throughout the wing blade. Furthermore, while hinge contraction contributes to PD tissue stresses in the blade, cells in the wing blade might also contribute autonomously to tissue flows and stresses. Thus, to understand the mechanical basis of pupal wing morphogenesis, we must understand the emergence of PD-oriented stresses in the wing blade, and distinguish stresses autonomously generated by wing epithelial cells from the response of epithelial cells to these stresses.
Here, we combine several quantitative methods to investigate how cell flows and global tissue shape changes emerge from the collective behavior and mechanical properties of many wing epithelial cells. We develop image analysis methods to track the majority of cells in the wing throughout morphogenesis, and analyze cell shapes and rearrangements of the junctional network. Furthermore, we develop theoretical methods to quantify the cellular contributions to tissue shear and area homeostasis in the wing blade. eLife digest The individual cells in a developing animal embryo organize themselves into tissues with specific and reproducible shapes, which requires the cells to communicate with one another. Cells in tissues exert forces on their neighbors, and respond to being pushed and pulled by the cells around them.
In the fruit fly Drosophila melanogaster, each wing consists mainly of a framework of proteins and other molecules that is built by epithelial cells. These epithelial cells divide and grow during the life of a fly larva, and then reorganize themselves into the shape of the wing after it forms into a pupa. During this reshaping, epithelial cells in some regions of the wing experience powerful contractions. Previous work had suggested that and these forces produced tension in the rest of the wing to pull it into its final elongated shape. But it wasn't clear what exactly these contractions were pulling against to produce the tension. Nor was it understood exactly how wing epithelial cells responded to tension to reorganize themselves into a different wing shape. Now, Etournay, Popovic, Merkel, Nandi et al. have analyzed the forces acting across the entire wing blade and how these forces shape the wing. All cell divisions, cell neighbor exchanges and changes in cell shape in the developing wing blade were tracked under a microscope; this revealed how each one of them contributed to the change in wing shape. Further experiments revealed that localized contractile forces produce tension in the wing because it is connected around its edge to surrounding structures via an extracellular protein called Dumpy. Releasing these contacts, by severing them with a laser or by mutating Dumpy, caused the wing to develop into abnormal shapes, showing that the tension in the wing blade has an important role in determining wing shape. Furthermore, by tracking cells in wings that had been severed by a laser, or mutated for Dumpy Etournay, Popovic, Merkel, Nandi et al. could figure out exactly which cellular processes were guided by epithelial tension.
Etournay Popovic, Merkel, Nandi et al. also present a theoretical model that describes how the interplay between active force generation and the response of cells to the resulting tension shapes the wings of fruit flies. They propose that epithelial tension provides a mechanism through which cells can communicate with each other to ensure that together the combined behavior of these cells generates reproducible shapes. Further studies are required to analyze how active force generation is patterned and cells sense and respond to external forces during development.
We show that localized apical extracellular matrix connections to the cuticle at the wing margin provide the counterforce to hinge contraction, and are required for the development of normal stresses in the wing blade. These stresses are essential to reshape the pupal wing while maintaining wing area homeostasis. We distinguish autonomously controlled from stress-driven cellular events, and present a continuum mechanical model that quantitatively explains wing shape changes on the basis of the relationship between tissue stress, cell elongation and cell rearrangements.

Results
Dumpy-dependent physical constraints at the margin maintain epithelial tension in the wing The emergence of two-dimensional stresses in the plane of the wing blade suggests that there are physical constraints on the movement of wing epithelial cells near the margin. We wondered whether there might be a matrix connecting the wing epithelium to the overlying pupal cuticle in this region. To investigate this, we used a laser to destroy the region between the margin of the E-Cadherin:GFP expressing wing epithelium and the cuticle after the two had separated as a consequence of molting. Although this treatment does not apparently damage either the wing or the cuticle, it causes the wing epithelium to rapidly retract away from the cuticle within seconds ( Figure 1A-B′′, Video 1). Laser ablation causes epithelial retraction when performed at any region along the wing blade margin-anteriorly, posteriorly or distally. During tissue flows, the now disconnected margin moves even further away from the cuticle, producing abnormal wing shapes ( Figure 1C-F). This shows that the wing is physically restrained by apical extracellular matrix connections to the overlying cuticle, and that these connections are required to shape the wing during tissue flows.
We wondered whether the large apical extracellular matrix protein Dumpy might contribute to these connections. Dumpy is a 2.5 MDa protein that is predicted to form filaments at least 1 μm long (Wilkin et al., 2000). It forms an elastic matrix in the embryonic tracheal lumen, and provides mechanical resilience of tendon cell attachments to the overlying cuticle (Dong et al., 2014). While dumpy null mutations are lethal, some hypomorphs produce wings that are short and misshapen-a defect that arises during pupal development (Waddington, 1939(Waddington, , 1940. To ask whether shape defects in dumpy wings might arise during pupal tissue flows, we imaged dumpy ov1 pupal wings that expressed E-Cadherin:GFP. The shape of dumpy ov1 wings is normal at 16 hr after puparium formation (APF), before molting occurs (Figure 2A,B). Shortly afterwards, when hinge contraction begins, the shape of the dumpy ov1 mutant wing blade begins to differ from wild type (WT). The wing blade epithelium retracts abnormally far from the distal cuticle and fails to elongate in the PD axis. By the time tissue flows have ended, the characteristic abnormal shape of the dumpy ov1 wing is apparent (Video 2 and Figure 2A-B′′).
To examine Dumpy distribution, we imaged wings from flies harboring a protein trap construct that expresses YFP:Dumpy from the endogenous chromosomal locus. YFP:Dumpy is present on the apical surface of epithelial cells throughout the wing, and within the overlying cuticle (Figure 2-figure supplement 1A, Video 3). Interestingly, Dumpy is also present in a fibrous-appearing matrix that connects the wing to the overlying cuticle in specific places. This matrix lies between the cuticle and the margin of the wing ( Figure 2C To investigate the extent to which wing margin constraints had been relieved by the dumpy ov1 mutation, we performed laser-severing experiments in the dumpy ov1 mutant background. Cutting a dumpy ov1 wing between the cuticle and the distal wing blade revealed almost undetectable retraction, suggesting that distal attachments of the wing blade to the cuticle are severely compromised. However retraction was still observed when severing was performed between the cuticle and the anterior or posterior margins (Figure 2-figure supplement 2). Furthermore, severing the matrix around the entire margin causes dumpy ov1 mutant wings to develop even more dramatic wing shape abnormalities during tissue flows (compare Figure 2-figure supplement 2G to Figure 2A′′). This suggests that apical matrix connections to the cuticle are not completely abrogated in dumpy ov1 wings.
To ask how the dumpy ov1 mutation influenced PD tension in the wing blade, we performed circular laser cuts covering about 5-10 cells in different regions of WT and dumpy ov1 wings ( Figure 2I,I′, Figure 2-figure supplement 3). We observed a recoil of the ablated region, indicating that the blade epithelium is under tension. From the recoil, we can compare both the isotropic and the anisotropic components of epithelial stress in WT and dumpy ov1 mutant wings ( Figure 2I,I′ and Figure 2-figure supplement 3). These stress patterns differ between WT and dumpy ov1 wings. The orientation of anisotropic stress in dumpy ov1 is somewhat splayed and not as well aligned with the PD axis. Furthermore, anisotropic tension in dumpy ov1 wings tends to be reduced in the central region and increased anteriorly and posteriorly.
Overall, Dumpy-dependent elastic connections are key to the emergence of the stress pattern during morphogenesis. This suggests that these stresses play an important role in guiding tissue flows. Dashed doublesided arrows depict the proximal-distal (PD) and anterior-posterior (AP) axes. The PD axis is defined by a regression line passing through selected sensory organs (red dots) that are easily identifiable in Ecad::GFP expressing wings. The x axis is defined to correspond to the PD axis pointing distally, and the y axis is defined to correspond to the AP axis pointing anteriorly. L2-L5 indicate longitudinal veins 2-5. Brown dashed line outlines the cuticular sac surrounding the wing epithelium. Scale bar 20 μm. (B, B′) Show the distal end of a wild-type (WT) Ecad::GFP-expressing wing at 24 hAPF (greyscale in B, B′) and the same wing 3.5 min after laser ablation in the space between wing margin and cuticle (magenta in B′). The blue dashed line indicates the site of laser ablation. (B′′) Shows wing margin displacement measured with respect to the cuticle (brown dashed line in B′) along the white dotted line in (B′). Experimental points (magenta) were interpolated by a polynomial (blue line). (C-F) Show 32 hAPF wings that were unperturbed (C) or subjected to laser ablation at 22 hAPF (D-F). Ablation of the connections between the wing margin and the cuticle were performed in different regions, indicated by blue dashed lines in (D-F), and lead to altered wing shapes at 32 hAPF compared to the unperturbed control (C). Scale bar 100 μm. DOI: 10.7554/eLife.07090.003

Quantifying wing morphogenesis
What are the cellular events that shape the wing blade during tissue flows? To quantitatively address this question, we developed methods to quantify cell shape changes, cell divisions, cell rearrangements and cell extrusions during wing morphogenesis. We imaged three E-Cadherin: GFP-expressing wings at cellular resolution every 5 min between 16 and 32 hr APF. We then extracted and projected the planes containing the apical adherens junctions, automatically detected cell edges, and tracked each cell in the wing over the course of the videos ( Figure 3A,A′ and Videos 4, 5). We designed a relational database (DB) to store information pertaining to all cells in a given video ('Materials and methods', Longterms time-lapse imaging and Data handling and image processing). Querying these DBs provides information about individual cellular properties such as shape, area, and associated cell boundaries. It also provides information about neighbor and lineage relationships, identifying neighbor exchanges (T1 transitions), cell divisions and cell extrusions (T2 transitions).

Cellular contributions to wing area changes
We showed previously that the area of the wing blade remains fairly constant during pupal morphogenesis, despite the fact that cells are dividing (Aigouy et al., 2010). To understand how the wing blade maintains area homeostasis, we quantified wing blade area, cell divisions, cell area changes and cell extrusions between 16 and 32 hAPF in three different WT wings ( Figure 3A-E). We analyzed a large region of the wing blade in which all cells could be tracked during the whole video (green in Figure 3A,C).
The tissue area expansion rate is v = 1 A dA dt , where A is the tissue area and d/dt denotes the total time derivative. The tissue expansion rate can be decomposed into its cellular contributions as where a is average cell area, k d and k e are cell division and extrusion rates, respectively. The cumulative area expansion rate of the whole wing blade is ln(A/A 0 ), where A 0 is the tissue area when recording starts (∼16 hAPF), can be obtained by integrating the area expansion rate v over time. Figure 3D shows both the tissue area expansion rate (dark blue) and the contributions to this expansion rate from cell area changes, cell divisions and extrusions. These are averages over three WT wings. Figure 3E shows the corresponding cumulative quantities. The dynamics of wing area changes in the 3 WT wing blades are extremely similar-after contracting slightly during the first half of morphogenesis, blade area gradually returns to very close to its original value (dark blue lines in Figure 3D,E). This almost constant area reflects a balance between cell divisions (orange) on the one hand, and cell area changes (green) and extrusions (light blue) on the other. Interestingly, blade area in three analyzed wings is more reproducible than would be expected from the variation in each cellular contribution, if they were independent of each other (see shaded regions depicting standard deviations in Figure 3D,E). To quantify this observation we compared the variance of overall relative area change with the sum of the variances of the cumulative cellular contributions. We find that sum of variances is about 20 times larger than variance of the sum. This shows that cellular contributions are not independent and that normal variations in the rate of cell division can be compensated by changes in cell area and/or extrusion to maintain wing blade area (  To ask whether connections to the cuticle were required to maintain wing blade area, we examined blade area changes and the underlying cellular contributions in dumpy ov1 mutant wings ( Figure 3F,G). As a complementary approach, we quantified the cellular contributions to area changes in lasersevered wings. We severed the wing either between the hinge and the blade, or at the very distal tip before hinge contraction occurred (Figure 3-figure supplement 2). We also severed connections between the wing margin and the cuticle ( Figure 1D,F) at about 22 hAPF. In contrast to unperturbed WT wings, total wing area decreases dramatically when connections at the margin are weakened by dumpy ov1 mutation or by laser severing ( Figure 3G and Figure 3-figure supplement 2, see dark blue curve). Thus, connections to the cuticle are required to maintain wing blade area during morphogenesis. These connections provide mechanical linkages that permit the buildup of tensile stresses while maintaining wing blade area.
How does epithelial stress influence the cellular events contributing to area homeostasis? To answer this question, we first compared cellular contributions to area change during morphogenesis of WT and dumpy ov1 mutant wings. Wing blade area decrease in dumpy ov1 mutant wings is not a consequence of fewer cell divisions-cells actually divide more than in WT (Figure 3-figure supplement 2D-F,I). Cells in dumpy ov1 mutant wings have a similar maximum division rate but divide over a longer period of time, resulting in more cells at the end of morphogenesis ( Figure 3G and Figure 3-figure supplement 2D-F, yellow curve). The reduced wing blade area in dumpy ov1 is quantitatively explained by reductions in cell area and by cell extrusions, which more than compensate the increased proliferation. Thus, reduced epithelial stresses in dumpy ov1 wings perturb the balance between cell divisions, area changes and extrusions seen in WT. All laser-severing perturbations decrease the final wing area, similar to dumpy ov1 mutant wings (Figure 3-figure supplement 2). In these wings, the analysis of cellular contributions to wing area changes is complicated by the delay between cutting and the initiation of time-lapse imaging (about 45 min). During this intervening time, the wing responds acutely to reduced tension, and both wing and cell area decrease to values below those expected for WT wings of the same stage (Figure 3-figure supplement 2D-F,H). While we can estimate changes to cell area during this time, we cannot know the rates of cell division and extrusion. Nevertheless, several interesting conclusions can be drawn by analyzing final cell Video 2. Synchronized time-lapses of wild-type (WT) and dumpy ov1 wings. The synchronization is based on the time when histoblast nests merge at ∼26.5 hAPF. DOI: 10.7554/eLife.07090.009 area, and the rates of division, area change and extrusion after time-lapse imaging begins. Wings that have been severed before hinge contraction (whether at the hinge-blade interface or at the distal tip) behave similarly to dumpy ov1 mutant wings. After an initial delay, the rate of cell division increases and cells divide more than in unperturbed wings. However, cell extrusions and decreasing cell area more than compensate for increased cell division to produce smaller wings. When cuticle connections are severed later at 22 hAPF, most cell divisions have already occurred, and this treatment does not increase proliferation. In this case, the number of cell extrusions increases, and the final cell area is smaller than that of unperturbed wings. Taken together, analyzing dumpy ov1 and laser-severed wings shows that epithelial stresses are required to balance cell divisions with cell extrusions and cell area changes to maintain area homeostasis during morphogenesis. This is consistent with observations in the thorax, where overcrowding drives delamination (Marinari et al., 2012).

A method to quantify cellular contributions to wing blade deformation
In the preceding section, we discussed how cellular processes contributed to wing blade area changes. These area changes correspond to the isotropic component of a tensor characterizing the tissue strain ( Figure 4A). Now, we discuss shape changes of the wing blade, which correspond to the anisotropic part of this tensor and characterize the process of elongation along an axis (i.e., pure shear). The rate of change of pure shear is described by the pure shear rate tensor v ∼ ( Figure 4B) with v ∼ xx characterizing the rate of elongation along the PD axis ( Figure 1A). Note that pure shear, that is, convergence-extension flow, is different from so-called simple shear, which results from a superposition of pure shear and a rotation ( Figure 4C). In the following, we use the term shear to refer to pure shear. We now discuss how tissue shear can be decomposed into contributions from cell shape changes and topological changes. These include cell divisions, cell neighbor exchanges (T1 transitions), and cell extrusions (T2 transitions) ( Figure 4D-G).
To understand the cellular contributions to the overall shear of the wing blade, we developed a method to distinguish and quantify shear caused by cell shape changes and shear caused by topological changes. In a piece of tissue where no topological changes occur within the cellular network, the deformation of the individual cells defines the deformation of the whole piece of tissue. However, when topological changes occur, deformation of individual cells no longer completely accounts for the overall shear ( Figure 4D-H). The triangle method we outline below represents an exact geometrical formalism to decompose large-scale deformations of the wing blade into contributions by cell deformation and by each kind of topological change.
First, we tile the cellular network with triangles as follows. Each vertex of the cellular network that touches three cells (red dot in Figure 4J) gives rise to a single triangle (red), whose corners are defined by the centers of the three cells (green dots). Vertices that touch more than three cells are treated as described in Appendix 1, 'Triangulation procedure'. The resulting set of triangles tiles the cellular network without gaps or overlaps ( Figure 4K). We choose a tiling into triangles, because the deformation of a single triangle between two frames of the video can be uniquely characterized by a single 2 × 2 tensor describing a linear transformation (see Appendix 1, 'Triangle deformation'). Note that such a characterization by a 2 × 2 tensor is in general not possible for polygons with more than three sides. For each triangle and time point, this tensor describes relative area changes, rotation, and shear of the triangle. The average shear rate of all triangles in the tissue corresponds to the overall tissue shear rate. To connect the tissue shear rate to cell elongation changes, we define a nematic tensor Q characterizing the state of triangle elongation (see Appendix 1, 'Triangle elongation'). Then, the change of triangle elongation corresponds exactly to triangle shear. Cell elongation is obtained as the average of the elongation tensors Q of those triangles that belong to a given cell. Hence, in the absence of topological changes, we find an exact relation between cell elongation change and overall tissue shear (Appendix 1, 'Large-scale shear in the absence of topological transitions').   If topological changes occur, then average cell elongation change is not equal to tissue shear. To include the effect of topological transitions, we write (Appendix 1, 'Contributions to shear by topological changes'): where v ∼ is the overall tissue shear rate tensor, DQ/Dt is a corotational rate of change in average triangle elongation, and R is the shear rate tensor due to topological changes. Contributions to R include T1 transitions (T ), cell divisions (C ), T2 transitions (E ), as well as correlated cell shape changes and cell rotations (D): R = T + C + E + D (see Figure 4I). How can we define the contributions by topological changes (T, C, and E ) to tissue shear? During a topological transition, the triangulation changes and thus the average triangle elongation changes ( Figure 4L-O). However, at the moment the topological change occurs there is no tissue shear. Therefore, tissue shear and triangle elongation are no longer the same. This can be compensated by introducing a contribution to tissue shear by topological transitions. This contribution corresponds to the negative change in the average triangle elongation caused by the change in the triangulation ( Figure 4M and Appendix 1, 'Intermediate network states' and Appendix 1, 'Contributions to shear by topological changes').
In the definition of R in Equation 2, we have also introduced the contribution D to tissue shear, which accounts for collective cellular events that combine to increase average triangle elongation in the absence of tissue shear and topological transitions. This occurs when several triangles have fluctuating shapes, such that the instantaneous elongation and the rotation rate or area expansion rate of triangles are correlated. Note that this effect does not occur when several triangles undergo equal deformations and rotations. One example of a cellular network deformation that produces the contribution D to tissue shear is shown in Figure 4H. Here, cells in neighboring rows slide relative to each other in alternating directions, such that no net pure shear occurs. However, there are alternating rows of simple shear and a net change in triangle elongation. We introduce the contribution D in the definition of R in order to compensate for the increase in the average triangle elongation DQ/Dt stemming from such correlations, should they exist in the wing blade.

Patterns of cellular contributions to tissue shear
We first used the triangle method to calculate the patterns of tissue shear and cellular contributions to this tissue shear in WT ( Figure 5, Video 6, and Appendix 1, 'Spatially resolved shear patterns') and dumpy ov1 mutant wings ( Figure 5-figure supplement 1 and Video 7). To visualize these patterns we averaged all quantities within squares of about 26 × 26 μm 2 . Figure 5 shows shear patterns in WT at early, intermediate, and late time points during pupal wing morphogenesis. Shear is indicated by a line whose orientation represents the shear axis and whose magnitude corresponds to the shear rate. Tissue shear in the WT wing blade is oriented in a fan-shaped pattern with a strong PD component ( Figure 5A-A′′). At about 21 hAPF, the shear pattern develops a sharp reorientation between veins L3 and L4, where shear is oriented along the AP axis. This region corresponds to a stripe of Dumpy-containing matrix that attaches the blade What are the patterns of cellular contributions to the tissue shear? These patterns reveal a surprising complexity that changes with time. Shear caused by cell elongation and cell rearrangements (T1 transitions) display significant contributions that are antagonistic. Unexpectedly, T1 transitions cause shear along the AP axis early in the process ( Figure 5B-B′′). This shear is opposed by increasing cell elongation along the PD axis ( Figure 5C-C′′). At intermediate and late stages, T1 transitions shift their average orientation to cause PD shear-at the same time, cells reduce their elongation along the PD axis, causing AP shear. The contribution of cell extrusions to tissue shear is negligible (not shown) and cell divisions result in significant shear only during the early stages of wing morphogenesis ( Figure 5D-D′′). Finally, plotting the contribution D to shear due to correlation effects reveals that these effects do exist in the wing ( Figure 5E-E′′). Thus, patterns of cell elongation, cell rearrangement, cell division and correlation effects make dynamically changing contributions to tissue shear that are sometimes antagonistic.
To investigate the cause of altered wing shape in dumpy ov1 , we performed a similar analysis. The patterns of cellular contributions to tissue shear in a dumpy ov1 mutant wing display subtle abnormalities ( Figure 5-figure supplement 1). However, a more quantitative analysis is required to understand the origin of the altered dumpy ov1 mutant wing shape.

Total cellular contributions to tissue shear
To better understand the quantitative relationships between the cellular processes contributing to tissue shear, we studied spatially averaged shear in the wing blade projected on the PD axis. We quantified average tissue shear, and shear caused by each cellular contribution over time in three different WT videos. These averages were taken over the tracked region shown in Figure 3A and the resulting quantities were further averaged over the three videos. These averages and the standard deviations between the videos are shown in Figure 6A,B. Positive values indicate shear along the PD axis and negative values indicate shear orthogonal to the PD axis (i.e., shear in the AP axis). Adding the contributions of shear caused by cell divisions, cell rearrangements, cell shape changes, and correlation effects (light-pink line in Figure 6-figure supplement 1A) reproduces the independently calculated total shear curve (blue line in Figure 6A and Figure 6-figure supplement 1A). Small differences between these two curves (about 3%) stem from small inaccuracies (see Appendix 1, 'Decomposition of the large-scale tissue shear rate'). Thus, we can decompose tissue shear into its individual cellular contributions.
Does the cumulative tissue shear we calculate account for the shape change of the wing blade? To verify this, we characterize the blade shape by a nematic determined by the outline of the tracked region (see Appendix 1, 'Characterization of wing blade anisotropy'). The change of this nematic with respect to its initial value over time is shown together with the cumulative tissue shear obtained by the triangle method ( Figure 6-figure supplement 1B). These quantities agree well, indicating that the shear projected on the PD axis accounts for the main features of tissue shape changes. Thus, we can now use the shear decomposition method to discuss how different cellular events contribute to shape change of the wing over time.
On average, the wing shears smoothly along its PD axis between 17 and 32 hAPF as the hinge contracts. In contrast, the cellular processes that combine to produce tissue shear change over time. During the first 6 hr of our videos, shear caused by cell elongation in the PD axis (green curves in Figure 6A,B) is even larger than PD tissue shear (blue curves). This shows that cell elongation increases more than the tissue elongates suggesting that active cellular processes also contribute to PD cell elongation. Subsequently, starting at about 22 hAPF, cell elongation decreases although the blade continues to elongate. These discrepancies between cell shape changes and tissue shape changes require topological changes in the cell network. To more clearly discuss these events, we define two distinct phases of wing morphogenesis (phases I and II) that are separated by the peak of cell elongation occurring at about 22.5 hAPF. Quantifying shear caused by T1 transitions and by cell elongation reveals that they change dynamically throughout morphogenesis with a striking reciprocal relationship (green and red curves in Figure 6A,B). This reciprocal relationship accounts to a large extent for the discrepancies between cell elongation and tissue shear. It further suggests that the active contribution to cell elongation (i.e., the amount of cell elongation that exceeds tissue shear) may be linked to AP-oriented T1 transitions-the orientation of these T1s, which work against the observed tissue shear, suggests that they are autonomously controlled. Active AP-oriented T1 transitions could produce PD cell elongation if mechanical constraints prevent the wing from shearing. In principle, it is also possible that active PD-oriented cell shape changes could produce AP-oriented T1 transitions under the same constraints. Cell divisions also contribute to PD shear The shear rate and shear rate contribution tensors were locally averaged within 26 × 26 μm 2 square elements (25-50 cells) of a fixed grid. A 45 min time window was used to smooth the shear values within each grid element. The resulting nematic tensors are represented by line segments whose length and direction correspond to the norm and orientation of the tensor, respectively. Scale bar: 100 μm. DOI: 10.7554/eLife.07090.017 The following figure supplement is available for figure 5: in the wing blade. Although cell divisions initially cause a small amount of AP shear, their direction changes during phase I such that their net contribution shears the wing in the PD axis. In addition, correlation effects produce significant shear in the AP direction and contribute most strongly at the time that T1 transitions are changing their orientation.
In summary, the continuous large-scale deformation of the wing blade emerges from complex patterns of cell dynamics on small scales. During phase I, cells undergo AP-oriented T1 transitions while elongating in the PD axis. Cell divisions during phase I contribute shear along the PD axis. During phase II, the orientation of T1 transitions shifts to the PD axis and cells relax their shape. Correlation effects contribute AP shear, and peak roughly at the time that T1's shift their orientation.
To ask whether the cellular contributions to tissue shear were independent of each other, we compared the sum of the variances of the final cellular contributions to tissue shear with the variance of final tissue shear itself. The ratio of these values is about 25, indicating that the cellular contributions to tissue shear-like the contributions to area change-are not independent of each other. Thus, the overall tissue shear is more reproducible than would be expected from the variations of the cellular contributions.

Alternating directions of shear and rotation in neighboring cell rows contribute to tissue shear
What exactly are cells doing that results in correlation effects? We found that shear due to correlation effects was mainly generated by correlations between local elongation and rotations ( Figure 6C-I).
To investigate this further, we determined the magnitude and orientation of shear and the rotation rate associated with each triangle for each frame of the video. We observed that triangles rotated and sheared in striking spatial patterns that rapidly fluctuate in time ( Figure 6C,F,G). These patterns correspond to rows of correlated shear and rotation that are distributed throughout the wing blade. To characterize these correlated patterns, we calculated the spatial power spectrum of the local tissue rotation rate ( Figure 6D, Appendix 1, 'Power spectrum of local tissue rotation'). This revealed that shear and rotation are correlated in regions corresponding to PD-oriented rows that were about 3-7 triangles long. These rows consist of triangles that all rotate and shear in the same direction. The rows are interspersed with other rows of similar length with mirrored patterns of shear and rotation (note blue and red rows in Figure 6F,G). Such a pattern of rotation and pure shear is characteristic of neighboring rows of triangles and cells undergoing simple shear in alternating directions ( Figure 6H). This would occur if PD-oriented rows of cells slide with respect to each other. As discussed above, such rearrangements can indeed contribute to correlation terms (see cartoon in Figure 4H).
If rows of cells slide past each other, cells typically engage in T1 transitions. Since the peak of AP correlation effects coincided with a shift in the net orientation of T1 transitions (i.e., when the red curve in Figure 6A crosses 0), we wondered whether correlation effects could be associated with T1 transitions at this time (see also Figure 6-figure supplement 1C). We therefore examined whether correlation effects were associated with a particular type of topological change. Indeed, correlation effects are mostly accounted for by those cells that are about to undergo a T1 transition within the next 9 frames, although they cover only a small fraction of the total area ( Figure 6I, Appendix 1, 'Role of T1 transitions in the correlation-induced shear'). For rows of cells to slide past each other, cells would have to undergo a peculiar type of T1 Video 6. Dynamic patterns of tissue shear and of its cellular contributions in a WT wing. Line segments are nematic representations of shear. DOI: 10.7554/eLife.07090.019 Video 7. Dynamic patterns of tissue shear and of its cellular contributions in a dumpy ov1 mutant. Line segments are nematic representations of shear. DOI: 10.7554/eLife.07090.020 Figure 6. Total cellular contributions to tissue shear throughout the WT wing blade. (A) Shows the tissue shear rate (blue) over time, and shear rates contributed by cell rearrangements (red), cell shape changes (green), cell divisions (orange), and correlation effects (magenta), averaged throughout the wing blade. These averages were taken over the tracked region shown in Figure 3A by averaging nematic tensors throughout the wing blade. The resulting quantities were further projected onto the PD axis and averaged over the three WT videos. Ribbons indicate the standard deviation between wings. The sign of the shear rate defines its orientation (>0 is PD-oriented and <0 is AP-oriented). (B) Shows the accumulated tissue shear over time throughout the blade, and the accumulated contributions of each cellular process (color code as in A). (C) Pattern of local tissue rotation rate at 21 hr after puparium formation (APF). The local tissue rotation rate ω m is plotted separately for each triangle m. Red circles correspond to a counter-clockwise rotation and blue circles correspond to a clockwise rotation. The area of each circle scales with the absolute rotation rate. (D) The spatial power spectrum of the local tissue rotation rate corresponding to the pattern in panel C (see Appendix 1, 'Power spectrum of local tissue rotation'). The power spectrum is a function of a wave vector q = (q x , q y ), which is measured in units of a typical cell diameter d 0 = 4 μm. The two peaks in the power spectrum at q peak ≈ (0, ±0.3d 0 /2π) correspond to the existence of horizontal bands of alternating tissue rotation that are separated by about 1.5 cell diameters (compare panels F, G). (E) Correlation effects contributing to shear along the PD axis, D xx (magenta curve). D xx can be decomposed into an area expansion part D e xx (green curve), which corresponds to a correlation between the local area expansion rate v m and local triangle elongation Q m xx : D e xx = −ðAEv m Q m xx ae − AEv m aeAEQ m xx aeÞ, and into a rotational part D r xx (blue curve), which corresponds to a correlation between the local tissue rotation rate ω m and local triangle elongation: D r xx ≃2ðAEω m Q m xy ae − AEω m aeAEQ m xy aeÞ (see Appendix 1, 'Large-scale shear in the absence of topological transitions'). The rotational part dominates the shear by correlation effects. (F) Enlargement of the rotation pattern in panel C with an additional indication of the pattern of the local shear rate tensor by green bars. Length and orientation of a bar correspond to magnitude and axis of the local shear rate, respectively. The axis of local shear is correlated with the sign of local rotation (indicated by red and blue circles). (G) The same region of the wing in the subsequent frame (about 5 min later). Three corresponding triangles in panels F and G are colored in cyan, yellow and orange, respectively. The patterns of local shear and rotation change on time scales of minutes. (H) A correlation of local rotation and local shear within bands as shown in panels F, G corresponds to bands of alternating simple shear. (I) Contribution to the shear due to correlation effects of the group of triangles that are going to disappear due to a T1 transition within nine video frames (<45 min) (Appendix 1, Figure 6. continued on next page transition in which the orientation of the boundaries gained and lost are similar. Boundary orientations around T1 transitions are difficult to measure, because boundary length is small. The triangulation method provides us with a more robust measure of a bond orientation. As a proxy for cell boundary orientation, we use the orientation of the line connecting the centers of the two corresponding triangles ( Figure 6-figure supplement 1D). We calculated the proportion of connection losses or acquisitions occurring in different directions over the course of morphogenesis, and compared them with the distribution of connection angles in general. We observe that at early times, when T1-induced shear is mainly along the AP axis, connections are typically lost along the PD axis and gained along the AP axis ( Figure 7A-D). At intermediate times, when the correlation term is maximal, both lost and gained connections were oriented along the PD axis ( Figure 7B). This is consistent with the unusual T1 transitions associated with sliding cell rows. Finally, as T1 transitions begin to cause net PD shear, AP connections are preferentially lost and PD connections are gained ( Figure 7D). Thus, as T1 transitions shift from an AP-to a PD-oriented shear, they pass through an intermediate state where connection gains and losses are still oriented but do not cause shear. Interestingly, at the same time, the correlation term has maximal magnitude. This suggests that the correlation effects are related to these unusual patterns of T1 transitions at intermediate times.

Interplay of cell dynamics and tissue stresses during tissue shear
To investigate which cellular events depended on epithelial stresses, we quantified shape changes and cellular contributions to tissue shear in the dumpy ov1 mutant wing blade ( Figure 8A,B). As a complementary approach, we also studied these events in wings that had been subjected to laser severing ( Figure 8C-F and Figure 8-figure supplement 1). In dumpy ov1 wings, cells experience hinge contraction but the counterforces exerted by cuticle connections seem to be reduced (Figure 2A-A′′). Tissue shear is dramatically altered in dumpy ov1 mutant wings as compared to WT-instead of shearing in the PD axis, these wings shear on average along the AP axis ( Figure 8A,B). Examining the different contributions to tissue shear in dumpy ov1 wings shows that the rates of AP shear caused by T1 transitions and by correlation effects are similar to WT and persist for longer times. Thus, these processes are likely to be autonomously driven. By the end of the video, they cause more accumulated AP-oriented shear than in WT. Analogously, cell divisions cause more cumulative PD shear than in WT-consistent with the increased number of cell divisions in the dumpy ov1 mutant wing. In contrast, cell elongation during phase I causes less PD shear than in WT. Thus, PD-oriented epithelial stresses must contribute to PD cell elongation. Interestingly, the increase of cell elongation in the PD axis still exceeds the increase of elongation of the blade in the dumpy ov1 mutant wing. This suggests that autonomous cellular processes cause the residual PD cell elongation in dumpy ov1 mutant wings. Finally, PD shear by T1 transitions in phase II is smaller than in WT. This is not due to a premature cessation of T1's-indeed quantifying the rate of T1 transitions (regardless of orientation) shows that they occur at a higher rate and for a longer time than in WT wings (Figure 8-figure supplement 2). Rather, T1 transitions fail to orient as effectively with the PD axis in the dumpy ov1 mutant wing (see shear patterns in Figure 5-figure supplement 1B-B′′).
As a complementary approach, we asked how tissue shear and cell dynamics were altered in wings that had been subjected to laser severing. We first examined wings that had undergone laser severing before hinge contraction at the distal tip of the wing blade ( Figure 8C,D and Figure 8-figure supplement 1A,B). These wings experience mechanical conditions similar to those in dumpy ov1 wings in that they undergo hinge contraction but cannot attach properly to the distal cuticle. Consistent with this, the final shape of distally severed wings resembles that of dumpy ov1 wings. Furthermore, distally severed wings show similar tissue shear and cellular contributions to shear as dumpy ov1 mutant wings. These observations suggest that weakened connections to the distal cuticle are key to the altered wing shape of dumpy ov1 mutants.
To perturb PD stresses even more strongly, we severed the region between the hinge and the blade before hinge contraction occurred ( Figure 8E,F). These wings are not subjected to externally generated PD stresses at all, although they may still experience autonomously-induced stresses. We expected that cellular events that depended on externally generated stresses would be even more strongly perturbed in proximally severed wings than in dumpy ov1 or distally-severed wings. The proximally severed wing undergoes significant but short-lived PD shear during the molting process until 18 hAPF. Subsequently, PD shear stops and becomes even negative by 20 hAPF. Cells do not elongate in the PD axis as much as in WT wings, although PD cell elongation still exceeds PD tissue shear. Again, this suggests that only a fraction of PD cell elongation is normally caused by external stresses, and that autonomous cellular events must produce the residual PD cell elongation in proximally severed wings. PD shear due to cell division is larger than in WT (as it is in dumpy ov1 and distally severed wings) confirming that these divisions do not depend on external stresses. Furthermore, like dumpy ov1 and distally severed wings, proximally severed wings undergo greater AP shear resulting from correlation effects. Thus the cellular events underlying correlation effects produce even more shear when stresses are reduced. However later, T1 transitions fail to generate significant PD shear in proximally severed wings. The reduction in T1-dependent PD shear is much stronger than in either dumpy ov1 or distally severed wings, confirming that reorientation of T1 transitions in phase II is dependent on externally generated PD stress. At the very beginning of the video, PD-oriented connections are preferentially lost and AP-oriented connections are preferentially gained, consistent with the AP shear caused by T1 transitions at this time (Figure 8-figure supplement 2C-F). However, unlike WT, the preferential loss of connections never shifts towards the AP axis. Loss of connections remains biased towards the PD axis throughout the video-despite the fact that net shear caused by T1 transitions becomes very small. Shear caused by T1 transitions becomes small because the preferred orientation of gained connections gradually shifts from the AP to the PD axis. By the end of the video, both the assembly and disassembly of cell boundaries are preferentially oriented along the PD axis. These observations suggest that PD-oriented cell boundaries To calculate these effective proportions, we subtracted the angular distribution of all cell boundaries from the angular distribution of cell-cell connections that were lost or gained by cell rearrangements, revealing the orientation of cell boundaries with a disproportionate tendency to be lost or gained. Rose diagrams show angles of cell boundaries that are more likely to be gained (red) or lost (blue) at specific times during the video corresponding to important changes in cell dynamics: (A) 18.9 hAPF (peak of negative shear rate by cell rearrangements), (B) 21 hAPF (peak of correlation effects), (C) 21.5 hAPF (shear rate by cell rearrangements crosses zero) and (D) 24.5 hAPF (peak of positive shear rate by cell rearrangements). DOI: 10.7554/eLife.07090.023 , and wing blades severed either distally (C, D) or proximally (E, F) before hinge contraction (∼16 hAPF). Blue = total tissue shear, Red = shear due to T1 transitions, Green = shear due to cell elongation change, Orange = shear due to cell division, Magenta = shear due to correlation effects. Corresponding plots for WT wings (identical to those in Figure 6) are inset in the upper right corners of (A, B) for the purposes of comparison. Insets in left corners of (A, C, E) show the tracked regions of each wing in green at 32 hAPF. To plot accumulated tissue shear in laser ablated wings, we offset the calculated accumulated tissue shear (blue) by a value corresponding to the difference in blade elongation before ablation and at start of recording (see Figure 3-figure supplement 2). All videos were aligned in time by taking the histoblast nests fusion time as a reference at about 26.5 hAPF. All cumulated shear curves start at 16.2 hAPF, which is the earliest common time point registered in all compared videos, including the dist-sev#2 video shown in have a greater tendency to disassemble than those oriented at other angles, and that this is an autonomous, planar polarized feature of wing epithelial cells.
To disturb connections to the overlying cuticle without damaging the wing epithelium we disrupted apical extracellular matrix between the wing margin and the cuticle shortly before the onset of phase II, when this region becomes accessible (see Videos 8,9). When anterior connections are severed, the wing blade shears even more in the PD axis while the area decreases slightly (Figure 8-figure supplement 3). This suggests that these connections restrain the narrowing of the wing blade in the AP axis. Increased PD tissue shear is mainly a consequence of T1 transitions. Since severing anterior connections reduces AP stress while PD stress does not change, this supports the idea that T1 events during phase II are oriented by anisotropic stress.
When distal connections are severed, cells rapidly reduce their elongation and area ( Figure 8-figure  supplement 3). Surprisingly, the rate of PD shear caused by T1 transitions continues to increase as it does in WT at this stage. However, after about 4 hr the PD shear rate due to T1 transitions decreases prematurely. Since the analysis of dumpy ov1 and proximally severed wing shows that PD T1 transitions are stress dependent, this suggests that there is a time delay between the change in tissue stress and the resulting T1 transitions.
Taken together, the analyses of dumpy ov1 wing blades and laser-severed wing blades distinguish autonomously driven cellular processes from those induced by tissue stresses. PD cell elongation and PD-oriented T1 transitions clearly depend on tissue stresses, whereas PD-oriented cell divisions are autonomously driven. AP-oriented T1 transitions and a corresponding fraction of PD-oriented cell elongation are also driven autonomously, as are the cellular events underlying correlation effects.
The overall picture that emerges is that changes of wing blade shape arise due to force balances that involve stresses exerted at the boundary of the tissue, and internal tissue stresses. Boundary stresses are due to hinge contraction and to the resistance of extracellular matrix attachments to the cuticle. Internal tissue stresses are generated by cell autonomous processes, like T1 transitions, and by elastic cell deformations. Tissue mechanics depends strongly on elastic connections of the wing to the cuticle. We must now understand how these isotropic and anisotropic mechanical stresses in the tissue, combined with boundary stresses, lead to cell and tissue remodeling. The interplay between boundary stresses and forces generated in the tissue is complex and requires a physical approach. We now present a continuum mechanical theory to understand these force balances and to calculate both tissue and cell shape changes.
We first define tissue material properties, starting with elastic properties, adding cell autonomous stresses and tissue shear due to topological changes. We then introduce elastic linkers to the surrounding cuticle. Finally, we compare predictions of this theory to the experimentally measured cell and tissue shape changes and determine key biophysical parameters characterizing tissue material properties.

Relationship between tissue stress and cell elongation
We first investigated stresses present in the tissue. In an elastic tissue, tissue deformations stem from cell shape changes and cell elasticity is responsible for tissue elasticity. For small deformations, Hooke's law states that the mechanical stress in the material is proportional to its deformation. We write the isotropic part of the stress where P denotes two-dimensional tissue pressure, a and a 0 are cell area and preferred cell area respectively, and K is the area compressibility. The preferred cell area a 0 changes when cells divide. For small (a − a 0 )/a 0 , Equation 3 corresponds to Hooke's law. Equation 3 implicitly contains the active contribution to pressure, which influences the preferred cell area a 0 (see Appendix 2, 'Constitutive equation for the tissue stress').
We now focus on the anisotropic part of the stress σ ∼ , also called the shear stress. For simplicity, we write the elastic anisotropic stress in the form of Hooke's law σ ∼ e = 2K Q, where K is a shear elastic modulus and Q is the cell elongation. In this expression, the cell shape is isotropic in the absence of anisotropic stresses. However, in tissues, planar polarized cells may spontaneously elongate. Therefore, we postulate the following constitutive equation for the anisotropic tissue stress where ζ is a tensor that can be interpreted as a cell autonomous active stress related to spontaneous cell elongation.
To test whether Equation 4 accounts for anisotropic stresses present in the tissue, we further analyzed the circular laser ablations performed in different regions of WT and dumpy ov1 wings at 22 hAPF ( Figure 2I,I′). We defined a shear rate v ∼ cut that characterizes the anisotropic recoil of the circular cut boundary into an elliptic shape (see 'Materials and methods', Analysis of circular laser ablations). We plotted the projection of this shear rate on the PD axis v ∼ cut xx as a function of the projected average cell elongation Q xx ( Figure 9A). Cell elongation was determined as the average cell elongation within the corresponding region in unperturbed wings. We found that the anisotropic part of the shear rate varied linearly with cell elongation ( Figure 9A). The positive slope of this linear relation indicates that the shear modulus K is positive. In this argument, we use the recoil shear rate v ∼ cut as a measure proportional to tissue stress. A linear fit to this data also has a positive intercept, corresponding to a positive ζ xx in Equation 4. This implies that wing blade cells would spontaneously elongate along the AP axis in the absence of stress. Equivalently, in the absence of cell deformation (Q = 0) wing blade cells are exerting higher contractile stress in the PD axis than in the AP axis. The relative contributions of elastic stress and stress associated with spontaneous cell elongation can be quantified from this linear fit. The ratio of the intercept and the slope of the linear fit equals ζ xx /2K, which leads to the estimate ζ xx /K = 0.333 ± 0.003 in WT and ζ xx /K = 0.316 ± 0.004 in dumpy ov1 wing. Experimental data obtained in WT and dumpy ov1 wings fall on similar lines ( Figure 9A), suggesting that internal mechanical properties of the tissue are not perturbed in dumpy ov1 mutant wings. Since dumpy mutant cells are less elongated than those of WT (Figure 9-figure supplement 1A), their similar mechanical properties imply that they are under less anisotropic stress consistent with the loosened connections to the overlying cuticle.
If we perform the same analysis for the isotropic part of the tissue stress, we use the area expansion rate of the cut circle during recoil as a proxy for negative tissue pressure −P. Plotting this rate as a function of average cell area in a given region, we find that a 0 is smaller than cell area a, which reveals that the tissue is under contractile tension or negative P ( Figure 9B). This is consistent with the observed tissue area contraction in dumpy ov1 mutant and severed wings. However, parameter values a 0 and K =K cannot be reliably estimated by this method because the average cell area a varies too little, and because cell divisions may introduce heterogeneity in the preferred cell area a 0 .
Overall, laser ablation experiments indicate that there are two contributions to anisotropic stress in the PD axis. First, blade cells are elastically deformed along the PD axis, in response to hinge contraction and margin attachments. Second, these cells would tend to spontaneously shorten in the PD axis even in the absence of external stresses.  Model for the dynamics and the orientation of topological changes So far, we have characterized the elastic properties of the tissue. On long time scales, elastic stresses can be relaxed by topological changes of the cellular network generating tissue shear viscosity. We therefore now develop a model describing the average rate and orientation of tissue shear R due to such topological changes. Although R includes the shear caused by T1 and T2 transitions, cell divisions and correlation effects, it is dominated by T1 transitions (Figure 6A,B). We asked how the shear rate R due to topological changes is regulated in the wing blade. If cells in the tissue are elongated, this elongation could drive topological changes that give rise to the shear rate R = Q/τ r . Here, τ r can be interpreted as the time during which cell elongation relaxes. Indeed, τ r is also the time beyond which the tissue behaves as a viscous fluid with viscosity Kτ r , due to elastic stress relaxation by topological changes including cell rearrangements and cell divisions. In addition, our observations of lasersevered wings suggest that the response of R to cell elongation is not instantaneous, but follows with a time delay τ d of a few hours (see previous discussion and Figure 8-figure supplement 3), which also must be incorporated in the theory.
In addition to elongation-induced rearrangements, planar polarized tissues may undergo oriented topological changes even in the absence of cell elongation. For instance, topological changes could be driven autonomously by increased contractility of certain bonds. Taking into account all these contributions, we postulate the following relation between the shear rate R and cell elongation ( Figure 9C): Here, the delayed response to cell elongation is introduced via the time derivative of the shear R. This way of implementing a delay implies that the current value of cell elongation has the strongest impact on shear due to topological changes and the effect of recent values of cell elongation fades exponentially over time, disappearing beyond the time τ d . Shear driven by polarity dependent processes is characterized by the tensor λ. Thus, even if cells are not elongated (Q = 0), λ drives shear due to topological changes. For shear along the AP axis λ xx < 0, while for shear along the PD axis λ xx > 0.
In order to verify whether Equation 5 captures the dynamics of the shear created by topological changes during pupal wing morphogenesis, we plotted R xx vs the cell deformation component Q xx for different times ( Figure 9D). For an instantaneous response of R (no delay, τ d = 0) these points would fall on a straight line. For non-zero delay, the history of the process matters and these points follow a curve spiraling towards a fixed point (see Appendix 2, 'Effect of delay in topological changes'). We find indeed that experimental data for each WT wing follow a spiral, confirming the existence of a delay. A fit of the theory to the data allows us to estimate the coefficients in Equation 5 (τ d = (3.7 ± 0.9) h, τ r = (1.8 ± 0.6) h and λ xx = (−0.10 ± 0.04) h −1 ). Interestingly, λ xx is negative, indicating that polarity driven topological changes create AP shear-consistent with conclusions from laser ablation experiments. Polarity driven AP-oriented topological changes may be related to increased tension of PD-oriented as compared to AP-oriented cell bonds. This anisotropy of cell bond tension could account for the positive value of ζ in Equation 4 because it would tend to elongate cells in the AP-axis.
The fact that a constant negative value of λ xx accounts for the experimental data suggests that the tendency to undergo polarity-driven topological changes exists during both phase I and II. The transition from phase I to phase II occurs when cell elongation-driven topological changes begin to exceed the effect of polarity-driven topological changes. During phase I topological changes are largely polarity driven along the AP axis. Equation 2 then implies that these topological changes, together with hinge contraction, contribute to the buildup of cell elongation in the PD axis.
Once flows have stopped, the value of λ xx determines the final value of cell elongation. Indeed, Equations 2, 5 predict that in a steady state without shear flows (R = 0, dQ/dt = 0), the cell elongation is given by Q = −τ r λ. We can therefore use the cell elongation at the end of the videos, where the tissue is almost stationary, to estimate τ r λ xx . In other words, the final cell elongation should be independent of external stresses. Instead final cell elongation is governed only by the internal dynamics of the tissue, that is, the tendency to undergo AP-oriented topological changes λ xx as well as the stress relaxation time-scale τ r . Interestingly, examining different perturbed conditions, we find that the final cell elongation and hence λ xx depend on cell area (Figure 9-figure supplement 1B and Appendix 4, 'Area dependence of λ xx ').
Overall, we find that the shear created by topological changes is driven in part by cell elongation and in part by cell polarity-dependent processes. Furthermore, topological changes respond to cell elongation with a delay of several hours.
Equations 1-5 constitute a full theory for tissue mechanics taking into account cell shape changes. We now want to test whether this tissue description quantitatively accounts for cell and tissue shape changes during wing morphogenesis.
Cell shape changes and tissue flows during pupal development can be understood by a continuum mechanical model We ask here whether the tissue and cell properties described by our equations can quantitatively account for the observed changes in the shapes of hinge and blade as well as the elongation and topological changes of their constituent cells. We use a simplified description of the wing where the hinge and blade are represented by two rectangles that are attached to each other. These rectangles have the same areas as the hinge and blade and they undergo pure shear that is the average shear in the respective parts of the wing. To describe tissue flow, Equations 1-5, which characterize local tissue properties, have to be complemented by the condition of force balance. At the boundary of the tissue, elastic linkers connected to the cuticle impose external forces that have to balance tissue stresses. In our model the rectangles are therefore connected to an external frame by elastic elements ( Figure 10A and Figure 10-figure supplement 1). These external elements correspond to the extracellular matrix and the frame corresponds to the cuticle. The elastic elements provide resistance to extension of blade and hinge along the PD and AP axes (see Appendix 3, 'Boundary stresses').
In each rectangle, anisotropic stresses and shear are described by Equations 4, 5, with different tissue parameters in hinge and blade. Isotropic internal stresses in the blade are described by Equation 3. Preferred cell area a 0 changes rapidly when cells divide and is modulated by the cell contractility ζ (see Appendix 2, 'Constitutive equation for the tissue stress'). In the hinge, we do not solve the full mechanical problem but rather impose the observed hinge area contraction. We do so by adjusting the pressure in the hinge to match the observed hinge area.
We start with rectangles whose areas and aspect ratio are consistent with those of hinge and blade at 16 hAPF (see Appendix 3). The initial conditions for Q and R are the observed average cell elongation and the initial shear rate due to topological changes. To trigger cell flows in the model, we turn on active stresses ζ, ζ and polarity dependent topological changes λ. We solved Equations 1-5 in two rectangles that correspond to hinge and blade. Note that the pressure in the hinge is not determined by Equation 3, but by imposing hinge area. In the blade, we used the measured rates of cell division and T2 transitions to calculate the blade area changes (Equation 1).
We normalized all elastic moduli and friction coefficients to the blade shear modulus K. The solutions Q(t), a(t), v(t), and v ∼ ðtÞ of Equations 1-5 depend on a set of parameters characterizing the tissue both in hinge and blade (see Tables 1, 2). In addition, elastic coefficients describe constraints imposed by linkers at the boundary. Friction due to motion with respect to the cuticle is captured by a friction coefficient γ. Furthermore, to fully account for the observed flows, we found that we needed to add a contribution of tissue area viscosity η to Equation 3 (see Appendix 2, 'Constitutive equation for the tissue stress', Equation 52). We include this term in all subsequent calculations. Laser ablation of the extracellular matrix is introduced in our model by removing the elastic linkers at the tissue boundary at the side of ablation. We first analyzed WT wings in both unperturbed and mechanically perturbed conditions where only the extracellular matrix was ablated. Because blade tissue was not damaged by matrix ablation we expected that most parameters characterizing the hinge and blade would be identical to those in unperturbed wings. However, λ xx can differ in wings characterized by different cell areas, as noted above. We therefore used the values of τ r , τ d and λ xx determined by jointly fitting Equation 5 to the data (see Table 1 and Figure 9-figure supplement 2). In addition, we used experimentally determined values of ζ/K estimated from circular laser ablation experiments in WT wing blades ( Figure 9A). A similar fit is performed for the hinge and corresponding parameters are superscripted with an index H. This leaves 10 parameters unspecified. These characterize the isotropic stresses in the blade, the anisotropic stresses in the hinge, the stiffness of the elastic linkers at the boundary, and the coefficient describing external friction ( Table 2). We estimated these parameters by performing joint fits of Equations 1-5 to the quantified time dependence of tissue area changes and shear in unperturbed Within each rectangle, the tissue is subjected to cell-autonomous anisotropic and isotropic stresses ζ and ζ, and to topological changes driven by cell polarity-dependent processes λ. The complex elastic material connecting the wing to the cuticle is represented by AP-oriented elastic links (green and red springs on the cartoon) and PD-oriented springs (see Figure 10-figure supplement 1). In WT wings, the blade distal end is fixed, while it is free to move in the dumpy ov1 mutant. (B) (Left) Experimental (solid line) and theoretical (dashed line) time courses of tissue shear rate (blue curves), cell elongation change (green curves) and shear due to topological changes (red curves), in the blade and along the PD axis. (Middle) Experimental and theoretical time courses of cumulative tissue shear (blue curve), cell elongation (green curve) and cumulative shear due to topological changes (red curves), in the blade and along the PD axis. (Right) Experimental and theoretical cumulative relative blade area change. Model parameters were obtained by a fitting procedure to experimental data (Tables 1, 2). The continuum mechanical model recapitulates cell shape changes and tissue flow during wing morphogenesis. DOI: 10.7554/eLife.07090.033 and mechanically perturbed WT wings. We found that anisotropic stresses in the hinge are small compared to isotropic stresses and therefore we could set ζ H and K H to zero. A single set of the remaining 8 parameters accounted well for all time dependent data (see Figure 10B and  Table 2). It is remarkable that this very simplified model can account for the main features of the complex shape changes that occur during pupal wing morphogenesis. The choice of rectangles indeed captures the average shear and area change that determine tissue shape change. Displaying the time evolution of the rectangles within a time-lapse video of the pupal wing reveals the close agreement between the calculated and observed shape changes of hinge and blade (see Videos 10-13).
Our experimental data suggest that distal elastic connections to the cuticle are weakened in the dumpy ov1 mutant wing. We therefore asked whether we could fit the rectangle model described above using WT tissue parameters but allowing parameters describing connections to the cuticle to change. We removed distal boundary connections to the blade (Figure 10-figure supplement 1) as suggested by laser ablation experiments and allowed other parameters characterizing external linkers to change. We used these values and the measured value of ζ/K (which is almost the same as in WT, see Figure 9A) in fits of the full dynamics of the rectangle model to the dumpy ov1 mutant data (see Figure 10B and Figure 10-figure supplement 2). These calculations show that distal and lateral connections to the wing blade margin are weakened, as expected (see Table 2).
Overall, the continuum theory of epithelia outlined in Equations 1-5 together with appropriate boundary conditions recapitulates the dynamics of cell and tissue shape over 16 hr of morphogenesis. By fitting the theory to experimental measurements, we determine the values of tissue parameters characterizing intrinsic tissue time-scales, cell autonomous stresses, tissue elastic moduli and elastic moduli of external linkers. Using this theory to investigate mechanical parameters in a dumpy ov1 mutant wing, we confirm the existence of Dumpy-dependent elastic connections between the wing and cuticle. Thus, this theoretical approach is a powerful tool for studying how mutations influence specific aspects of cell and tissue mechanics.

Discussion
Developing tissues are active viscoelastic materials that generate and respond to mechanical stresses to change their shapes. How tissue material properties emerge from the properties and behavior of  their many constituent cells, and how these properties quantitatively account for tissue shape change is a major question in developmental biology. These questions have been difficult to answer because of both experimental and theoretical limitations. Methodologies for large-scale imaging and image analysis of entire tissues at the cellular level were insufficient, and direct measurement of forces and stresses in vivo have been rare (Hutson et al., 2003;Farhadifar et al., 2007;Rauzi et al., 2008;Bosveld et al., 2012;Bambardekar et al., 2015). Furthermore, theoretical approaches are required to understand how tissue shear emerges from different kinds of cellular processes. By developing new experimental and theoretical approaches, we have been able for the first time to account for the shape change of a developing tissue on the basis of its active and its viscoelastic properties. Our approach allows us to quantitatively understand how such tissue properties emerge from the interactions of its constituent cells. Our work is stimulated by previous approaches to tissue dynamics. For example, previous work on tissue tectonics highlighted the decomposition of tissue shear in cell shape changes and cell intercalation (Brodland et al., 2006;Graner et al., 2008;Marmottant et al., 2008;Blanchard et al., 2009;Butler et al., 2009). Our triangle method refines these ideas using a geometrical scheme that defines and distinguishes exact cellular contributions to tissue shear. These are contributions due to cell deformation, cell division, T1 and T2 transitions. Cell autonomous shear stress in wing blade of WT and dumpy ov1 are determined from circular laser cut experiments. Unperturbed and mechanically perturbed WT wings are first simultaneously fitted using results listed in Table 1.
Then, the dumpy ov1 wing is fitted keeping the values of hinge and blade tissue parameters the same as in WT. The effective anterior-posterior (AP) and PD elastic constants describe effects of external elastic elements providing resistance to changes in size of blade and hinge along the AP and PD direction. All quantities are normalized by the elastic shear modulus of the blade tissue K. Quantities containing spatial dimensions are also normalized by the initial length L 0 of the WT wing. Uncertainties reported for the parameters in this table (expect for the cell autonomous shear stress ζ xx ) were determined by the fit. Note that they do not reflect uncertainties arising from approximations made in the rectangle model (supplement section 4) and from pre-processing of experimental data (supplement section 1.6 This approach has also led to the discovery of a new correlation effect that contributes significantly to tissue shear. This collective effect can result from local correlations between cell elongation and tissue rotation, or between cell elongation and cell area changes. For instance, this correlation effect can be produced by relative sliding of neighboring cell rows ( Figure 4H). In our work we combine this approach with a novel continuum mechanical model based on patterns of cell elongation. Our model combines visco-elastic material properties with additional cell autonomous so-called active stresses. These active stresses contribute to force balances and also can drive topological rearrangements. This model allows us to quantitatively discuss tissue deformations that emerge from the dynamic interplay of externally applied forces, internally generated stresses and the resulting collective cell rearrangements. Our model is related to continuum models of active gels and tissues (Kruse et al., 2005;Brodland et al., 2006;Ranft et al., 2010) that can be obtained by coarse-graining cell-based models (Honda et al., 1984;Graner and Glazier, 1992;Drasdo et al., 1995;Chen and Brodland, 2000;Brodland et al., 2007;Farhadifar et al., 2007;Basan et al., 2011). However, our model differs from existing continuum models in several ways: it is directly based on measurable cellular contributions to tissue shear, it includes a new time scale corresponding to a delay for the generation of T1 transitions, and it introduces a cell-autonomous contribution to T1 transitions. Such cell-autonomous T1 transitions (Irvine and Wieschaus, 1994) may for example be generated by orientation-dependent contractility of cell boundaries (Bertet et al., 2004;Skoglund et al., 2008;Rauzi et al., 2010;Sawyer et al., 2011). Overall, our approach can bridge scales from individual cellular events to cell flows and large-scale tissue shape changes.
We find that the wing tissue shapes itself through patterned contractions and shear that occur in the context of patterned extra-cellular matrix (ECM) connections to a surrounding cuticular scaffold. Thus, by pulling on these connections, the tissue forces itself into the right shape. In dumpy mutants, these connections are weakened and the final tissue shape is dramatically altered as compared to WT (Waddington, 1940;Carmon et al., 2010). Dumpy plays a similar role in shaping the Drosophila trachea in embryonic development (Dong et al., 2014). The mechanical function of such ECM connections may be generally important during epithelial morphogenesis in vivo.
Our quantitative analysis has provided new insights into the biology of pupal wing morphogenesis. We have discovered that cell-autonomous planar polarized T1 transitions occur early, during the first phase of the process. Surprisingly, they occur along the AP axis and actively increase cell elongation in the PD axis. Why were these AP-oriented T1 transitions not detected in our previous analysis (Aigouy et al., 2010)? Comparing the rate of T1 transitions with the amount of shear they generate shows that the orientation of T1 transitions at this time is less focused than it is later when T1 transitions cause shear in the PD axis. However, because their overall rate is higher in phase I, a slight bias in their orientation causes significant AP shear. The fact that there is only a slight AP bias in these T1 transitions likely explains why they have not been detected before, and highlights the importance of the large-scale quantitative analysis that we employ here. What mechanisms might underlie AP-oriented T1 transitions that occur in phase I? Quantitative analysis of laser ablation data indicates that the preferred shape of wing epithelial cells in the absence of external forces is not isotropic, but rather elongated in the AP axis. Increased contractility of PD-oriented cell boundaries could explain this preference, and would also tend to favor AP-oriented T1 transitions. We note that components of both the Fat and Core planar cell polarity pathways are enriched on PD-oriented cell boundaries during phase I . These systems are known to regulate both cell boundary tension and Cadherin turnover (Classen et al., 2005;Rogulja et al., 2008;Mao et al., 2011;Bosveld et al., 2012;Warrington et al., 2013;Nagaoka et al., 2014), and it will be interesting to investigate their involvement in orienting T1 transitions at this stage. Our theoretical analysis suggests that APoriented T1 transitions are unlikely to strongly influence the final shape of the wing. However, they have a strong influence on peak cell elongation and final cell shape. What function could this serve during morphogenesis? One possibility is that AP-oriented T1 transitions influence the morphology of cuticular ridges on the adult wing surface. The cuticle of the adult wing is shaped in a reproducible pattern of ridges that form between cell rows (Doyle et al., 2008;Merkel et al., 2014). These ridges likely influence anisotropic mechanical properties of the wing during flight. Orderly packing of the wing epithelium into hexagons is necessary for the long-range order of these ridges , and their spacing presumably depends on the spacing of cell rows at the time the adult cuticle is secreted. AP-oriented T1 transitions could influence both these features. First, they are predicted to increase the final magnitude of PD cell elongation and may therefore influence ridge spacing. Second, increased cell elongation contributed by AP-oriented T1 transitions may help focus the direction of the subsequent PD-oriented T1 transitions and increase hexagonal order. Interestingly, PCP mutant wings are characterized by irregular cell packing geometry (Classen et al., 2005), and less long-range order in the pattern of cuticular ridges .
Our analysis has revealed unexpected properties of the cell elongation-dependent shift of T1 transitions towards the PD axis. Separately quantifying the orientations of cell boundary loss and cell boundary expansions during this shift shows that these two components of a T1 transition respond differently to cell elongation. The average orientation of cell boundary expansion shifts from the AP to the PD axis at 21 hAPF before the complementary change in orientation of cell boundary loss (Figure 7). This suggests that PD-oriented epithelial stresses promote expansion of cell boundaries along the PD axis at a lower threshold than is required to block shortening of PD-oriented boundaries. Thus, the mechanisms that operate to contract cell boundaries and to expand new ones respond differently to epithelial stresses. This difference means that cells tend to both gain and lose boundaries along the PD axis during the shift from AP to PD-oriented T1 transitions. Interestingly, this does not always simply represent the contraction and re-expansion of the same cell boundaries. Rather, contractions and expansions of different cell boundaries in the same direction occur as rows of cells move with respect to each other-a process that is associated with the maximal contribution of measured correlation effects (compare Figure 6A magenta line and Figure 7).
PD-oriented cell divisions also contribute to the shape change of the wing blade during pupal morphogenesis. Interestingly, while epithelial stresses promote cell division in cultured epithelial cells (Puliafito et al., 2012;Streichan et al., 2014), this does not appear to be the case in the wing. In fact cell divisions occur at a similar rate and for longer times in laser ablated and dumpy ov1 mutant wings. Thus, there are more cell divisions overall in situations where stresses are reduced. This suggests that the cell divisions that occur normally are not initiated by the stresses associated with hinge contraction but are controlled by other factors. Examining their pattern suggests that signals from veins may play a role (data not shown, [Garcia-Bellido et al., 1994]). It is also possible that hormonal signals such as ecdysone could control their timing. We do not yet understand why perturbations that reduce stress result in slightly more cell divisions. These extra divisions may occur by different mechanisms that respond to wounding.
Furthermore, the magnitude of PD shear caused by cell divisions is proportional to the rate of cell divisions, suggesting that cell division orientation is also unaffected by reducing PD stresses.
Video 13. Time evolution of the rectangles obtained from the rectangle model in a WT wing where the extracellular matrix was ablated anteriorly at the onset of phase II. Note that length of the hinge is increased by a correction term described in Appendix 4, 'Fitting of the rectangle model to cell and tissue shape in the hinge and blade'. DOI: 10.7554/eLife.07090.041 Thus, it may be that cell divisions in the pupal wing are autonomously controlled by planar polarized cortical cues. Alternatively, because the AP-oriented T1 transitions still cause some PD cell elongation in laser severed and dumpy ov1 mutant wings, it may be that residual cell elongation is still sufficient to orient cell divisions.
During development, tissues are shaped with extreme reproducibility. It has been estimated that the shape of the Drosophila wing is precise to within 1 cell diameter (Abouchar et al., 2014). How is such reproducibility achieved? The cell lineage in the wing is indeterminate, and both wing shape and wing area can be preserved in the face of a variety of developmental insults-including cell death and differential growth rates. This requires that all cells within the wing be able to sense its overall size and shape. How do cells acquire this information? Based on our findings, we would propose that stresses in the wing epithelium could be one cue.
Interestingly, the total area and shear of WT wings is more reproducible than would be expected from the variation in individual cellular contributions (Figures 3, 6). Thus, these cellular events influence and can compensate for each other. The fact that epithelial tension is required to maintain area homeostasis and to control tissue shear suggests a mechanism for this compensation. Variations in cell division, cell shape change, cell rearrangements and cell extrusions could be detected by cells as changes in epithelial tension. The ability of cells to sense and respond to tension could underlie this compensation and produce wings of reproducible sizes and shapes.

Fly strains and crosses
Flies were raised at 25˚C under standard conditions unless stated otherwise. Pupae were collected for imaging as described previously (Classen et al., 2008). Ecad::GFP flies (Huang et al., 2009) were used as control for all live imaging experiments. The dumpy ov1 mutation (Sturtevant et al., 1929) was recombined with Ecad::GFP for imaging, and for quantifying tissue flows and cell behaviors in a dumpy ov1 mutant (Bloomington, reference number 276). The Dumpy::YFP protein trap line (DGRC) was used to describe Dumpy distribution in the pupal wing. To inhibit cell divisions in the pupal wing while imaging it, we used the thermosensitive allele cdc2 E1-E24 of cdk1 that we recombined with Ecad:: GFP (Stern et al., 1993). Cells expressing two copies of cdc2 E1-E24 and shifted to 30˚C arrest in G2 phase just prior to entering mitosis (Weigmann et al., 1997).

Long-terms time-lapse imaging
Pupae were prepared for live imaging as previously described (Classen et al., 2008). Long-term time-lapses were acquired with a Zeiss spinning disk microscope driven by the Axiovision software and equipped with an inverted Axio Observer stand, a motorized xyz stage, and a Zeiss LCI Plan-Neofluar 63× 1.3 NA Imm Corr lens associated to an objective heater set to 25˚C. The fly pupa was placed in a temperature-controlled chamber set to 25˚C and equipped with a humidifier to prevent desiccation. Images were recorded with an AxioCamMR3 camera (2 × 2 binning). Laser power was measured through a 10×/0.45 NA lens using a power-meter (PT9610, Gigahertz-Optik). A power of 0.980 mW was found to be optimal to prevent noticeable bleaching during 24 hr of continuous acquisition with an exposure time of 265 ms. Briefly, the dorsal cell layer of the pupal wing was scanned within 5 min in the (x, y, z) dimensions, over 24 overlapping positions of about 30 z-sections each. This scan procedure was continuously repeated about 260 times to cover more than 20 hr of development. Custom Fiji macros helped keeping the tissue in focus in a semi-automated way.

Data handling and image processing
We benefited from the flexible architecture provided by our computer department to handle several TB of data using custom unix bash scripts to archive, compress and store data on tape and retrieve them easily. Image pre-processing steps were mostly performed using custom Fiji macros called from a master bash script to enable parallelization using GNU parallel (Tange, 2011). A low-pass filter was first applied to images to remove high frequency noise. For each z-stack, the signal of highest contrast was projected using a C++ algorithm as previously described . Tiles were stitched using Fiji (Preibisch, 2009;Schindelin et al., 2012). Stitched images were then loaded in Packing Analyzer v8.5 (PA8.5) for cell edge detection by using a seeded-watershed algorithm. To facilitate cell tracking, cross-correlation of subsequent images was performed in PA8.5 to calculate local tissue displacement beforehand. Resulting cell-tracking masks were parsed using a custom C++ parser that extracted cell shape properties, cell lineages and cell neighbor relationships for further storage in a SQLite relational DB. Topology was kept in the DB by ordering cell neighbors counter-clockwisely around each cell, namely by following a directed path of cell-cell junctions around each cell. One DB file was created per video. All DB queries were carried out using R (R Development Core Team, 2014) or Python (www.python.org ;Hunter, 2007;Pérez and Granger, 2007;van der Walt et al., 2011).
We restricted our analyses to a subset of cells that were trackable throughout the entire course of the video, excluding cells that became visible only at later stages due to the apposition of the two cell layers. At 32 hAPF, most cells of the dorsal cell layer are visible. Therefore, we manually drew regions of interest (ROIs) on 32 hAPF wings using the veins and hinge-blade interface as landmarks, where the Ecad::GFP signal is more intense. These ROIs include the blade region delimited by the hinge-blade interface and the wing margin, but also the hinge-blade interface itself. In addition, we defined a triangular portion of the hinge delimited by the hinge-blade interface and the most anterior sensory organ located in the bulk of the hinge. In other hinge regions, cells were too small and elongated to be consistently tracked. Starting from these ROIs at 32 hAPF, we developed a backward-tracking algorithm that uses information about cell divisions and cell extrusions to reconstitute the corresponding group of cells at start of recording (usually ∼16 hAPF), discarding margin cells that were not present through the course of the video. In the resulting set of cells, cells in contact with the wing margin at start of recording (about 16 hAPF) were not perfectly segmented and therefore further discarded. This gave rise to a group of cells that is fully consistent in time.
Since part of the hinge could not be analyzed at cellular resolution, we complemented our analysis of the small portion of the hinge by using particle image velocimetry (PIV) (Adrian and Westerweel, 2011) to extract information about the entire hinge deformation as described thereafter. PIV was implemented by Benoit Aigouy and described elsewhere .

Laser ablations
Laser ablations were performed using an ultraviolet laser microdissection apparatus as described elsewhere (Grill et al., 2001). For tissue severing, a C-apochromat-40×/1.2 water immersion lens was used to focus the beam along AP-oriented line segments to ablate both dorsal and ventral cell layers. Subsequent long-term imaging of the whole wing was then carried out on our dedicated Zeiss spinning disk. Extracellular-matrix ablations were done using a 63×/1.4 oil lens to focus the beam between the tissue and the cuticle and to cut over 10 μm in depth. Circular cuts were also conducted with the 63×/1.4 oil lens but over a depth of 5 μm. Each circular cut experiment was repeated on five distinct pupae.

Analysis of circular laser ablations
Fly wings expressing Ecad::GFP were used in all circular cut experiments. Circular cuts were performed to disconnect a small subset of cells present inside the circle from the rest of the tissue. The interface of the tissue with the ablated circular region was always visible and well spatially defined after laser ablation. Therefore we quantified its displacement to calculate the initial velocity gradient that describes the immediate tissue deformation in response to the ablation. To do so, we first fitted an ellipse to this manually segmented interface, at 50 s after ablation. The major and minor axes of the ellipse define the orientation of two orthogonal kymographs, each intersecting the ellipse in two points corresponding to two opposite sides of the interface between the tissue and the ablated region. Each kymograph depicts the displacement as a function of time of the two opposite sides of this interface, which were manually segmented using sub-pixel resolution ROI. Thus, we quantified the displacement of four points corresponding to the intersections of the minor and major axes with the ellipse (see Figure 2-figure supplement 3). This procedure was semiautomated in a custom Fiji macro. The initial response of the tissue reflects the orientation and amplitude of tissue stresses. This can be captured by the initial velocity gradient within a few seconds after the 4 s ablation. Indeed, the initial displacement normal to the cut boundary was approximately linear as a function of time within the first 5 s after ablation. A linear fit to the data provided the initial normal velocities V || and V ⊥ along the major and minor axis of the ellipse, respectively. We then define the velocity gradient tensor in the coordinate system of the ellipse and rotate it into the image coordinate system: where θ is the angle between the major axis of the ellipse and the PD axis (x axis). The radii r || and r ⊥ are the half lengths of major and minor axes of the ellipsoidal shape at the time the velocities are determined. This velocity gradient can be decomposed into trace and a traceless part as: where v cut kk = 2C is the isotropic expansion rate, v ∼ cut xx is the shear rate projected onto the PD-axis. While in the main text, tensors are denoted as bold face symbols, in the supplement we often use an explicit index or matrix notation for tensors. Latin indices i, j denote the x and y coordinates of a cartesian coordinate system of the tissue.
Circular cuts were performed in nine different regions of the wing where cell elongation differed ( Figure 2I,I′). Due to the fast imaging settings optimized to catch the initial velocity gradient, the image quality was not sufficient for cell segmentation. Therefore we used wings of different animals at the same stage (22.5 hAPF) to estimate cell elongation Q xx for each of the nine ablation regions that were easily located using morphological landmarks such sensory organs and veins (see Figure 9A). Average cell area a cell was estimated using a similar approach (see Figure 9B).
We find that the shear rate projected onto the PD axis v ∼ cut xx increases with cell elongation. We find that data of v The correlation between the isotropic expansion rate v cut kk and the logarithm of cell area is less apparent. However, following Equation 8, we performed a linear fit to the data and found the following best fit parameters: c WT: 1 2 v kk ½h −1 = ð0:008 ± 0:007Þlnða cell =a ref Þ + ð0:017 ± 0:002Þ c Dp: 1 2 v kk ½h −1 = ð0:005 ± 0:008Þlnða cell =a ref Þ + ð0:015 ± 0:002Þ Parameters obtained in isotropic expansion rate fits have high uncertainties and we do not use them in the rest of the paper.

Measurements of wing dimensions
In order to compare theory to the experimental data we need to approximate the observed hinge and blade shapes by rectangles. We first show how a characteristic height and length can be associated to an arbitrary two-dimensional shape. We then specify how hinge and blade regions were selected to obtain the associated height and length. Note that for the purpose of estimating the wing dimensions, we do not used tracked regions of a subpart of the wing as in 'Data handling and image processing', but we use all available segmented data. Finally, we introduce corrections that account for changes of the visible part of the tissue in the field of view.

Rectangle approximation to an arbitrary shape
To approximate an arbitrary tissue shape with a rectangle, we have to define a height h and a length L along the x direction that properly represent that region. To do this, we take into account the region area and elongation. First, we impose the condition that the rectangle representing the shape has the same area: hL = A: (8) Additionally, we can use the elongation tensor Q t defined in Appendix 1, 'Characterization of wing blade anisotropy' to obtain the rectangle aspect ratio. We then define the aspect ratio of the rectangle to be: These two conditions finally give us:

Boundary between hinge and blade
To separate the wing into hinge and blade, we used a tracked region of cells on the boundary between the hinge and the blade (Data handling and image processing). At each time point we fit a fourth order polynomial through the coordinates of cells in this boundary region. This curve is used as a boundary separating hinge and the blade.

Extrapolation of region dimensions
During the early time of recording, the tissue flows in and out of view. In order to correct for this, we performed the following measurements: c We first obtain h S and L S from the shape of the segmented hinge and blade as described above. c We then calculate the average shear in the hinge and blade from PIV measurements (see Data handling and image processing). From the average shear, an other estimate of the hinge and length can be obtained, up to unknown factors β L and β h : c After a time t f = 22.7 hAPF for WT and t f = 24.5 hAPF for dumpy ov1 , no flow is visible at the boundary of the tissue, and we therefore expect the two measurements of height and length to coincide. We therefore determine the coefficients β L and β H by minimizing the difference between L PIV and L S and h PIV and h S for t > t f . A good agreement was indeed found for the two measurements for t > t f .
Wings where the extracellular matrix was ablated distaly and anteriorly start much later than unperturbed WT and dumpy ov1 and we do not observe significant inflow of the tissue. The same procedure was applied for these wings, with t f the first time where measurements were started.
Here we define the average pure shear rate of a tissue region and its decomposition into cellular contributions (Merkel, 2014). The cellular network of the tissue region is triangulated ('Triangulated cellular networks'). Averaging triangle deformations over all triangles within the region yields expressions for the large-scale shear rate and the shear rate contributions by different cellular processes ('Decomposition of the large-scale tissue shear rate'). The determination of spatial patterns of shear ( Figure 5, Figure 5-figure supplement 1) and rotation ( Figure 6C,D,F,G) are discussed in 'Spatially resolved shear patterns' and 'Power spectrum of local tissue rotation'.

Triangulated cellular networks
Here, we first describe the tiling of the cellular network in detail before defining deformation and elongation of a single triangle.

Triangulation procedure
To define a triangulation of a given cellular network, we consider all vertices of the network that are not located at the margin. Each of these inner vertices abuts at least three cells. If a given inner vertex m abuts exactly three cells α, β, and γ, we define a single triangle m that has its corners at the centers R α , R β , and R γ of the abutting cells. Here, R α denotes the area center of cell α: R α = ð R rd 2 rÞ=A α , where the integral runs over the projected area A α of cell α.
For an inner vertex m that abuts N > 3 cells α 1 ,…,α N , we define N triangles m 1 ,…,m N . The corners of each of these triangles are defined as follows. One corner of each triangle is defined by the average position C = ðR α 1 + … + R α N Þ=N of all N cell centers. The other two corners of a given triangle m i with i = 1,…,N are defined by the cell centers R α i and R α ði + 1Þ (R α ðN + 1Þ corresponds to R α 1 ). Applying these two rules to all inner vertices, we obtain a tiling of the whole network without gaps or overlaps, leaving out only a small stripe at the margin of the network, which is ca. half a cell diameter wide.

Triangle deformation
We characterize the deformation of a given triangle m occurring during a time interval Δt by a tensor M m ij . The Latin indices i, j,… denote dimension indices of vectors, tensors, and matrices. For matrices, the first index is the row index and the second index is the column index. To define M m ij , we denote the initial and the final positions of the triangle corners by R 1 , R 2 , R 3 and R′ 1 , R′ 2 , R′ 3 , respectively. Then, the triangle side vectors before and after the deformation read ΔR βα = R β − R α and ΔR′ βα = R′ β − R′ α for α, β ∈ {1, 2, 3}. Now, there is a unique affine transformation defined by the tensor M m ij that maps the initial side vectors to the final side vectors: for any combination α, β ∈ {1, 2, 3}. Equal indices on one side of the equation imply a summation over the two spatial dimensions. From standard linear algebra follows that the tensor M m = ðM m ij Þ has the following unique matrix representation: In this equation the dot denotes the matrix product.
Based on the affine transformation tensor M m ij , we compute a discrete velocity gradient tensor v m ij for the triangle m: Here, the δ ij symbol denotes the Kronecker symbol. Note that the tensor M m ij enters transposed into this equation. This leads to a discrete velocity gradient v m ij that approximates the local velocity gradient v ij = ∂ i v j in the continuum limit, where v(r) denotes the velocity field.
Finally, we decompose the discrete triangle velocity gradient into trace v m kk , traceless, symmetric part v ∼ m ij , and antisymmetric part ω m : with the generator of rotations defined by The trace v m kk represents the relative area expansion rate of the triangle, the traceless, symmetric part v ∼ m ij represents a discrete triangle shear rate, and the antisymmetric part ω m represents a triangle rotation rate.

Triangle elongation
We characterize triangle shape by an affine transformation tensor S m ij with respect to a universal reference triangle. This reference triangle is equilateral, and it has a fixed orientation and a fixed area A 0 . The corners of the reference triangle are denoted by C 1 , C 2 , and C 3 . At some fixed time point, we denote the corners of a given triangle m by R 1 , R 2 , and R 3 . Like above, we denote the respective triangle side vectors by ΔC βα = C β − C α and ΔR βα = R β − R α for α, β ∈ {1, 2, 3}. Then, we define the tensor S m = ðS m ij Þ via the following matrix representation: The so-defined tensor S m ij maps the side vectors of the reference triangle to the side vectors of the triangle m, analogously to Equation 14.
To extract shape properties from the tensor S m = ðS m ij Þ, we use the following polar decomposition (Kaballo, 2013): Here, A m denotes the area of triangle m. This equation uniquely defines the symmetric, traceless 2 × 2 tensor Q m , where the exponential of a tensor is defined by the Taylor series of the exponential function, and Rot(Θ m ) denotes the matrix describing a counter-clockwise rotation by the angle Θ m . The symmetric, traceless tensor Q m denotes the elongation of triangle m and the angle Θ m denotes a triangle orientation angle.

Decomposition of the large-scale tissue shear rate
In this subsection, we consider the deformation of the cellular network during a finite time interval Δt between two consecutive video frames. We denote the observed states of the cellular network corresponding to these two video frames by O n and O n + 1 , where n is an integer. We explain how we quantify the overall tissue shear rate between these two video frames and how we decompose it into cellular contributions. To separate the shear by topological changes from the shear by pure cellular deformation, we first introduce virtual intermediate network states between the two observed states O n and O n + 1 .

Intermediate network states
Intermediate network states are introduced for the following reasons. First, topological changes occur instantaneously, and as a consequence of this, no tissue shear can occur during a topological transition. Second, topological transitions involve the addition/removal of triangles to/from the network ( Figure 4L-O). Third, we cannot determine the exact time at which a topological transition occurred during the interval Δt between two observed network states O n and O n + 1 .
To separate different types of changes between O n and O n + 1 , we introduce three virtual intermediate network states denoted I 1 , I 2 , and I 3 . This is illustrated by the following scheme: We first define the intermediate state I 3 by starting from the observed state O n + 1 . To this end, we fuse in I 3 all pairs of daughter cells in O n + 1 . This corresponds to undoing the cell divisions occurring during the time interval Δt. The position of the mother cell center in I 3 is defined as the midpoint between the daughter cell centers in O n + 1 . We then define the intermediate state I 1 from the state O n by removing all cells that will undergo a T2 transition during the time interval Δt. As a consequence, I 1 and I 3 contain the same set of cells. However, the topology and the cell center positions are different between these two states. Finally, we construct the state I 2 from the state I 1 by displacing each cell center to the position it attains in state I 3 . Thus, the intermediate state I 2 has the same topology as the state I 1 , but cell center positions correspond to those in I 3 . As a consequence of all these definitions, all T2 transitions occur between states O n and I 1 , any deformation (without topological changes) occurs between I 1 and I 2 , all T1 transitions occur between I 2 and I 3 , and all cell divisions occur between I 3 and O n + 1 (compare scheme (21)). Furthermore, since topological transitions occur instantaneously, the pure deformation between states I 1 and I 2 coorresponds to the time interval Δt.
The scheme in (21) is the basis for the computation of the shear rate and all cellular contributions to it. In the following, we will separately discuss shear in the absence of topological transitions (i.e., between the states I 1 and I 2 ) and the shear by topological transitions.

Large-scale shear in the absence of topological transitions
The overall tissue shear rate v ∼ ij between the states I 1 and I 2 is computed as the average triangle shear rate: In this equation, the average corresponds to an area-weighted average over all triangles, where the area weights are taken from state I 1 . More generally, for any triangle-based quantity q m , we define such an average by Here, the sum runs over all triangles m of the network, A m X denotes the area of triangle m in state X ∈ {O n , I 1 , I 2 , I 3 , O n + 1 }, and A X is the total area of the triangulation in state X: The index X is necessary, because the triangle areas generally change between two states. The average quantity q m may be evaluated for a network state Y different from X if both contain the same set of triangles. Interestingly, the large-scale shear rate v ∼ ij can equally be computed from margin displacements alone (Merkel, 2014;Merkel et al., in preparation). Furthermore, since there is no shear occurring during topological transitions, the average shear rate v ∼ ij between the intermediate states I 1 and I 2 also corresponds to the average shear rate between the observed states O n and O n + 1 .
The overall tissue shear rate between the intermediate states I 1 and I 2 can be expressed in terms of the average triangle elongation change during that time: Here, is a corotational term and we have introduced the correlation contribution to shear D ij , which can be approximated as AEω m Â Q kj ðI 1 Þ + Q kj ðI 2 Þ Ã ae I 1 − AEω m ae I 1 AEQ kj ðI 1 Þ + Q kj ðI 2 Þae I 1 : In Equation 26, the first term captures effects of area expansion rate correlation with cell elongation, while the second term corresponds to correlations of rotation rate with cell elongation ( Figure 6E). Here, we neglected higher order terms in the deformation and in the cell elongation. An exact decomposition of v ∼ ij in the limit of small time intervals Δt can also be obtained (Merkel, 2014;Merkel et al., in preparation).

Contributions to shear by topological changes
Between the two states O n and O n + 1 , we define the contributions of topological changes to the shear rate as follows: The total change of the average triangle elongation between the observed video states is ΔQ ij ðO n ; O n+1 Þ = AEQ ij ðO n + 1 Þae O n + 1 − AEQ ij ðO n Þae On . Putting everything together, we obtain Equation 2 from the main text: with R ij = T ij + C ij + E ij + D ij . The corotational derivative DQ ij /Dt is defined by and the quantities v ∼ ij , J ij , and D ij are defined in the previous subsection.

Spatially resolved shear patterns
To visualize the spatial patterns of shear contributions, we divided the entire wing into a fixed grid of non-overlapping boxes, sized 26 μm × 26 μm. Then, each triangle is assigned to the box that contains the triangle center. This allows for the computation of average triangle quantities within each of these boxes. As in the previous section, for a given time interval Δt between two consecutive video frames O n and O n + 1 , we created intermediate states I 1 , I 2 , and I 3 (compare scheme (21)). Note that the set of triangles associated to a given box depends not only on the time interval n, but also on the state considered (O n , I 1 , I 2 , I 3 , or O n + 1 ).
The five shear patterns shown in Figure 5 and Figure 5-figure supplement 1 were computed as follows. The total shear rate v ∼ ij , the corotational term J ij , and the correlations D ij in a given box were computed using Equations 22, 25, 26, respectively, where for each average appearing in these expressions, only those triangles associated to the given box in state I 1 were included. The remaining quantities were computed as follows: Here, b denotes the box index. The averaging bracket AEQ m ij ðXÞae b denotes an average of the triangle elongation in state X over all triangles m that are assigned to box b in state X with area weights from state X, where X ∈ {O n , I 1 , I 2 , I 3 , O n + 1 }: Here, the sums runs over all triangles m that are associated to box b in state X. The contributions by T2 transitions were negligible.
For the pattern of corotational cell elongation change, we neglected convective contributions.
To evaluate, how much this affected the shown patterns, we compared our results to an alternative implementation that tracks patches of tissue over time. Here, the fixed grid was only used for the initialization of the patches at the beginning of the video. For this alternative implementation, including a convective contribution is not necessary since each patch of tissue is followed over time. However after all, visual comparison of both implementations yielded a good agreement (not shown). Thus, the patterns in Figure 5 and Figure 5-figure supplement 1 were not essentially affected by neglecting the convective term.

Power spectrum of local tissue rotation
Here, we explain how we computed the spatial power spectrum of the local tissue rotation rate shown in Figure 6D. For a given time interval between two consecutive video frames, the triangle rotation rates ω m are computed as described by Equation 17 and 'Decomposition of the large-scale tissue shear rate'. This yields a piecewisely constant field of local rotation rates ω(r): at a given position r that lies within a triangle m, the value of the field is defined by the triangle rotation rate ω(r) = ω m .
The power spectrum jF ½ωj 2 ðqÞ of this field is then computed as the squared norm of following Fourier transform: The integration in this expression was carried out over a squared region with a side length of L x = L y = 800 pixels (164 μm) that is located in the center of the wing, covering parts of the longitudinal veins 2 to 4. The product in the exponential is a dot product. The components of the wave vector q = (q x , q y ) could take values q x = 2πn x /L x and q y = 2πn y /L y , where n x and n y are integers with −N/2 ≤ n x , n y < N/2. We set N = 100. The function w 2D (r) is a window function defined on the squared region. With the position r = (r x , r y ) measured relative to the bottom-left corner of the squared region, it reads w 2D ðrÞ = w BH r x L x w BH r y L y : Here, w BH (z) denotes the one-dimensional Blackman Harris Window, which is defined by w BH ðzÞ = a 0 − a 1 cosð2πzÞ + a 2 cosð4πzÞ − a 3 cosð6πzÞ; with a 0 = 0.35875, a 1 = 0.48829, a 2 = 0.14128, and a 3 = 0.01168. The scalar W is defined by It corresponds to the average value of the squared window function W = AEw 2 2D ae. W is used in the normalization factor in Equation 36 in order to ensure that the second moment of ω corresponds to the integrated power spectrum in the limit N → ∞ and in the absence of any correlation between the averaging window w 2D (r) and ω(r) such that AEw 2 2D ω 2 ae = AEw 2 2D aeAEω 2 ae. Here, the average of any spatially varying function q(r) is defined by the area average over the squared region: AEqae = ð R qðrÞd 2 rÞ=L x L y .
We numerically evaluated the Fourier transform using the following expression: where we used the fact that the window function w 2D (r) varied only little over the size of a single triangle. Here, the vector r m c denotes the center position of triangle m and the integral is carried out over the area of triangle m. For a triangle m with corners r 0 = (x 0 , y 0 ), r 1 = (x 1 , y 1 ), and r 2 = (x 2 , y 2 ) in counter-clockwise order, the integral over triangle area A m can be expressed as: In this expression, indices are taken modulo three such that i = 3 corresponds to i = 0 and i = −1 corresponds to i = 2. Figure 6D results from averaging the power spectra jF ½ωj 2 computed from the 11 time intervals in between 12 consecutive video frames. From Figure 6D, we can extract the width of the observed rotation bands. This width corresponds to the position of the dominant mode in the power spectrum, which we observe at q peak ≈ (0, ±0.3d 0 /2π), where d 0 = 4 μm denotes the typical cell diameter at 21 hAPF. This mode corresponds to a wave length of ca. 3d 0 , and thus to a width of the bands of ca. 1.5d 0 .

Role of T1 transitions in the correlation-induced shear
We have divided all triangles in the blade region in a WT wing into two groups. The first group contains all triangles that will disappear due to a T1 transition during the next 9 frames of the video, corresponding to about 45 min. The second group contains all other triangles. The correlation contribution to the total tissue shear can be approximately written as Here, A T1 and A T1 are the areas and D T1 ij and D T1 ij are the correlation contributions to shear of the first and the second group of triangles, respectively. The area of the blade region is denoted by A.
In the inset of Figure 6I we show that the area of the first group of cells A T1 (green line) is only a small fraction of the total blade area A (blue line). On average, the ratio is A T1 /A = 15.3%. Despite this, we show in Figure 6I that the measured contribution of the first group projected R ij = T ij + C ij + E ij + D ij ; where T ij is the shear rate tensor created by T 1 transitions, C ij is the shear rate created by cell division, E ij is the shear rate created by cell extrusion, and D ij denotes a tensor accounting for correlation effects.

Constitutive equation for the tissue stress
The tissue stress can be decomposed into an anisotropic traceless part, σ ∼ ij , and an isotropic part corresponding to a two dimensional tissue pressure P: In an elastic tissue, the anisotropic part of the stress and the pressure depends on cell deformation. We therefore write the following constitutive equations for the stress: Here, K and K are a shear and an area elastic modulus. The traceless symmetric tensor ζ ij ∼ q ij describes a cell autonomous active stress which can also be interpreted as a spontaneous cell elongation. Here, η is a cell area viscosity, related to the dissipation associated with changes of cell area. Note that this viscosity was not included in Equation 3 in the main text, but is used in the fitting procedure described in 'Fitting of the rectangle model to cell and tissue shape in the hinge and blade'.
Cell division is associated with a corresponding rapid reduction of cell area ( Figure 3C,D) suggesting that the preferred cell area a 0 changes during cell division. We thus express the preferred cell area a 0 as: where a 0 is a cell intrinsic preferred cell area that changes rapidly due to cell division and ζ corresponds to an isotropic contractile tension. The contribution a 0 follows the equation which corresponds to halving of a 0 for each round of division. In addition to Equation 59 we can choose the value of a 0 at an arbitrary time point. For simplicity we choose a 0 ðt 0 Þ = aðt 0 Þ.
With the decomposition of the preferred cell area Equation 58 and using the isotropic strain rate decomposition Equation 50, the isotropic pressure can be expressed as where A = A 0 expð R t t 0 vdtÞ is the change of tissue area between t 0 and t, A 0 is the tissue area. Equation 60 clarifies that the definition of the reference area Equation 58 corresponds to a decomposition of the tissue pressure into an elastic part (first term in the right-hand side of Equation 60), a part corresponding to the stress generated by loss of cells by extrusion, isotropic contractility ζ driving tissue contraction or expansion, and a viscous part associated with the rate of cell area change.

Constitutive equations for the shear created by topological changes
The shear created by topological changes, R ij , is a traceless symmetric tensor, whose value must depend on other traceless symmetric tensors in the tissue. Physically, topological changes occur preferentially along a direction set by an external cue. In the hydrodynamic description we propose, such a cue can arise from cell elongation Q ij or cell polarity q ij . In the spirit of linear response theory, we write that R ij depends linearly on the cell deformation tensor Here, D/Dt is the corotational convected derivative (see Equation 53) and τ r is a timescale, characterising the rate of shear corresponding to topological changes induced by cell deformation. Note that for simplicty, we neglect nonlinearities arising from the corotational derivative in Equation 61 in the main text and in these supplements. The traceless symmetric tensor λ ij ∼ q ij describes the effect of polarity-induced topological changes. Tissue stability requires τ r > 0. The time τ d characterizes the delay in the response of topological changes to cell elongation and cell polarity.

Effect of delay in topological changes
We discuss here the delay term for topological changes introduced in Equation 4 of the main text. Neglecting non-linearities introduced by convection and rotations, the shear decomposition Equation 52 and constitutive Equation 61 can be written: which appear as a dynamical system for the evolution of Q ij and R ij . This dynamical system evolves towards the fixed point R Assuming a constant shear rate v ∼ ij , deviations around the fixed point follow the equation The dynamics of the system around the fixed point depends on the timescales τ r and τ d . For a short enough delay τ d < τ r /4, the eigenvalues of the matrix in Equation 64 are real negative, and the system relaxes exponentially to the stable fixed point. For a long enough delay τ d > τ r /4, the matrix has complex eigenvalues −1=2τ d ± i ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4τ d τ r − τ 2 r p =2τ r τ d . In that case, the system behaves as a damped oscillator and relaxes to the fixed point following a spiral in the phase space R ij , Q ij . The system relaxes towards the fixed point on a timescale 2τ d , with a period of oscillations A relaxation resembling a damped oscillator in the phase space R xx , Q xx is indeed observed experimentally ( Figure 9D and Figure 9-figure supplement 2). Fitting the relaxation in the phase space R xx , Q xx to Equations 62, 63 yields a delay τ d of the order of a few hours, thus playing a significant role during pupal wing development.
In this section, we use the theoretical framework described in the previous section to obtain a simplified description of the morphogenesis of the pupal wing. We represent the hinge and blade by rectangles, connected by elastic springs to an external fixed frame, corresponding to the cuticle ( Figure 10A and Figure 10-figure supplement 1).

Geometry and boundary conditions Geometry
For simplicity, we choose to represent the hinge and the blade by two rectangles. The x axis labels the proximal-distal direction, while the y axis labels the anterior-posterior direction. The two rectangles only deform uniformly, and are therefore subjected to a homogeneous velocity gradient with components v xx and v yy . The position of the interface between the hinge and the blade is denoted x HB , the heights of the hinge and of the blade h H and h. In the WT situation, we consider the total length of the wing (including the hinge and the blade) to be constant and equal to L.

Velocity
Because the proximal distal velocity field v x has to be continuous at the contact point x HB (Figure 10-figure supplement 1B), and cancels at the proximal and distal ends, the gradient of flow in the hinge and in the blade can be directly obtained as a function of the velocity of the interface v HB = dx HB /dt: Similarly, the gradient of flow in the y direction can be related to the change of heights of the hinge and of the blade:

Boundary stresses WT wing
From force balance, the stresses acting in the hinge σ H xx and in the blade σ xx , integrated along the height, have to be balanced at the interface: Here, γ is an effective friction coefficient, reflecting dissipative processes external to the tissue and limiting its velocity. The two elastic moduli constraining the length of the hinge and blade are denoted k H PD and k PD . They reflect external attachments of the tissue to the cuticle (see Figure 10A and Figure 10-figure supplement 1A). The initial PD lengths of the hinge and the entire wing are denoted x 0 and L 0 , respectively.
In addition, we consider the effects of the elastic connections at the anterior and posterior boundaries. These connections resist changes of heights of the two rectangles, and can have a different elastic modulus along the hinge and the blade, denoted k H and k, respectively. Thus, force balance at the anterior and posterior boundaries in the hinge and in the blade can be written: Here, h 0 and h H 0 are the heights where the springs are relaxed in the blade and in the hinge. We take these reference heights to be equal to the initial blade and hinge heights.

Extracellular-matrix ablation at the anterior margin
Laser cut experiments severing the blade from the cuticle at the anterior margin modify the elastic modulus resisting the deformation of the blade along the y direction ( Other boundary conditions are not modified and are taken as in the WT case.

Extracellular-matrix ablation at the distal margin
The distally severed blade retracts distally away from the cuticle (Figure 8-figure supplement 3A and Figure 10-figure supplement 1A). As a result, the total length of hinge and blade L can change. The dynamics for the evolution of the length of the blade is given by the dynamic equation: where γ is an effective friction coefficient characterising dissipative processes external to the tissue. We use the same coefficient γ in Equations 72, 74, because we expect that both frictions arise from the same frictional forces acting between the wing and the surrounding cuticle.

dumpy ov1 mutants
In dumpy ov1 mutants the distal margin of the wing detaches from the cuticle similar to laser ablation of ECM on the distal margin ( Figure 8A,B and Figure 10-figure supplement 1A). We therefore use the boundary condition Equation 74 for dumpy ov1 mutants.