Self-organized patterning of cell morphology via mechanosensitive feedback

Tissue organization is often characterized by specific patterns of cell morphology. How such patterns emerge in developing tissues is a fundamental open question. Here, we investigate the emergence of tissue-scale patterns of cell shape and mechanical tissue stress in the Drosophila wing imaginal disc during larval development. Using quantitative analysis of the cellular dynamics, we reveal a pattern of radially oriented cell rearrangements that is coupled to the buildup of tangential cell elongation. Developing a laser ablation method, we map tissue stresses and extract key parameters of tissue mechanics. We present a continuum theory showing that this pattern of cell morphology and tissue stress can arise via self-organization of a mechanical feedback that couples cell polarity to active cell rearrangements. The predictions of this model are supported by knockdown of MyoVI, a component of mechanosensitive feedback. Our work reveals a mechanism for the emergence of cellular patterns in morphogenesis.


INTRODUCTION 31
During morphogenesis, tissues with complex morphologies are formed from the 32 collective interplay of many cells. This process involves spatial patterns of signaling activity, 33 and previous work has discovered mechanisms for generating tissue-scale patterns of 34 activity in signaling pathways such as Hedgehog, TGF , and Wnt (Green & Sharpe, 2015). 35 In addition, patterns of cellular morphology arise during morphogenesis. Such patterns can 36 be important for ensuring the function of the resulting tissue. For example, the compound eye 37 of Drosophila consists of hundreds of ommatidia organized in a precise hexagonal array that 38 is required to fully sample the visual field (Kumar, 2012). Patterns of cellular morphology that 39 arise during morphogenesis can also guide the morphogenetic processes itself. For example, 40 spatial patterns of cell morphology emerge during growth of the Drosophila larval imaginal 41 discs, which are precursors of adult tissues (Aegerter-Wilmsen et al., 2010;Condic et al., 42 1991;LeGoff et al., 2013;Mao et al., 2013). These patterns have been proposed to be 43 involved in the eversion process, during which these flattened epithelial sacs turn themselves 44 inside out when the animal transitions from larva to pupa (Condic et al., 1991). While 45 extensive work has studied the emergence of biochemical signaling patterns, how patterns of 46 cellular morphology arise during tissue development is poorly understood. 47 Here, we investigate tissue-scale patterning of cell morphology in the Drosophila larval 48 wing imaginal disc, which has a geometry that is ideal for studying spatial patterns of 49 epithelial cell morphology. We focus specifically on the cell shape patterns in the central 50 "pouch" region, which is the precursor of the adult wing blade. To a good approximation, this 51 region is planar and ellipsoidal. Cells near the center have smaller cell areas and are more 52 isotropic in shape, whereas cells near the periphery have larger cell areas and are elongated 53 tangentially (Aegerter-Wilmsen et al., 2010;LeGoff et al., 2013;Mao et al., 2013). Cell shape 54 has been correlated with mechanical stress: tangentially oriented bonds of elongated cells in 55 the periphery are under higher tension than radially oriented bonds (LeGoff et al., 2013). 56 It has been previously proposed that this pattern of cell morphology in the wing pouch 57 could stem from differential proliferation: if the center grows faster than the rest, the resulting 58 area pressure could stretch peripheral cells tangentially (Mao et al., 2013). Indeed, there is 59 evidence to suggest that cells divide slightly faster closer to the center during very early 60 stages (before 80hr after egg laying, AEL). It was suggested that this early growth differential 61 is sufficient to account for the persistence of the cell morphology pattern through the 62 remaining ~40hrs of development. However, it has since been shown that cell 63 rearrangements occur (Dye et al., 2017;Heller et al., 2016), which could relax stress patterns 64 once growth has become uniform. Furthermore, stress patterns may even relax during 65 homogeneous growth in the absence of cell rearrangements (Ranft et al., 2010). Thus, it 66 remains unclear how cell morphology patterns generated early by differential growth could be 67 maintained through later stages, and alternative mechanisms for the establishment of these 68 patterns must be considered. 69 Here, we measure the spatial patterns of cell morphology, cell divisions, and cell 70 rearrangements during the middle of the third larval instar (starting at 96hr AEL). We quantify 71 the pattern of tangential cell elongation and show that it becomes stronger over time, even 72 though growth is spatially uniform and cell rearrangements are frequent. Strikingly, this 73 change in tangential cell elongation is coupled to a radially biased pattern of cell neighbor 74 exchanges. Using a physical model of tissue dynamics, we show that active patterning of 75 radial cell neighbor exchanges can account for the observed morphology patterns in the 76 absence of differential growth. Lastly, using a combination of experiment and theory, we 77 provide evidence that this active patterning is self-organized by mechanosensitive feedback. 78

80
Cell morphology patterns can persist and strengthen in the absence 81 of differential growth 82 Cell morphology patterns in the wing disc have been previously analyzed using static 83 images (Aegerter-Wilmsen et al., 2010;LeGoff et al., 2013;Mao et al., 2013). However, 84 relating cell morphology patterns to patterns of growth, cell divisions and cell rearrangements 85 requires dynamic data. We therefore performed long-term timelapse imaging of growing 86 explanted wing discs using our previously described methods (Dye et al., 2017), starting at 87 96 ℎ AEL and continuing for ~13 ℎ of imaging. We used E-cadherin-GFP as an apical 88 junction marker ( Fig 1A). 89 To quantify morphology, we averaged apical cell area and cell elongation locally in 90 space, using data from 5 wing discs, and in time using ~2 ℎ intervals. We present the 91 spatial patterns calculated for the middle timepoint in Fig 1B-C. Cell area is represented as a 92 color code. Cell elongation is characterized by a tensor, which defines an axis and a strength 93 of elongation and is represented by bars in Fig 1C (see Methods, Fig S1B) (Etournay et al.,94 2015; Merkel et al., 2017). To quantify the radial symmetry of this pattern, we first determined 95 the center of symmetry (Fig S1 and methods). We then introduce a polar coordinate system 96 at the center with the radial coordinate and present the radial projection of cell 97 elongation tensor as a color code in Fig 1C (see also Fig S1). This figure highlights the 98 pattern of tangential cell elongation, with cells elongating on average perpendicular to the 99 radial axis (blue in Fig 1C). It also reveals that this pattern is interrupted around the Dorsal-100 Ventral (DV) boundary, where cells are elongated parallel to this boundary (red in Fig 1C). 101 We quantify this region separately ( Fig S2) and exclude it from our analysis of the circular 102 patterns (Fig 1D-F). The Anterior-Posterior (AP) boundary also affects cell morphology 103 (Landsberg et al., 2009), but the effect is weaker and more variable at this stage; thus, we do 104 not quantify it separately. 105 The spatial maps of cell morphology reveal that both cell area and cell elongation 106 magnitude are largest at the periphery and decrease toward the center (Fig 1B-C). We 107 quantified this radial gradient in cell area and observe that cell area ranges from ~3 − 7 108 when moving toward the periphery (Fig 1E). We also observe a radial gradient in cell 109 elongation starting from ≈ 0 at = 10 and extending to about ≈ −0.1 at = 110 40 ( Fig 1F). The negative value corresponds to tangential elongation ( Fig 1C). When 111 evaluated over the timelapse, we find that these radial gradients grow slightly more 112 pronounced over time (Fig S3C-D). 113 As previously proposed, differential growth can generate such patterns of cell 114 elongation (Mao et al., 2013). However, indirect metrics of tissue growth in vivo do not 115 indicate that differential growth still occurs at this later stage (Mao et al., 2013). We directly 116 measured the spatial pattern of growth during timelapse. Cell division rate has been 117 previously used as an indicator of growth; however, tissue growth actually results from a 118 combination of cell division, cell area changes, and cell extrusions ( Fig 1G). Thus, we 119 evaluated the spatial pattern (Fig 1H, I, K, L) and radial profiles (Fig 1J, M) of total tissue 120 growth and its cellular contributions. These data show that tissue area growth, as well as cell 121 division rate, are to a good approximation independent of the distance to the center. 122 In summary, we quantified cell morphology patterns in mid-third instar wing explants 123 during live imaging. We confirm previous static observations of the pattern and further 124 identify a region around the DV boundary with a morphological pattern that differs from the 125 rest of the wing pouch. We quantify the radial gradients in cell area and cell elongation 126 existing outside of the DV boundary region and show that they strengthen in time in the 127 absence of differential growth, raising the question of what mechanism underlies the 128 persistence of these cell morphology patterns during mid to late stages of wing growth. 129 130 Radially-oriented cell rearrangements balance tangential cell 131 elongation 132 To directly relate the observed cell morphology patterns with cell rearrangements, we 133 next analyzed the spatial patterns in cellular contributions to tissue shear (Fig 2A). Tissue 134 shear can be decomposed into contributions from cell divisions, cell elongation changes, T1 135 transitions, and so-called correlation effects (Etournay et al., 2015;Merkel et al., 2017). Here, 136 correlation effects result mainly from correlated fluctuations in cell elongation and cell rotation 137 (see Theory Supplement part III and (Merkel et al., 2017)). 138 We find that the spatial patterns of tissue growth and its cellular contributions exhibit 139 overall anisotropies perpendicular to the DV boundary, as reported previously (Dye et al., 140 2017). In addition, the patterns of cell elongation changes and T1 transitions can be 141 described as a superposition of a uniformly oriented pattern and an approximately radial or 142 tangential pattern (Fig 2C-D). To determine the magnitude of the radial or tangential patterns 143 in all quantities, we quantified their average radial projections as cumulative plots over time 144 ( Fig 2B). The radial component of tissue shear is small, and cell divisions do not contribute to 145 radial tissue shear. In contrast, we observe a pronounced buildup of a tangential pattern of 146 cell elongation accompanied by a radial pattern of T1 transitions and of correlation effects. 147 As shown above (Fig 1J), tissue area growth does not have a radial gradient and thus 148 does not contribute to the increase in tangential cell elongation that we observe at this time 149 (Fig 2A-B, Fig S3C-D). Furthermore, we observe numerous T1 transitions (on average 150 1.0 ℎ ), and their radially-biased orientation increases rather than relaxes tangential 151 cell elongation (Fig 2E-F). Thus, we are not observing the relaxation of a pattern of cell 152 elongation caused by early differential growth. Rather, our data support a model whereby a 153 radially patterned morphogenetic cue actively biases the direction of T1 transitions and 154 consequently the complementary pattern of cell shape changes. 155 156 Polarity-driven cell rearrangements can create the observed cell 157 morphology patterns in the wing disc 158 We next apply a biophysical model to determine whether radially patterned T1 159 transitions could account for the observed cell morphology patterns in the wing disc. We 160 propose that the wing disc has a patterning cue that biases T1 transition orientation. This cue 161 could be represented as a cell polarity that is nematic in nature, having a magnitude and an 162 orientation axis, which is the same symmetry as that of an elongated cell. To capture the 163 mechanics and dynamics of the epithelium, we adapt a physical model of the tissue material 164 properties that takes into account the interplay of T1 transitions, cell shape changes, and 165 tissue shear in a continuum description (Etournay et al., 2015;Popović et al., 2016). Such a 166 model can predict the patterns of cell morphology that arise from patterned T1 transitions. 167 In our model, we consider the spatial patterns of tissue shear rate , the patterns of 168 cell shape , and the patterns of cell rearrangements , which have nematic symmetry and 169 are represented by traceless symmetric tensors that describe the magnitude and orientation. 170 Tissue shear results from a combination of cell shape changes and rearrangements 171 according to: 172 = + (1),

173
where D /D is a co-rotational time derivative of , and the shear due to rearrangements 174 includes contributions from T1 transitions, cell divisions and extrusions, and also correlation 175 effects. Eq 1 follows from geometric considerations (Etournay et al., 2015;Merkel et al., 176 2017).

181
Here, the shear stress tensor is the sum of the elastic stress associated with cell 182 elongation and the active stress associated with the cell polarity cue (Fig 3A). The elastic 183 stress is characterized by the shear elastic modulus . The cell polarity cue is provided by 184 the nematic tensor , and the coefficient describes the anisotropic stress associated with 185 the polarity cue. The shear due to cell rearrangements given by Eq 3 is driven in part by 186 shear stress, and therefore depends on the cell elongation , and in part by active processes 187 that are oriented by the nematic cell polarity (Fig 3B). is a relaxation time for cell 188 rearrangements, and is the rate of cell rearrangements driven by polarity. 189 To discuss the wing disc, we consider a radially symmetric geometry and average the 190 oriented quantities after projection onto the radial axis. Radial tissue shear is small compared 191 to that associated with cell shape changes and T1 transitions during our observed time 192 window ( Fig 2B). We therefore consider a steady state with = 0, D /D = 0, and = 193 0. In this case, cell elongation becomes: 194 = − (4).

195
Thus, we find that the steady state cell elongation pattern is a result of cell rearrangements 196 that are oriented by the cell polarity cue (Fig 3C-D). Note that our data show that the wing 197 disc is not exactly at steady state: cells slowly change their shape and rearrange radially ( Fig  198  2A-B). However, as we show in Theory Supplement part I, Eq. 4 holds to a good 199 approximation. 200 Can the radial pattern of T1 transitions defined by also explain the observed radial 201 profile of cell area (Fig 1B, E)? To answer this question, we then considered force balances 202 in the tissue. We consider tissue area pressure 203

204
where is the difference in pressure from a reference value, is tissue area 205 compressibility, is the average cell area, and is a reference cell area. As pressure 206 increases, cell area decreases. To calculate the cell area profile, we again approximate the 207 wing pouch as a radially symmetric disc. In the radially symmetric geometry, force balance 208 can be expressed as: 209 = + (6).

210
A radial profile of pressure determined from this equation implies a radial pattern of cell area 211 via Eq 5 (see Fig 3C-E and Theory Supplement part I). To fit this equation to our data, we 212 first estimate the profile of by measuring the radial profile of cell elongation using Eq 4 213 ( Fig 3F). Then, we solve for the cell area profile using Eqs 2, 5, and 6 ( Fig 3G and Theory 214 Supplement part I). We obtain a fit that accounts for the observed pattern of cell area. 215 From this analysis, we conclude that the cell morphology patterns observed in the 216 wing disc could be generated by radially-biased cell rearrangements. Next, we test whether 217 the stress profile predicted by the model (Eq 2) exists in the tissue, and we measure key 218 mechanical parameters of the model. Later, we address the potential molecular origin of the 219 cell polarity cue orienting the cell rearrangements. 220 221 Circular laser ablation reveals patterns of tissue stress 222 Our model predicts a stress pattern in the wing disc that results from active processes 223 that are radially oriented by a cell polarity cue. To compare this prediction to experiment we 224 infer tissue stress using laser ablation. Tissue stress has been estimated previously by laser 225 ablation techniques that are based on determining the initial retraction velocity (Bonnet et al., 226 2012;Etournay et al., 2015;Farhadifar et al., 2007;LeGoff et al., 2013;Mao et al., 2013;227 Shivakumar & Lenne, 2016). However, to compare theory and experiment, ideally one 228 should measure quantities that are well-captured by the model. Therefore, instead of using 229 initial retraction velocity, we perform circular cuts and analyze the final, relaxed position of 230 the inner and outer elliptical contours of tissue formed by the cut (Fig 1A, Theory Supplement 231 part II). From the size and the anisotropy of the cut, we can infer anisotropic and isotropic 232 tissues stress, normalized by the respective elastic constants, as well as the ratio of elastic 233 constants (see Theory Supplement part II). Furthermore, we can also infer the existence of 234 polarity-driven stress. We name this method ESCA (Elliptical Shape after Circular Ablation). 235 We perform local measurements of tissue mechanics by cutting the tissue in the 236 smallest possible circle that would still allow us to measure the shape of the inner piece left 237 by the cut (radius = 7 , encircling ~ 5 − 15 cells, Fig 4A). To relate measurements of 238 tissue stress to cell elongation, we calculate the average cell elongation in the ablated region 239 before it is cut ( Fig 4A). In Fig 4B, we present the measured cell elongation and shear stress 240 tensors at the position of ablation. We find that local average cell elongation correlates well 241 with the direction of shear stress. Also, we observe that the cells in the band around the DV 242 boundary have different mechanical properties than elsewhere in the tissue. Near the DV 243 boundary, cells elongate less than outside this region for comparable amounts of stress ( Fig  244  4B and Fig S4B, E). The ratio of elastic constants in this region is also smaller: near the DV 245 boundary 2 / = 2.3 ± 0.3 , whereas outside this region, 2 / = 3.4 ± 0.4 (see also Fig  246  S4C, F). We focus hereafter on the radial patterns of elongation and stress outside of the DV 247 boundary region. 248 The relationship between cell elongation and stress normalized by the elastic 249 modulus has a slope 1 in the absence of polarity-driven stress (see Eq 2). We observe a 250 much smaller slope for this relationship in our data (Fig 4C), indicating that polarity-driven 251 stress is significant. We now use these data to estimate the parameters of our mechanical 252 model. We write the shear stress defined in Eq 2 in terms of cell elongation and cell 253 rearrangements, eliminating the orientational cue using Eq 3. For the radial components, 254 we have: 255 estimate for the tissue relaxation time = 2 ± 2 ℎ , which is roughly consistent with that 261 found during pupal morphogenesis (Etournay et al., 2015). From our data, we can also infer 262 the radial profile of tissue area pressure, revealing that pressure increases towards the 263 center ( Fig 4D and Theory Supplement). This finding is consistent with the observed cell 264 area profile, with smaller cell areas towards the center (Fig 1B,E). 265 In sum, we find a stress profile in the wing disc that is consistent with the observed 266 measurements of both cell elongation and area. Further, we use these data to measure 267 certain parameters of our biophysical model, including the tissue relaxation timescale and the 268 effective shear elastic coefficient. In our model, the radial orientation cue is required to generate the observed patterns of 273 cell morphology, cell rearrangement, and tissue stress. Candidates for such an orientational 274 cue are the planar cell polarity pathways (PCP), which are groups of interacting proteins that 275 polarize within the plane of the epithelium. There are two well-characterized PCP pathways: 276 Fat and Core (Butler & Wallingford, 2017;Eaton, 2003). In the wing, these systems form 277 tissue-scale polarity patterns during growth (Brittle et al., 2012;Merkel et al., 2014;Sagner et 278 al., 2012) and are required to position the hairs and cuticle ridges on the adult wing (Adler et 279 al., 1998;Doyle et al., 2008;Eaton, 2003;Gubb & Garcia-Bellido, 1982;Hogan et al., 2011). 280 To determine whether either of these pathways could function as the orientational cue 281 described in our model, we analyzed cell elongation patterns after their removal. 282 We perturbed the Fat pathway using Nub-Gal4 to drive the expression of RNAi 283 constructs targeting both Fat and Dachs (D) in the pouch region throughout the third larval 284 instar. Consistent with previous work showing that loss of both Dachs and Fat suppresses 285 wing growth (Cho & Irvine, 2004), we observe a significantly smaller wing pouch upon Fat/D 286 double knockdown ( Fig 5A). Nonetheless, the pattern of tangential cell elongation persists to 287 the end of larval development (Fig 5A-C). Using scaled coordinates, we find that the radial 288 profiles of cell elongation in Fat/D RNAi and control wings are similar ( Fig 5D). 289 We perturbed the Core PCP pathway using a previously characterized null mutation in 290 Prickle (pk 30 ), which causes defects in adult wing hair orientation (Gubb et al., 1999). We 291 found that the cell elongation pattern in the Pk mutant is similar to the wild type control ( Fig  292  5E-I). In the Pk mutant, the region of tangential cell elongation extends even further into the 293 center than in control wings. 294 We conclude that the tangential cell elongation pattern persists in the absence of 295 either PCP pathway. This result excludes these pathways as orienting cues for the cell 296 elongation patterns. 297

298
Mechanosensitive feedback generates self-organized patterns of 299 cell morphology 300 We have shown that perturbing PCP pathways does not affect the radial patterns of 301 morphology, raising the question of how orientational cues might arise. In previous sections, 302 we have considered the orientational cue to be provided by a cell polarity system that is 303 independently patterned. However, cell polarity in general would be affected by stresses in mechanics can give rise to spontaneous emergence of the cell polarity cue (Fig. 6). 308 Mechanosensitivity is incorporated into our model of tissue mechanics through a 309 dynamic equation for the orientational cue that becomes stress-dependent: 310

311
Here, is a relaxation timescale for , is a mechanosensitive feedback strength, the 312 coefficient > 0 ensures stability, and is a coupling strength locally aligning orientational 313 cues. 314 Now, Eq 8 provides a mechanosensitive feedback to Eqs 1, 2 and 3. These combined 315 equations show a novel behavior. In particular, the orientational cue can emerge 316 spontaneously by self-organization ( Fig 6A-B). Beyond a critical value of the 317 mechanosensitive feedback strength , an isotropic tissue with = 0 is no longer stable, and 318 a state with an orientational cue ≠ 0 emerges instead ( where a positive coefficient is needed to stabilize the polarized state. By this 321 mechanism, the anisotropic cue introduced earlier in our model can be locally generated by 322 mechanosensitive self-organization and does not require the existence of pre-patterned 323 polarity cues. To generate a large-scale pattern from locally generated anisotropic cues, they 324 need to be aligned in neighboring regions. This local alignment is captured in Eq 8 by the 325 orientation coupling term with strength , which is similar to alignment terms found in 326 anisotropic physical systems, such as liquid crystals (Gennes & Prost, 1993;Jülicher et al., 327 2018; Marchetti et al., 2013). 328 To discuss cell morphology profiles in the wing disc, we consider a simplified tissue 329 model with radial symmetry, where the rate of radial cell rearrangement is given (as 330 estimated in Theory Supplement part IB) and the cell shape pattern and tissue stress pattern 331 are calculated. Using a fit of cell elongation to the experimental data, we find a set of 332 parameter values that accounts for the observed cell elongation patterns in the wing disc 333 (Table 1, Theory Supplement). From this cell elongation pattern also follows the cell area 334 pattern (as described above, Fig 3). 335

Fit parameter
Parameter

341
We conclude that our mechanosensitive model can account for the radial pattern of 342 cell morphology in the wing disc. Note that because of the relatively large number of 343 parameters used to fit a single experimental curve, there are large uncertainties when 344 estimating parameter values. Nonetheless, the qualitative prediction that reduction in 345 mechanosensitivity would lead to less polarization and thereby reduced cell elongation is 346 rather robust and thus not likely to be affected by these uncertainties (see Theory 347 Supplement part ID). We next test this prediction of our model experimentally. In order to test our prediction that the reduction of mechanosensitivity will reduce the 352 magnitude of cell elongation, we used RNAi to knockdown Myosin VI, a molecular motor 353 implicated in mechanosignaling. Myosin VI is an upstream component of a Rho-dependent 354 signaling pathway that reorganizes the actin-myosin cytoskeleton in response to mechanical 355 stress (Acharya et al., 2018). Experiments in wing discs also indicate that mechanosensation 356 involves Rho polarization and signaling (Duda et al., 2019). We performed RNAi against 357 MyoVI in the wing pouch using Nub-Gal4 and evaluated cell morphology at the end of larval 358 development (~120 ℎ AEL). We observe a clear reduction in the magnitude of tangential cell 359 elongation as compared to wild type at this stage ( Fig. 7A-D, Fig S5). In addition, our model 360 predicts that such reduction of cell elongation would result in an increase of cell area in the 361 central parts of the wing (Eqs 1-6). The observed pattern of increased cell area in MyoVI 362 knockdown experiments is consistent with this prediction (Fig 7E-F). Therefore, the 363 qualitative predictions of our model upon reducing the mechanosensitive feedback strength 364 are confirmed by the experimental downregulation of the mechanosensitive motor, MyoVI. 365

367
Here, we have shown that patterns of cell shape and stress in the mid-third 368 instar Drosophila wing disc do not rely on PCP pathways or differential growth. Instead, 369 radially-oriented T1 transitions and tangential cell elongation emerge via mechanosensitive 370 feedback in a self-organized process. We have presented a continuum model of tissue 371 dynamics for this self-organization based on a mechanosensitive nematic cell polarity that 372 accounts for the observed patterns of cell area, T1 transitions, and cell shape. Our work 373 highlights a mechanism for the self-organized emergence of cellular patterns in 374 morphogenesis. Our work shows that the spatial pattern of T1 transitions is an integral part of the 379 emergence of tissue organization during wing development. In contrast to situations such as 380 germband extension, where T1 transitions exhibit clearly discernible patterns, the patterns of 381 T1 transitions in the wing disc have been elusive. Many T1 transitions occur in the tissue in 382 seemingly random orientations. However, on average, they exhibit a spatial pattern. We 383 revealed these patterns by quantifying the nematics of T1 transitions and cell shape changes 384 using the previously-described triangle method (Merkel et al., 2017) and then quantified them 385 with radial averaging (Fig 2A-B). In this way, we revealed that a radial pattern of T1 386 transitions is linked to a tangential pattern of cell elongation. 387 Given this radial pattern of T1 cell rearrangements, the observed cell morphology 388 pattern follows from a continuum tissue model based on a radially-oriented nematic cell 389 polarity field (Fig 3). The polarity-oriented radial T1s create a cell shape pattern with 390 corresponding patterns of tissue stress and tissue area pressure. The 2D area pressure is 391 higher in the center and is lower towards the periphery. Note that this pressure profile does 392 not rely on differential proliferation, as was previously proposed (Mao et al., 2013) but 393 instead relies on a radial pattern of T1 transitions. We test our model using a novel circular 394 laser ablation method. This method allows us to determine specific combinations of tissue 395 parameters. In particular, we estimate the ratio of elastic constants 2 / = 3.4 ± 0.4 396 and * / = 0.05 ± 0.02, as well as the cell shape relaxation timescale = 2 ± 2ℎ . 397 This analysis raised the question of which nematic cell polarity cues guide the cell 398 rearrangement and cell elongation patterns. PCP pathways are required for the proper 399 orientation of T1 transitions in other contexts . However, we found that 400 neither of the two known PCP pathways in the wing are required for the observed tangential 401 cell elongation patterns (Fig 5). We instead show that an orientation cue can arise through 402 self-organization via mechanosensitivity and identify MyoVI as a key molecular player.

404
Mechanosensitive feedback can create self-organized patterns of 405 cell morphology 406 Cell polarity cues can emerge via mechanosensitive feedback by transforming 407 mechanical cues into chemical anisotropies. Nematic cell polarity can then orient active 408 stresses and thereby amplify the mechanical stimulus (Fig 6). We introduce this 409 mechanosensitive feedback in our continuum theory, which quantitatively describes the 410 emergence of patterns of cell shape and cell rearrangements. The strength of this 411 mechanosensitive feedback is described by a parameter . If exceeds a critical value c, an 412 orientation cue and elongated cell shapes spontaneously emerge, while for < c, the tissue 413 remains isotropic (Fig 6B). This model can account for the observed patterns of cell area and 414 cell elongation in the wing disc and predicts that the reduction of mechanosensivity will result 415 in reduced cell elongation. To test this prediction, we perturbed a RhoA-dependent 416 mechanotransduction pathway by knocking down an upstream component, MyoVI, using 417 RNAi (Fig 7). We find a clear phenotype of reduced cell elongation and increased cell areas 418 in the center region, as predicted by our model. 419 This result, together with the fact that MyoVI is involved in Rho-dependent activation of 420 actin-myosin cytoskeleton (Acharya et al., 2018), suggests that MyoVI is a molecular 421 component of the mechanosensitive feedback we describe in our self-organized model. 422 However, the molecular nature of the cue that defines the nematic cell polarity is unknown. 423 This cue may organize the structure or dynamics of the actin-myosin cytoskeleton or the 424 actin-myosin system itself could define nematic cell polarity. Indeed, it has been shown that 425 myosin II localizes to long cell boundaries in the wing (LeGoff et al., 2013), corresponding to 426 a nematic polarity aligned with the nematic cell polarity . Also, wing disc stretching 427 experiments have shown that MyoII can polarize in response to exogenous stress in a Rho-428 dependent manner (Duda et al., 2019). Precisely how the actin-myosin cytoskeleton is 429 affected by MyoVI in this system and how these cytoskeletal elements together guide cell 430 rearrangements in response to anisotropic tissue stresses and cell shape changes remain 431 open questions for future research. 432 Lastly, our laser ablation analysis shows that the region around the DV boundary has a 433 different ratio of elastic constants than the rest of the tissue, which could affect the self-434 organized pattern formation we describe. Therefore, it will be interesting to study how 435 Wingless/Notch signaling, which defines the DV boundary, may influence the mechanical 436 properties that lead to mechanosensitive self-organization of polarity and morphology. In 437 addition, we observe a richer pattern emerging very late in development (see Fig S5, Fig 5), 438 including a region anterior to the AP boundary that is radially elongated. Future research will 439 expand upon the model presented here to explore the dynamics of these patterns. 440 In summary, we used the Drosophila wing disc to identify a mechanism by which tissue 441 morphology can arise from the self-organization of a mechanical feedback coupling cell 442 polarity to active cell rearrangements. This mechanism is general and could be employed in 443 other tissues and organisms to generate patterns of cell shape and cell area.  Cumulative tissue deformation and its cellular contributions, measured in a grid centered on 502 the compartment boundaries (dotted lines) and averaged over all 5 movies. Bars represent 503 nematic tensors, where the length is proportional to the magnitude of deformation and the 504 angle indicates its orientation. The contribution from cell extrusion is small and thus not 505 shown. (B) The radial projection of cumulative tissue deformation and its cellular 506 contributions are plotted as a function of time after the first 2 ℎ of adaption to culture (Dye et 507 al., 2017). Solid lines indicate the average over all five movies; shading indicates the 508 standard deviation. (C) Schematics indicating how a uniformly oriented pattern would 509 combine with a tangential pattern to produce a pattern resembling that of the cell elongation 510 change (left) or with a radial pattern to produce a pattern resembling that of the T1 transitions 511 (right). (E-F) Illustrations demonstrating radially-oriented T1 transitions coupled to tangential 512 cell elongation. For simplicity, we diagram only the posterior-dorsal quadrant, but the pattern 513 is radially symmetric. In (E), we define the radial coordinate system, where velocity in the 514 radial direction is positive and that in the negative direction is negative. Patterns in A indicate 515 that T1 transitions are biased to grow new bonds in the radial direction, and cell elongation is 516 biased to increase tangentially. (F) In a radially oriented T1 transition, cells preferentially 517 shrink tangentially oriented bonds and grow new bonds in the radial direction. When oriented 518 in this direction, T1 transitions do not dissipate tangential elongation but increase it (green 519 bars in each cell represent cell elongation). 520 elongation is averaged in this region before the cut. After the tissue relaxes (2 ), we fit 537 ellipses to the cut region. We then infer mechanical properties using our model, which inputs 538 deformation and outputs stress as a function of elastic constant and the ratio of the isotropic 539 and anisotropic elastic constants (see also Fig S4). genotype was quantified by averaging in radial bins and plotting as a function of . Since 565 the fat/d RNAi wings are smaller, we present this profile as a relative distance to the center. 566

(B) Map of stress and cell elongation in
The band of cells around the DV boundary was removed, and because the fat/d RNAi discs 567 are smaller, this region of exclusion is a larger relative distance. Plots in (C-D) represent 568 averages of = 7 − 8 wing discs per genotype; for plots in (G-I), = 11 − 14 per genotype. 569 the orientation coupling term . We find a set of parameters that can account for the 576 experimental data on cell elongation in the wing disc (See Theory Supplement part ID and 577 Table 1). 578 analyzed, and the sex of the animals was not recorded, as we have no reason to believe 601 there is any sexual dimorphism in the studied phenomenon. To synchronize development, 602 we collected eggs deposited within a defined time window on apple juice agar plates. To do 603 so, we transferred the flies from standard food vials to cages covered by apple juice agar 604 plates containing a dollop of yeast paste for food. After at least 2 ℎ , the plates were 605 replaced, and the timing of collection started. Eggs laid within a ~2 ℎ time window were 606 collected by cutting out a piece of the agar and transferring it to a standard food vial. We 607 limited the number of eggs per vial to < 15 to avoid crowding. The middle of the time window 608 for egg collection was considered to be 0 ℎ after egg laying (AEL). Experiments from 609 timelapse imaging and laser ablation (Figs 1, 2 Wing explants were grown ex vivo as described previously (Dye et al., 2017). Briefly, 619 wing discs were dissected from larvae in growth media (Grace's cell culture media + 5% fetal 620 bovine serum + 20nM 20-hydroxyecdysone + Penicillin-Streptomycin) at room temperature. 621 Then, they were transferred to a Mattek #1 glass bottom petri-dish to the center of a hole cut 622 in a double-sided tape spacer (Tesa 5338) and covered with a porous filter. The dish was 623 then filled with fresh growth media. 624 Acquisition:

625
Data from movies 1-3 were used previously (Dye et al., 2017). Movie 4 was acquired 626 after publication of the first manuscript, in the same exact way as movies 1-3. Briefly, E-627 cadherin-GFP-expressing wing discs were imaged in growth media using a Zeiss spinning-628 disc microscope to acquire 0.5 spaced Z-stacks at 5 intervals with a Zeiss C-629 Apochromat 63X/1.2NA water immersion objective and 2x2 tiling (10% overlap). This 630 microscope consisted of an AxioObserver inverted stand, motorized xyz stage, stage-top 631 incubator with temperature control set to 25C, a Yokogawa CSU-X1 scanhead, a Zeiss 632 AxioCam MRm Monochrome CCD camera (set with 2x2 binning), and 488 laser illumination. 633 We circulated growth media during imaging using a PHD Ultra pump (Harvard Apparatus) at 634 a rate of 0.03 / . Time 0 ℎ of the movie was considered to be the start of imaging 635 acquisition, typically 45 − 60 from the start of dissection (time required for sample and 636 microscope preparation). 637 Movie 5  Processing: 648 Raw Z-stacks were denoised using a frequency bandpass filter and background 649 subtraction tools available in FIJI (Schindelin et al., 2012). Then, we used a custom algorithm 650 as described previously (Dye et al., 2017) to make 2D projections of the apical surface, 651 marked by Ecadherin-GFP. This algorithm also outputs a height-map image, in which the 652 value for each pixel corresponds to the level in the Z-stack of the identified apical surface.

653
For those movies that were tiled, we used the Grid/Collection Stitching FIJI plugin (Preibisch 654 et al., 2009) to stitch the 2D projections, and then used that calculated transformation to 655 stitch the height map images, so that we could correct the cell area and elongation values for 656 local curvature (see below). We focus in this work exclusively on the disc proper layer, which 657 is the proliferating layer of the disc that goes on to produce the adult wing. We did not 658 consider cell shape patterns in the peripodial layer, the non-proliferating layer of the wing 659 disc that is largely destroyed at the onset of pupariation. Acquisition:

666
All imaging of single timepoint data was performed with the same microscope described 667 above for Movie 5. We did not acquire timelapse data for these genetic perturbations, and 668 thus we chose to image the entire disc using 2x2 tiling, 0.5 spaced Z-stacks. 669 Approximately 6-10 discs were imaged per dish, within approximately 30 − 60 . 670 Processing:

671
Tiled images were stitched together as Z-stacks; then we obtained the apical surface 672 projection and its corresponding height map as described above and in (Dye et al., 2017). 673 Circular Laser ablation 674 Acquisition:

675
Wing discs from mid-third-instar larvae (96 ℎ AEL) were dissected and mounted exactly 676 as was done for the live imaging timelapses. Due to constraints on the speed of ablation, we 677 only cut regions of the wing disc pouch that were completely flat, so that we could cut in a 678 single plane at the apical surface (rather than having to cut in each plane of a Z stack). 679 Where this flat region lies depends on how the disc happens to fall on the coverslip during 680 mounting. Because the anterior compartment is larger and higher, most of our ablations are 681 in this compartment. To access other regions, we also mounted some wing discs on an 682 agarose shelf: stripes of 1% agarose (in water) were dried onto the surface of the coverslip 683 and wing discs were arranged their anterior half propped up on the agarose shelf prior to 684 adding the porous filter cover. 685 Ablations were performed using ultraviolet laser microdissection as described in (Grill et 686 al., 2001) using a Zeiss 63X water objective. First, we took a full Z-stack of the sample prior 687 to the cut. Then, we selected a 7 radius circle that would ablate in the flattest region of 688 the tissue. No imaging is possible during ablation, but we acquired a 2 timelapse 689 immediately after the cut in a single Z-plane. This timelapse data was not used except to 690 estimate whether or not the sample was fully cut. After, once the sample has finished 691 expanding but not started to heal (~2 ), we took another full Z-stack to image the 692 endpoint. We excluded a small number of data points if any of the following were true: (1) the 693 inner piece remaining after the ablation was no longer visible (sometimes it floats or is 694 destroyed); (2) the cut appeared to expand highly asymmetrically (rare); (3) the wing discs 695 were clearly too young to be considered 96 ℎ AEL (poor staging). 696 Immunofluorescence after ablation:

697
Due to a limited field of view on the microscope used for ablation, we performed 698 immunofluorescence after the ablation in order to better estimate the position of the cut in the 699 wing pouch. After all the discs in the dish were ablated (<10 discs/dish), the entire dish was 700 fixed through the filter by adding 4% PFA and incubating 20min at room temperature. After, 701 the dish was rinsed and kept in PBST (PBS + 0.5% Triton X-100) until all discs from that 702 image acquisition day were completed (2 − 4 ℎ ). All samples were then blocked using 1% 703 BSA in PBST+250mM NaCl for 45 − 1 ℎ , and then incubated in primary antibody 704 overnight (diluted in BBX: 1% BSA in PBST). Initially, we labeled samples with SRF (Active 705 Motif, 1:100 dilution), but later we switched to Patched and Wingless (DSHB, 1:100) to 706 identify the compartment boundaries. Primary antibody toward GFP (recognizing the E-707 cadherin-GFP) was also included. After overnight incubation at 4C, we washed with BBX, 708 followed by BBX+ 4% normal goat serum (NGS), for at least 1 ℎ . Secondary antibodies 709 were added for 2 ℎ at RT in BBX+NGS. Finally, samples were washed 4x10min in PBST 710 and imaged in this media. Imaging was performed on one of the two spinning disc 711 microscopes described above for live imaging. We matched the stained samples with the 712 ablation images by (1) keeping track of the position of each disc on the dish and where it was 713 ablated during acquisition, and (2) morphology of the disc before/after ablation. All antibody 714 information is listed in the Key Resources Using the 2D projections, we performed cell segmentation and tracking using the FIJI 720 plugin, TissueAnalyzer . We manually corrected errors in the automated 721 segmentation and tracking as much as possible and then generated a relational database 722 using TissueMiner . Images were rotated to a common orientation 723 (Anterior left, dorsal up). We then manually identified three regions of interest at the last point 724 in the timelapse, using the FIJI macros included with TissueMiner: the "blade" was defined 725 roughly as an elliptical region surrounded by the most distal folds; and the Anterior-Posterior 726 and Dorsal-Ventral boundaries were estimated using the brightness of E-cadherin-GFP and 727 apical cell size (Jaiswal et al., 2006). 728 For the timelapse data, we used only the cells within these manually defined regions 729 that were trackable during the entire course of the movies. An example of this region is 730 presented for Movie 1 in Fig 1A. Furthermore, we also excluded from all of our analysis the 731 first 2 ℎ of imaging, the so-called adaption phase, where cells uniformly shrink in response 732 to culture (Dye et al., 2017). 733 For the RNAi data (Fig 5 and 7), we did not acquire timelapses, rather a single timepoint 734 at late stages of development. We nevertheless chose to analyze these data using the same 735 TissueMiner workflow for simplification. Because TissueMiner was developed for timelapse 736 data, however, it requires at least two timepoints. Thus, we duplicated the data and labelled 737 the two images as if they were timepoints 0 and 1 of a timelapse. We then only analyzed 738 timepoint 0. The regions of interest in these data, therefore, are manually defined and not 739 simply the region that is trackable (since it is static data). 740 Alignment on compartment boundaries:

741
To average across all movies of timelapse data, or all discs of the same genotype of the 742 static RNAi data, we generated a disc-coordinate system by normalizing the position of each 743 cell to the AP and DV boundaries for that disc. To do so, we averaged the absolute xy 744 positions of all the cells in the AP or DV boundaries over all time after the adaption phase. 745 We then calculated the distance of each cell to the new X axis (average position of the AP 746 boundary) and Y axis (average position of the DV boundary). 747 Analysis of cell size and shape 748 Definition of the cell elongation tensor:

749
The cell elongation tensor is a traceless symmetric tensor that quantifies the 750 anisotropy of cell shapes in a region of the tissue. We define cell elongation using a 751 triangulation of the tissue obtained by connecting centroids of connected cells (Etournay et 752 al., 2015;Merkel et al., 2017). The state of each triangle is described by a tensor , defined 753 by a mapping an equilateral reference triangle to the triangles in the tissue. The state tensor 754 contains information about the triangle elongation tensor , orientation angle , and area a: 755 Here, is the area of the reference triangle, and is 2-dimensional rotation matrix. Cell 757 elongation in the tissue region is defined as an area weighted average of triangle elongation. 758 For details about the method see (Merkel et al., 2017). 759 Adjustment of cell area and elongation to account for tissue curvature:

760
The wing disc pouch has a slightly domed shape (Fig S1A). After projection onto a 2D 761 surface, the cell shapes and areas will be distorted. To ensure that the radial profiles we 762 measure in a 2D projection (Fig. 1B,C,E, and F and Fig. 3F,G) are not a result of tissue 763 curvature, we account for the distortion caused by projection. We use the height maps 764 generated by our projection algorithm, which identifies the apical surface of the pouch within 765 the 3D Z-stacks. We smooth this height map with a gaussian filter of width = 2 to find 766 the height field ℎ( , ), and then we calculate the height gradient field. We then smooth the 767 result again with the same gaussian filter to find ⃗ ℎ = ( ℎ, ℎ). The deprojected cell or 768 triangle area is obtained from the measured area as = 1 + | ⃗ ℎ| , where ⃗ ℎ is 769 evaluated at the cell center. To find the deprojected cell elongation of a triangle we 770 evaluate ⃗ ℎ at the triangle center. Then, we find the angle of steepest ascent = 771 arctan( ℎ/ ℎ) in the projected plane and define the tilt transformation: We apply this transformation to the triangle shape tensor , as defined in (Merkel et al.,774 2017), to determine the deprojected triangle shape tensor 775 = 776 The cell elongation tensor is obtained from the triangle shape tensor as the 777 corresponding area weighted average using the deprojected cell area, as described in 778 (Merkel et al., 2017). 779 Spatial maps of cell size and shape:

780
To generate the color-coded smoothed plots of area and elongation (Fig 1B-C), we 781 divided the aligned wings into a grid with boxsize = 10pixels (~2 ). For each position, we 782 averaged cell area or performed an area-weighted averaged of triangle elongation in a 783 neighborhood box = 20 pixels (~4 ). To create the nematic elongation pattern, we similarly 784 averaged elongation of triangles whose centers lie within each grid box, with grid box size = 785 30 pixels (~6 ). 786 Calculation of radial elongation center point:

787
We define the center of the wing pouch to be on the DV boundary. To determine the 788 location of the center along the DV boundary we divide the pouch into four regions defined 789 by the DV boundary and a line perpendicular to it located at some position , as shown in Fig  790  S1E. Then, we define a function 791 for all values of along the DV boundary. Here, are the areas of the four regions 793 and ⟨ ⟩ are the area weighted averages of the cell elongation component 794 calculated in the four regions. Finally, is the value that minimizes ( ) (Fig S1F). We find 795 that the center point lies just anterior to the intersection with the anterior-posterior (AP) 796 boundary (Fig 1, S1, consistent with (LeGoff et al., 2013)). 797 In time-lapse experiments, cell centers were determined using the last 20 timepoints. 798 In single timepoint image experiments (Figs 5 and 7), all images of a single genotype were 799 used together after alignment on the DV and AP boundaries. 800 Exclusion of the band of cells near the DV boundary:

801
In contrast to the other regions (blade, AP and DV boundaries), the definition of the 802 band around the DV boundary ( Fig S2) and the region of the blade that excludes this region 803 ( Fig 1D) were not defined manually using the FIJI scripts of TissueMiner. Rather, we defined 804 them after the TissueMiner databases were generated, using the position of cells relative to 805 the DV boundary in the last frame. Using the plots of the radial elongation on each disc, we 806 estimated the width of the stripe region, and then filtered for cells included in and excluded 807 from this region. Once these cells in the last frame were identified, we used the 808 UserRoiTracking.R of TissueMiner to backtrack these two regions, producing a list of all cells 809 belonging to this lineage that are traceable forward and backward in the movie. 810 For the static RNAi data, we used the average width of the band around the DV 811 boundary from timelapse data as an estimate and excluded this region from the static images 812 to analyze the radial gradients in cell area and elongation (Fig 5 and 7). 813 Plotting area and elongation versus distance to the center of elongation:

814
For the timelapse data, we defined the radial elongation center (described above) for 815 each movie and then calculated the distance away from this center for each cell. After 816 excluding the band of cells around the DV boundary, we binned cells by radius by rounding 817 to the nearest 10 . We also binned the movie across five equal time windows (~2 ℎ ), 818 excluding the adaption phase. Average cell area and area-weighted average of cell 819 elongation were calculated for each of these bins in each time group for each movie. 820 Including data from all 5 movies, we report the global average and its standard deviation ( Fig  821  1E-F). For Fig 3, we averaged over the last ~5 ℎ . For the band of cells around the DV 822 boundary (Fig S2), we binned cells in , rather than in radius, and report the gradient along 823 the axis, defined as the DV boundary. 824 The static RNAi data were analyzed similarly. We report the global average of all discs 825 in the genotype and the standard deviation (Figs 5 and 7). The number of discs analyzed in 826 each genotype is listed in the figure legend. 827 Regional analysis of isotropic tissue deformation 828 Spatial maps of cellular contributions to isotropic tissue deformation:

829
To analyze the pattern of isotropic deformation, we locally averaged cell behavior by 830 dividing the tissue into a grid centered upon the crossing of the AP and DV boundaries. First, 831 we determined the average position of cells in the AP and DV boundaries to generate a 832 common frame of reference across all movies. Second, we divided the tissue visible in the 833 last frame into a grid centered on these compartment boundary positions, with grid size = 834 8 . Grid boxes that were incompletely filled (less than 33% of the area of the box filled by 835 cells) were discarded to eliminate noise along the tissue border. Third, each grid box was 836 considered an "ROI" and then tracked through the entire movie using TissueMiner's ROI 837 tracking code. Fourth, the rate of deformation by each type of cellular contribution was 838 calculated as a moving average (kernal = 11) for each grid position for each timestep. Last, 839 we averaged these rates over all time points post-adaption period and plotted in space (Fig  840  1H,I,K,L). 841 Plotting the radial profile:

842
To show tissue deformation as a function of distance to the center, we first calculated 843 the distance to the center of symmetry for each cell. Second, we divided the tissue visible in 844 the last frame into radial bins, rather than a grid, by rounding the radial position of each cell 845 to its nearest 10 . Third, as we did for the grid, we defined each radial bin to be an "ROI" 846 and then tracked the region forward and backward using TissueMiner's ROI tracking code. 847 Fourth, the rate of deformation by each type of cellular contribution was calculated as a 848 moving average (kernal = 11) for each radial bin ROI at each timestep. Last, we averaged 849 these rates over all time points post-adaption period and plotted this rate as a function of 850 distance to the center. Note that we also show the spatial profiles of tissue growth in the 851 band around the DV boundary in Fig S2, but here we binned along , not . 852 Regional analysis of anisotropic tissue deformation 853 Spatial maps of cellular contributions to anisotropic tissue deformation:

854
We previously published patterns of cellular contributions to anisotropic tissue shear 855 from movies 1-3 (Dye et al., 2017); however here, we calculated these patterns in a more 856 accurate way and report averages over multiple movies (Fig 2A). Previously, we assigned a 857 grid at each timepoint as in (Etournay et al., 2015. While this method provides a 858 simple first approximation of the pattern, it is not completely accurate because we are not 859 tracking the box in time but reassigning it at each timepoint; thus, cell movement in/out of the 860 box is not counted. Here, we perform grid box tracking, as described above for the 861 calculation of cellular contributions to isotropic tissue deformation but with grid box size = 862 15 . Further, we present a global average of not just movies 1-3, but also the two new 863 movies. We calculated for each movie the rate of tissue deformation by each cellular 864 contribution as a moving average (kernal = 11) in time. Then, we averaged across all five 865 movies at each timepoint and presented the pattern as a cumulative sum, starting from the 866 end of the adaption phase (first 2 ℎ ) and continuing through the end of the timelapse. All 867 correlation terms were added together (see (Merkel et al., 2017)). The contribution to tissue 868 shear from cell extrusion is very small and thus not shown. 869 Accumulating cellular contributions to radial tissue shear: 870 We calculated a moving average (kernal size = 11) total value for each cellular 871 contribution to radial tissue shear, averaged across the blade (excluding the band around the 872 DV boundary). Then, we accumulated the contributions after the adaption phase (first 2 ℎ ) 873 and plotted over time (Fig 2B). In addition to the previously described types of 874 correlations (Merkel et al. 2017), the radial decomposition also involves a correlation 875 between cell area and shear (see Theory Supplement part III). We added all correlation 876 terms together in Fig 2B. 877 Quantification of Circular Laser Ablation 878 Measurement of Stress, Cell Elongation, and Cell Area:

879
For the ablation data, we first projected the Z-stack images of the wing discs before and 880 after ablation using our custom surface projection (Dye et al., 2017). 881 To quantify cell elongation and area, we used the images taken before the cut. Cells 882 were segmented using the FIJI plugin TissueAnalyzer and then processed with TissueMiner 883 to rotate the disc to a common orientation (anterior to the left; dorsal up) and to create a 884 triangle network and a database structure. The area-weighted average of triangle elongation 885 and the average cell area was calculated for those cells included in the center of a circle of a 886 radius = 1.3*cut radius (4.55 ) centered at the center of the cut. Varying the size of this 887 region only slightly affects the result: <1.0*cut radius, the distribution becomes more noisy 888 because we are averaging a smaller region with less cells, but regions that are too big 889 (~1.8*cut radius) may begin to include cells that span different regions of the wing (i.e., the 890 band around the DV boundary). In Fig 4C and Fig S4B,E, we plot cell elongation projected 891 onto the stress axis. 892 To measure the shape of the tissue after ablation, we fit ellipses to the inner and outer 893 piece left by the cut using the projected after-cut images. These images were first rotated 894 using the transformation performed on the before-cut images to orient the tissue. We used 895 Ilastik (v. 1.2.2) (Berg et al., 2019) to segment the cut region: we trained the classifier using 896 all the data from a single day's acquisition, delineating three regions: cells, membranes, and 897 dark regions. Using the trained classifier, we then generated a thresholded binary image of 898 the cut tissue's shape and cropped the image around the cut. We then fit inner and outer 899 ellipses to these images. 900 Grouping ablations into regions of the wing pouch:

901
To quantify the position of the cut in the wing pouch, the immunofluorescence post-cut 902 images were projected using maximum intensity. In FIJI, we manually measured (in ) the 903 distance from the center of the cut to the DV boundary and AP boundary. We noted that 904 fixation causes the tissue to shrink. We estimated the extent of shrinkage using a subset of 905 the samples in which the compartment boundaries were visible in the Z-stack taken of the 906 live sample before the cut. We used FIJI to measure distances from the center of the cut to 907 the compartment boundaries in both the pre-fix live images and in the post-fix 908 immunofluorescence images for this subset. We measure a discrepancy indicating that the 909 tissue shrinks by ~15%. Thus, we compensate for this shrinkage when measuring distances 910 to the compartment boundaries in our analysis. We also noted the compartment 911 (dorsal/ventral/anterior/posterior) in post-fix images and added a sign to the distances to the 912 boundaries (ventral and posterior getting negative "distances") to create the spatial map 913 shown in Fig 4B.  914 To classify the cuts by position, we first estimated the width of the band around the DV 915 boundary from the timelapse data to be ~22 (centered on the DV boundary). We then 916 classified a cut as in the DV boundary if its cut boundary extended no more than 1.4 (20% 917 of the 7 radius) outside this 22 horizontal stripe region. Likewise, cuts were classified 918 as outside this region if it did not extend more than 1.4 into the 22 horizontal stripe.

919
We excluded all cuts that appear to straddle the border between these regions (grey in Fig  920  4B). We also excluded those cuts with centers lying within 7 of the AP boundary, in case 921 material properties at this boundary region are also different. 922 ESCA: Determination of stress from response to circular ablation:

923
To obtain the normalized shear stress /(2 ), normalized isotropic stress / and 924 the ratio of elastic constant 2 / , we fit ellipses to the inner and outer tissue outlines 925 simultaneously. These three parameters determine the large and small semi-axes of the two 926 ellipses. Other fitting parameters are the center positions of the two ellipses and the angles 927 of major axes. 928 Fits were performed on all of the laser ablations in the same region (either within or 929 outside of the band around the DV boundary) for a range of fixed values of the ratio 2 / . 930 The optimal value was considered to be the one that minimizes the sum of fit residuals (Fig  931  S4C,F). We defined a threshold of fit residual, as shown in Fig S4A,    Position relative to AP (μm) Position relative to DV (μm) Anisotropic tissue shear and its decomposition into cellular contributions during mid-third instar ( Position relative to AP Boundary (μm) Position relative to DV Boundary (μm)

Before Ablation
After Ablation Before Ablation After Ablation Before Ablation After Ablation Before Ablation After Ablation 1. Select region 2. Measure cell elongation 3. Image 2' after ablation 4. Fit ellipses 5. Infer mechanical properties