Self-organized patterning of cell morphology via mechanosensitive feedback
Abstract
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.
eLife digest
During development, carefully choreographed cell movements ensure the creation of a healthy organism. To determine their identity and place across a tissue, cells can read gradients of far-reaching signaling molecules called morphogens; in addition, physical forces can play a part in helping cells acquire the right size and shape. Indeed, cells are tightly attached to their neighbors through connections linked to internal components. Structures or proteins inside the cells can pull on these junctions to generate forces that change the physical features of a cell. However, it is poorly understood how these forces create patterns of cell size and shape across a tissue.
Here, Dye, Popovic et al. combined experiments with physical models to examine how cells acquire these physical characteristics across the developing wing of fruit fly larvae. This revealed that cells pushing and pulling on one another create forces that trigger internal biochemical reorganization – for instance, force-generating structures become asymmetrical. In turn, the cells exert additional forces on their neighbors, setting up a positive feedback loop which results in cells adopting the right size and shape across the organ. As such, cells in the fly wing can spontaneously self-organize through the interplay of mechanical and biochemical signals, without the need for pre-existing morphogen gradients.
A refined understanding of how physical forces shape cells and organs would help to grasp how defects can emerge during development. This knowledge would also allow scientists to better grow tissues and organs in the laboratory, both for theoretical research and regenerative medicine.
Introduction
During morphogenesis, tissues with complex morphologies are formed from the collective interplay of many cells. This process involves spatial patterns of signaling activity, and previous work has discovered mechanisms for generating tissue-scale patterns of activity in signaling pathways such as Hedgehog, TGFβ, and Wnt (Green and Sharpe, 2015). In addition, patterns of cellular morphology arise during morphogenesis. Such patterns can be important for ensuring the function of the resulting tissue. For example, the compound eye of Drosophila consists of hundreds of ommatidia organized in a precise hexagonal array that is required to fully sample the visual field (Kumar, 2012). Patterns of cellular morphology that arise during morphogenesis can also guide the morphogenetic processes itself. For example, spatial patterns of cell morphology emerge during growth of the Drosophila larval imaginal discs, which are precursors of adult tissues (Aegerter-Wilmsen et al., 2010; Condic et al., 1991; Legoff et al., 2013; Mao et al., 2013). These patterns have been proposed to be involved in the eversion process, during which these flattened epithelial sacs turn themselves inside out when the animal transitions from larva to pupa (Condic et al., 1991). While extensive work has studied the emergence of biochemical signaling patterns, how patterns of cellular morphology arise during tissue development is poorly understood.
Here, we investigate tissue-scale patterning of cell morphology in the Drosophila larval wing imaginal disc, which has a geometry that is ideal for studying spatial patterns of epithelial cell morphology. We focus specifically on the cell shape patterns in the central ‘pouch’ region, which is the precursor of the adult wing blade. To a good approximation, this region is planar and ellipsoidal. Cells near the center have smaller cell areas and are more isotropic in shape, whereas cells near the periphery have larger cell areas and are elongated tangentially (Aegerter-Wilmsen et al., 2010; Legoff et al., 2013; Mao et al., 2013). Cell shape has been correlated with mechanical stress: tangentially oriented bonds of elongated cells in the periphery are under higher tension than radially oriented bonds (Legoff et al., 2013).
It has been previously proposed that this pattern of cell morphology in the wing pouch could stem from differential proliferation: if the center grows faster than the rest, the resulting area pressure could stretch peripheral cells tangentially (Mao et al., 2013). Indeed, there is evidence to suggest that cells divide slightly faster closer to the center during very early stages (before after egg laying, AEL). It was suggested that this early growth differential is sufficient to account for the persistence of the cell morphology pattern through the remaining of development. However, it has since been shown that cell rearrangements occur (Dye et al., 2017; Heller et al., 2016), which could relax stress patterns once growth has become uniform. Furthermore, stress patterns may even relax during homogeneous growth in the absence of cell rearrangements (Ranft et al., 2010). Thus, it remains unclear how cell morphology patterns generated early by differential growth could be maintained through later stages, and alternative mechanisms for the establishment of these patterns must be considered.
Here, we measure the spatial patterns of cell morphology, cell divisions, and cell rearrangements during the middle of the third larval instar (starting at AEL). We quantify the pattern of tangential cell elongation and show that it becomes stronger over time, even though growth is spatially uniform and cell rearrangements are frequent. Strikingly, this change in tangential cell elongation is coupled to a radially biased pattern of cell neighbor exchanges. Using a physical model of tissue dynamics, we show that active patterning of radial cell neighbor exchanges can account for the observed morphology patterns in the absence of differential growth. Lastly, using a combination of experiment and theory, we provide evidence that this active patterning could be self-organized by mechanosensitive feedback.
Results
Cell morphology patterns can persist and strengthen in the absence of differential growth
Cell morphology patterns in the wing disc have been previously analyzed using static images (Aegerter-Wilmsen et al., 2010; Legoff et al., 2013; Mao et al., 2013). However, relating cell morphology patterns to patterns of growth, cell divisions, and cell rearrangements requires dynamic data. We therefore performed long-term timelapse imaging of growing explanted wing discs using our previously described methods (Dye et al., 2017), starting at AEL and continuing for of imaging. We used Ecadherin-GFP as an apical junction marker (Figure 1A).
To quantify cell morphology, we averaged apical cell area and cell elongation locally in space, using data from five wing discs, and in time using intervals. Over the course of the time-lapse, the qualitative features of the morphological patterns do not change. Therefore, in Figure 1B–C, we present the spatial patterns calculated for the middle timepoint. Cell area is represented as a color code. Cell elongation is characterized by a tensor , which defines an axis and a strength of elongation and is represented by bars in Figure 1C (see Materials and methods, Figure 1—figure supplement 1B; Etournay et al., 2015; Merkel et al., 2017). To quantify the radial symmetry of this pattern, we first determined the center of symmetry (Figure 1—figure supplement 1 and Materials and methods). We then introduce a polar coordinate system at the center with the radial coordinate and present the radial projection of the cell elongation tensor as a color code in Figure 1C (see also Figure 1—figure supplement 1). This figure highlights the pattern of tangential cell elongation, with cells elongating on average perpendicular to the radial axis (blue in Figure 1C). It also reveals that this pattern is interrupted around the Dorsal-Ventral (DV) boundary, where cells are elongated parallel to this boundary (red in Figure 1C). We quantify this region separately (Figure 1—figure supplement 2) and exclude it from our analysis of the circular patterns (Figure 1D–F). The Anterior-Posterior (AP) boundary also affects cell morphology (Landsberg et al., 2009), but the effect is weaker and more variable at this stage; thus, we do not quantify it separately.
The spatial maps of cell morphology reveal that both cell area and cell elongation magnitude are largest at the periphery and decrease toward the center (Figure 1B–C). We quantified this radial gradient in cell area and observe that cell area ranges from when moving toward the periphery (Figure 1E). We also observe a radial gradient in cell elongation starting from at and extending to about at (Figure 1F). The negative value corresponds to tangential elongation (Figure 1C). When evaluated over the timelapse, we find that these radial gradients grow slightly more pronounced over time (Figure 1—figure supplement 3C–D).
As previously proposed, differential growth can generate such patterns of cell elongation (Mao et al., 2013). However, indirect metrics of tissue growth in vivo do not indicate that differential growth still occurs at this later stage (Mao et al., 2013). We directly measured the spatial pattern of growth during timelapse. Cell division rate has been previously used as an indicator of growth; however, tissue growth actually results from a combination of cell division, cell area changes, and cell extrusions (Figure 1G). Thus, we evaluated the spatial pattern (Figure 1H,I,K,L) and radial profiles (Figure 1J,M) of total tissue growth and its cellular contributions. These data show that tissue area growth, as well as cell division rate, are to a good approximation independent of the distance to the center.
In summary, we quantified cell morphology patterns in mid-third instar wing explants during live imaging. We confirm previous static observations of the pattern and further identify a region around the DV boundary with a morphological pattern that differs from the rest of the wing pouch. We quantify the radial gradients in cell area and cell elongation existing outside of the DV boundary region and show that they strengthen in time in the absence of differential growth, raising the question of what mechanism underlies the persistence of these cell morphology patterns during mid to late stages of wing growth.
Radially oriented cell rearrangements balance tangential cell elongation
To directly relate the observed cell morphology patterns with cell rearrangements, we next analyzed the spatial patterns in cellular contributions to tissue shear (Figure 2A). Tissue shear can be decomposed into contributions from cell divisions, cell elongation changes, T1 transitions, and so-called correlation effects (Etournay et al., 2015; Merkel et al., 2017). Here, correlation effects result mainly from correlated fluctuations in cell elongation and cell rotation (see Appendix 3 and Merkel et al., 2017).
We find that the spatial patterns of tissue growth and its cellular contributions exhibit overall anisotropies perpendicular to the DV boundary, as reported previously (Dye et al., 2017). In addition, the patterns of cell elongation changes and T1 transitions can be described as a superposition of a uniformly oriented pattern and an approximately radial or tangential pattern (Figure 2C–D). To determine the magnitude of the radial or tangential patterns in all quantities, we quantified their average radial projections as cumulative plots over time (Figure 2B). The radial component of tissue shear is small, and cell divisions do not contribute to radial tissue shear. In contrast, we observe a pronounced buildup of a tangential pattern of cell elongation accompanied by a radial pattern of T1 transitions and of correlation effects.
As shown above (Figure 1J), tissue area growth does not have a radial gradient and thus does not contribute to the increase in tangential cell elongation that we observe at this time (Figure 2A–B, Figure 1—figure supplement 3C–D). Furthermore, we observe numerous T1 transitions (on average ), and their radially biased orientation increases rather than relaxes tangential cell elongation (Figure 2E–F). Thus, we are not observing the relaxation of a pattern of cell elongation caused by early differential growth. Rather, our data support a model whereby a radially patterned morphogenetic cue actively biases the direction of T1 transitions and consequently the complementary pattern of cell shape changes.
Polarity-driven cell rearrangements can create the observed cell morphology patterns in the wing disc
We next apply a biophysical model to determine whether radially patterned T1 transitions could account for the observed cell morphology patterns in the wing disc. This model takes into account the interplay of T1 transitions, cell shape changes, and tissue shear in a continuum description (Etournay et al., 2015; Popović et al., 2017). Active anisotropic force-generating processes that bias cell rearrangements in the tissue, such as polarization of the actomyosin cytoskeleton, are captured by a nematic cell polarity , defined by a magnitude and an orientation axis. We propose that such a patterning cue leads to the radially oriented pattern of T1 transitions we observe in the wing disc.
In our model, we consider the spatial patterns of tissue shear rate , the patterns of cell shape , and the patterns of cell rearrangements , which have nematic symmetry and are represented by traceless symmetric tensors that describe the magnitude and orientation. Tissue shear is defined as a velocity gradient tensor that results from a combination of cell shape changes and cell rearrangements (Etournay et al., 2015; Merkel et al., 2017):
Here, is a co-rotational time derivative of , and the shear due to cell rearrangements includes contributions from T1 transitions, cell divisions and extrusions, and also correlation effects.
Tissue material properties are described by constitutive equations for the tissue shear stress and the shear due to cell rearrangements (Figure 3A–B):
Here, the shear stress tensor is the sum of the elastic and active stress. The elastic stress is associated with cell elongation and characterized by the shear elastic modulus . The active stress associated with the cell polarity cue (Figure 3A) is described by the coefficient . The shear rate due to cell rearrangements given by Equation 3 is driven in part by shear stress, and therefore depends on the cell elongation , and in part by active processes that are oriented by the nematic cell polarity (Figure 3B). is a relaxation time for cell rearrangements over which elastic stresses are relaxed, and is the rate of cell rearrangements driven by cell polarity. In this coarse-grained picture, subcellular processes are captured by effective coefficients. For simplicity, we do not include dissipative processes on cellular scales, such as cytoskeletal viscosity or cell-cell friction. Such dissipative processes are relevant on short timescales and are small in comparison with elastic stresses on tissue relevant timescales.
To discuss the wing disc, we consider a radially symmetric geometry and average the oriented quantities after projection onto the radial axis. Radial tissue shear is small compared to that associated with cell shape changes and T1 transitions during our observed time window (Figure 2B). We therefore consider a steady state with , , and . In this case, cell elongation becomes:
Thus, we find that the steady state cell elongation pattern is a result of cell rearrangements that are oriented by the cell polarity cue (Figure 3C–D). Note that our data show that the wing disc is not exactly at steady state: cells slowly change their shape and rearrange radially (Figure 2A–B). However, as we show in Appendix 1 part 1-2, Equation 4 holds to a good approximation.
Can the radial pattern of T1 transitions defined by also explain the observed radial profile of cell area (Figure 1B,E)? To answer this question, we then considered force balances in the tissue. We consider tissue area pressure:
where is the difference in pressure from a reference value, is tissue area compressibility, is the average cell area, and is a reference cell area. As pressure increases, cell area decreases. To calculate the cell area profile, we again approximate the wing pouch as a radially symmetric disc. In the radially symmetric geometry, force balance can be expressed as:
A radial profile of pressure determined from this equation implies a radial pattern of cell area via Equation 5 (see Figure 3C–E and Appendix 1 part 1-2). To test this implication, we first quantify the radial profile of cell elongation to estimate the profile of using Equation 4 (Figure 3F). We represent the cell elongation data in a functional form using an empirical power law that is fit to the data. Then, using this functional form in Equations 2, 5, and 6, we solve for the cell area profile (Figure 3G and Appendix 1 part 1-2). Finally, we show that this function can account for the observed pattern of cell area (see Figure 3G).
From this analysis, we conclude that the cell morphology patterns observed in the wing disc could be generated by radially biased cell rearrangements. Next, we test whether the stress profile predicted by the model (Equation 2) exists in the tissue, and we measure key mechanical parameters of the model. Later, we address the potential molecular origin of the cell polarity cue orienting the cell rearrangements.
Circular laser ablation reveals patterns of tissue stress
Our model predicts a stress pattern in the wing disc that results from active processes that are radially oriented by a cell polarity cue. To compare this prediction to experiment, we infer tissue stress using laser ablation. Tissue stress has been estimated previously by laser ablation techniques that are based on determining the initial retraction velocity (Bonnet et al., 2012; Etournay et al., 2015; Farhadifar et al., 2007; Legoff et al., 2013; Mao et al., 2013; Shivakumar and Lenne, 2016). However, to compare theory and experiment, ideally one should measure quantities that are well-captured by the model. Therefore, instead of using initial retraction velocity, we perform circular cuts and analyze the final, relaxed position of the inner and outer elliptical contours of tissue formed by the cut (Figure 4A, Appendix 2). From the size and the anisotropy of the cut, we can infer anisotropic and isotropic tissues stress, normalized by the respective elastic constants, as well as the ratio of elastic constants (see Appendix 2). Furthermore, we can also infer the existence of polarity-driven stress. We name this method ESCA (Elliptical Shape after Circular Ablation).
We perform local measurements of tissue mechanics by cutting the tissue in the smallest possible circle that would still allow us to measure the shape of the inner piece left by the cut (radius = , encircling cells, Figure 4A). To relate measurements of tissue stress to cell elongation, we calculate the average cell elongation in the ablated region before it is cut (Figure 4A). In Figure 4B, we present the measured cell elongation and shear stress tensors at the position of ablation. We find that local average cell elongation correlates well with the principal direction of shear stress. Also, we observe that the cells in the band around the DV boundary, which are exposed to high Wg and Notch signaling, have different mechanical properties than elsewhere in the tissue. Near the DV boundary, cells elongate less than outside this region for comparable amounts of stress (Figure 4B and Figure 4—figure supplement 1B,E). The ratio of elastic constants in this region is also smaller: near the DV boundary , whereas outside this region, (see also Figure 4—figure supplement 1C,F). We focus hereafter on the radial patterns of elongation and stress outside of the DV boundary region.
The relationship between cell elongation and stress normalized by the elastic modulus has a slope 1 in the absence of polarity-driven stress (see Equation 2). We observe a much smaller slope for this relationship in our data (Figure 4C), indicating that polarity-driven stress is significant. We now use these data to estimate the parameters of our mechanical model. We write the shear stress defined in Equation 2 in terms of cell elongation and cell rearrangements, eliminating the orientational cue using Equation 3. For the radial components, we have:
where * = is an effective shear elastic coefficient. The difference between and depends on the parameters and associated with the nematic cell polarity. We fit Equation 7 to the data and find and (Figure 4C, see Appendix 1 part 3). Combined with data from Figure 1, we find an estimate for the tissue relaxation time , which is roughly consistent with that found during pupal morphogenesis (Etournay et al., 2015). From our data, we can also infer the radial profile of tissue area pressure, revealing that pressure increases toward the center (Figure 4D and Appendix 1 part 3). This finding is consistent with the observed cell area profile, with smaller cell areas toward the center (Figure 1B,E).
In sum, we find a stress profile in the wing disc that is consistent with the observed measurements of both cell elongation and area. Further, we use these data to measure certain parameters of our biophysical model, including the tissue relaxation timescale and the effective shear elastic coefficient.
Reduction of planar cell polarity pathways does not reduce tangential cell elongation
In our model, the radial orientation cue is required to generate the observed patterns of cell morphology, cell rearrangement, and tissue stress. Candidates for such an orientational cue are the planar cell polarity pathways (PCP), which are groups of interacting proteins that polarize within the plane of the epithelium. There are two well-characterized PCP pathways: Fat and Core (Butler and Wallingford, 2017; Eaton, 2003). In the wing, these systems form tissue-scale polarity patterns during growth (Brittle et al., 2012; Merkel et al., 2014; Sagner et al., 2012) and are required to position the hairs and cuticle ridges on the adult wing (Adler et al., 1998; Doyle et al., 2008; Eaton, 2003; Gubb and García-Bellido, 1982; Hogan et al., 2011). To determine whether either of these pathways could function as the orientational cue described in our model, we analyzed cell elongation patterns after their removal.
We perturbed the Fat pathway using nub-Gal4 to drive the expression of RNAi constructs targeting both fat (ft) and dachs (d) in the pouch region throughout the third larval instar. This perturbation results in almost complete loss of Dachsous from the apical membrane, and any residual signal is no longer polarized, confirming the loss of PCP (Figure 5—figure supplement 1A–B). Furthermore, we observe a suppression of tissue growth upon ft+d double RNAi knockdown (visible at the end of larval development in Figure 5A and in the resulting adult wings in Figure 5—figure supplement 1C), consistent with previous work on the loss of both Dachs and Fat (Cho and Irvine, 2004). Nonetheless, the pattern of tangential cell elongation persists to the end of larval development (Figure 5A–C). Using scaled coordinates, we find that the radial profiles of cell elongation in ft+d RNAi and control wings are similar (Figure 5D).
-
Figure 5—source data 1
- https://cdn.elifesciences.org/articles/57964/elife-57964-fig5-data1-v3.zip
-
Figure 5—source data 2
- https://cdn.elifesciences.org/articles/57964/elife-57964-fig5-data2-v3.zip
-
Figure 5—source data 3
- https://cdn.elifesciences.org/articles/57964/elife-57964-fig5-data3-v3.zip
We perturbed the Core PCP pathway using a previously characterized null mutation in prickle (pk30), which causes defects in adult wing hair orientation (Gubb et al., 1999). We found that the cell elongation pattern in the pk mutant is similar to the wild-type control (Figure 5E–I). In the pk mutant, the region of tangential cell elongation extends even further into the center than in control wings.
We conclude that the tangential cell elongation pattern persists in the absence of either PCP pathway. This result excludes these pathways as orienting cues for the cell elongation patterns.
Mechanosensitive feedback generates self-organized patterns of cell morphology
We have shown that perturbing PCP pathways does not affect the radial patterns of morphology, raising the question of how orientational cues might arise. In previous sections, we have considered the orientational cue to be provided by a cell polarity system that is independently patterned. However, cell polarity in general would be affected by stresses in the tissue. Indeed, there are many examples of cells polarizing in response to mechanical stress (Duda et al., 2019; Hirashima and Adachi, 2019; Ladoux et al., 2016; Ohashi et al., 2017). Here, we show that introducing mechanosensitive feedback to the model of tissue mechanics can give rise to spontaneous emergence of the cell polarity cue (Figure 6).
Mechanosensitivity is incorporated into our model of tissue mechanics through a dynamic equation for the orientational cue that becomes stress-dependent:
Here, is a relaxation timescale for , is a mechanosensitive feedback strength, the coefficient ensures stability, and is a coupling strength locally aligning orientational cues.
Now, Equation 8 provides a mechanosensitive feedback to Equations 1, 2 and 3. These combined equations show a novel behavior. Specifically, the orientational cue can emerge spontaneously by self-organization (Figure 6A–B). Beyond a critical value of the mechanosensitive feedback strength , an isotropic tissue with is no longer stable, and a state with an orientational cue emerges instead (Figure 6B and Appendix 1 part 4). The magnitude of this spontaneous polarization is , where , where a positive coefficient is needed to stabilize the polarized state. By this mechanism, the anisotropic cue introduced earlier in our model can be locally generated by mechanosensitive self-organization and does not require the existence of pre-patterned polarity cues. To generate a large-scale pattern from locally generated anisotropic cues, they need to be aligned in neighboring regions. This local alignment is captured in Equation 8 by the orientation coupling term with strength , which is similar to alignment terms found in anisotropic physical systems, such as liquid crystals (Gennes and Prost, 1993; Jülicher et al., 2018; Marchetti et al., 2013).
To discuss cell morphology profiles in the wing disc, we consider a simplified tissue model with radial symmetry, where the rate of radial cell rearrangement is given (as estimated in Appendix 1 part 2) and the cell shape pattern and tissue stress pattern are calculated. Using a fit of cell elongation to the experimental data, we find a set of parameter values that accounts for the observed cell elongation patterns in the wing disc (Table 1, Appendix 1, Figure 6C). From this cell elongation pattern also follows the cell area pattern (as described above, Figure 3).
We conclude that our mechanosensitive model can account for the radial pattern of cell morphology in the wing disc. Due to the relatively large number of parameters used to fit a single experimental curve, there are large uncertainties when estimating parameter values. Nonetheless, these uncertainties do not affect the qualitative prediction that a reduction in mechanosensitivity would lead to less polarization and thereby reduced cell elongation (see Appendix 1 part 4). We next test this prediction of our model experimentally.
Suppression of mechanosensitivity weakens the gradients in cell elongation and cell size
In order to test our prediction that the reduction of mechanosensitivity will reduce the magnitude of cell elongation, we used RNAi to reduce the levels of Myosin VI (MyoVI), a molecular motor implicated in mechanosignaling. MyoVI, encoded by jaguar (jar) in Drosophila, is an upstream component of a Rho-dependent signaling pathway that reorganizes the actin-myosin cytoskeleton in response to mechanical stress (Acharya et al., 2018). Experiments in wing discs also indicate that mechanosensation involves Rho polarization and signaling (Duda et al., 2019). We performed RNAi targeting MyoVI in the wing pouch using nub-Gal4 and evaluated cell morphology at the end of larval development ( AEL). We observe a clear reduction in the magnitude of tangential cell elongation as compared to wild type at this stage (Figure 7A–D, Figure 7—figure supplement 1). In addition, our model predicts that such reduction of cell elongation would result in an increase of cell area in the wing (Equations 1-6). The observed pattern of increased cell area upon reducing MyoVI levels with RNAi is consistent with this prediction (Figure 7E–F). Therefore, the qualitative predictions of our model upon reducing the mechanosensitive feedback strength are confirmed by the experimental downregulation of the mechanosensitive motor MyoVI.
Discussion
Here, we have shown that patterns of cell shape and stress in the mid-third instar Drosophila wing disc do not rely on PCP pathways or differential growth. Instead, radially-oriented T1 transitions and tangential cell elongation emerge via mechanosensitive feedback in a self-organized process. We have presented a continuum model of tissue dynamics for this self-organization based on a mechanosensitive nematic cell polarity that accounts for the observed patterns of cell area, T1 transitions, and cell shape. Our work highlights a mechanism for the self-organized emergence of cellular patterns in morphogenesis, expanding our understanding of pattern formation emerging from mechanical feedbacks in active systems (Bois et al., 2011; Howard et al., 2011; Recho et al., 2019).
A pattern of T1 transitions is critical for cell morphology patterning in the Drosophila wing
Our work shows that the spatial pattern of T1 transitions is an integral part of the emergence of tissue organization during wing development. In contrast to situations such as germband extension, where T1 transitions exhibit clearly discernible patterns, the patterns of T1 transitions in the wing disc have been elusive. Many T1 transitions occur in the tissue in seemingly random orientations. However, on average, they exhibit a spatial pattern. We revealed these patterns by quantifying the nematics of T1 transitions and cell shape changes using the previously-described triangle method (Merkel et al., 2017) and then quantified them with radial averaging (Figure 2A–B). In this way, we revealed that a radial pattern of T1 transitions is linked to a tangential pattern of cell elongation.
Given this radial pattern of T1 cell rearrangements, the observed cell morphology pattern follows from a continuum tissue model based on a radially oriented nematic cell polarity field (Figure 3). The polarity-oriented radial T1s create a cell shape pattern with corresponding patterns of tissue stress and tissue area pressure. The 2D area pressure is higher in the center and is lower toward the periphery. Note that this pressure profile does not rely on differential proliferation, as was previously proposed (Mao et al., 2013) but instead relies on a radial pattern of T1 transitions. We test our model using a novel circular laser ablation method. This method allows us to determine specific combinations of tissue parameters. In particular, we estimate the ratio of elastic constants and , as well as the cell shape relaxation timescale .
This analysis raised the question of which nematic cell polarity cues guide the cell rearrangement and cell elongation patterns. PCP pathways are required for the proper orientation of T1 transitions in other contexts (Bosveld et al., 2012). However, we found that neither of the two known PCP pathways in the wing are required for the observed tangential cell elongation patterns (Figure 5). We instead show that an orientation cue can arise through self-organization via mechanosensitivity and identify MyoVI as a key molecular player.
Mechanosensitive feedback can create self-organized patterns of cell morphology
Cell polarity cues can emerge via mechanosensitive feedback by transforming mechanical cues into chemical anisotropies. Nematic cell polarity can then orient active stresses and thereby amplify the mechanical stimulus (Figure 6). We introduce this mechanosensitive feedback in our continuum theory, which quantitatively describes the emergence of patterns of cell shape and cell rearrangements. The strength of this mechanosensitive feedback is described by a parameter . If exceeds a critical value , an orientation cue and elongated cell shapes spontaneously emerge (Figure 6B). This model can account for the observed patterns of cell area and cell elongation in the wing disc and predicts that the reduction of mechanosensitivity will result in reduced cell elongation. To test this prediction, we perturbed a RhoA-dependent mechanotransduction pathway by lowering levels of an upstream component, MyoVI, using RNAi (Figure 7). We find a clear phenotype of reduced cell elongation and increased cell areas in the center region, as predicted by our model. In the spatially resolved model, which includes the term coupling polarity orientation of neighboring cells (Figure 6C), a polarity pattern can be induced even below by imposing polarity at the tissue boundaries. Thus, the residual pattern of cell elongation that we observe after removal of MyoVI in the wing pouch could be due to polarity existing outside of this perturbed region. In addition, the residual pattern could also indicate an incomplete knockdown from RNAi or the presence of other mechanosensitive elements in the tissue. Nevertheless, given the clear phenotype that is fully consistent with our model, we propose that the mechanosensitive feedback mechanism is a significant determinant of the cell shape patterns in the wing pouch.
Our data, together with the fact that MyoVI is involved in Rho-dependent activation of actin-myosin cytoskeleton (Acharya et al., 2018), suggest that MyoVI is a molecular component of the mechanosensitive feedback we describe in our self-organized model. However, the molecular nature of the cue that defines the nematic cell polarity is unknown. This cue may organize the structure or dynamics of the actin-myosin cytoskeleton or the actin-myosin system itself could define nematic cell polarity. Indeed, it has been shown that Myosin II (MyoII) localizes to long cell boundaries in the wing (Legoff et al., 2013), corresponding to a nematic polarity aligned with the nematic cell polarity . Also, wing disc stretching experiments have shown that MyoII can polarize in response to exogenous stress in a Rho-dependent manner (Duda et al., 2019). Furthermore, it has been suggested that MyoII polarity arises as a consequence of cell stretching and functions as a negative mechanical feedback (Legoff et al., 2013), consistent with the role of in our model. Precisely how the actin-myosin cytoskeleton is affected by MyoVI in this system and how these cytoskeletal elements together guide cell rearrangements in response to anisotropic tissue stresses and cell shape changes remain open questions for future research.
Lastly, our laser ablation analysis shows that the region around the DV boundary has a different ratio of elastic constants than the rest of the tissue, which could affect the self-organized pattern formation we describe. Therefore, it will be interesting to study how Wingless/Notch signaling, which defines the DV boundary, may influence the mechanical properties that lead to mechanosensitive self-organization of polarity and morphology. In addition, we observe a richer pattern emerging very late in development (see Figure 7—figure supplement 1, Figure 5), including a region anterior to the AP boundary that is radially elongated. Future research will expand upon the model presented here to explore the dynamics of these patterns.
In summary, we used the Drosophila wing disc to identify a mechanism by which tissue morphology can arise from the self-organization of a mechanical feedback coupling cell polarity to active cell rearrangements. This mechanism is general and could be employed in other tissues and organisms to generate patterns of cell shape and cell area. Thus, we hope our work inspires new avenues of research that integrate theory and experiment to understand biological self-organization.
Materials and methods
Lead contact and materials availability
Request a detailed protocolThis study did not generate any new unique reagents. All requests for further information and reagents may be directed to the lead author, natalie_anne.dye@tu-dresden.de.
Experimental model and subject details
Request a detailed protocolAll experiments were performed with Drosophila melanogaster, using lines that are publicly available and previously published. Our Drosophila lines were fed with a standard media containing cornmeal, molasses agar and yeast extract and grown under a light/dark cycle. All experiments were performed at . Both males and females were analyzed, and the sex of the animals was not recorded, as we have no reason to believe there is any sexual dimorphism in the studied phenomenon. To synchronize development, we collected eggs deposited within a defined time window on apple juice agar plates. To do so, we transferred the flies from standard food vials to cages covered by apple juice agar plates containing a dollop of yeast paste for food. After at least , the plates were replaced, and the timing of collection started. Eggs laid within a time window were collected by cutting out a piece of the agar and transferring it to a standard food vial. We limited the number of eggs per vial to to avoid crowding. The middle of the time window for egg collection was considered to be AEL. Experiments from timelapse imaging and laser ablation (Figures 1, 2, and 4) were captured after explanting at AEL, whereas those involving RNAi (Figures 5 and 7) were explanted at AEL to allow the maximal amount of time for the RNAi phenotype to emerge. The specific genotypes used for each experiment are indicated in the Key Resources Table and in the figures.
Method details
Timelapse image acquisition and processing
Request a detailed protocolSample preparation
Request a detailed protocolWing explants were grown ex vivo as described previously (Dye et al., 2017). Briefly, wing discs were dissected from larvae in growth media (Grace’s cell culture media + 5% fetal bovine serum 20-hydroxyecdysone + Penicillin-Streptomycin) at room temperature. Then, they were transferred to a Mattek #1 glass bottom petri-dish to the center of a hole cut in a double-sided tape spacer (Tesa 5338) and covered with a porous filter. The dish was then filled with fresh growth media.
Acquisition
Request a detailed protocolData from movies 1 to 3 were used previously (Dye et al., 2017). Movie 4 was acquired after publication of the first manuscript, in the same exact way as movies 1 to 3. Briefly, Ecadherin-GFP-expressing wing discs were imaged in growth media using a Zeiss spinning-disc microscope to acquire spaced Z-stacks at 5 min intervals with a Zeiss C-Apochromat 63X/1.2NA water immersion objective and tiling (10% overlap). This microscope consisted of an AxioObserver inverted stand, motorized xyz stage, stage-top incubator with temperature control set to , a Yokogawa CSU-X1 scanhead, a Zeiss AxioCam MRm Monochrome CCD camera (set with binning), and 488 laser illumination. We circulated growth media during imaging using a PHD Ultra pump (Harvard Apparatus) at a rate of . Time of the movie was considered to be the start of imaging acquisition, typically from the start of dissection (time required for sample and microscope preparation).
Movie 5 was acquired on a newer microscope with a larger field of view that eliminated the need to tile across the wing pouch. This microscope has an Andor IX 83 inverted stand, motorized xyz stage with a Prior ProScan III NanoScanZ z-focus device, a Yokogawa CSU-W1 with Borealis upgrade, and a Pecon cage incubator for temperature control at . We used an Olympus 60x/1.3NA UPlanSApo Silicone-immersion objective with 488 laser illumination and an Andor iXon Ultra 888 Monochrome EMCCD camera. We acquired spaced Z-stacks of a single tile at intervals, as for the other movies. Movie 5 was also acquired without the constant flow of new media.
For all movies, care was taken to limit light exposure, using laser power values of and exposure times less than per image.
Processing
Request a detailed protocolRaw Z-stacks were denoised using a frequency bandpass filter and background subtraction tools available in FIJI (Schindelin et al., 2012). Then, we used a custom algorithm as described previously (Dye et al., 2017) to make 2D projections of the apical surface, marked by Ecadherin-GFP. This algorithm also outputs a height-map image, in which the value for each pixel corresponds to the level in the Z-stack of the identified apical surface. For those movies that were tiled, we used the Grid/Collection Stitching FIJI plugin (Preibisch et al., 2009) to stitch the 2D projections, and then used that calculated transformation to stitch the height map images, so that we could correct the cell area and elongation values for local curvature (see below). We focus in this work exclusively on the disc proper layer, which is the proliferating layer of the disc that goes on to produce the adult wing. We did not consider cell shape patterns in the peripodial layer, the non-proliferating layer of the wing disc that is largely destroyed at the onset of pupariation.
Single timepoint image acquisition and processing for analysis of RNAi/mutant phenotypes
Request a detailed protocolSample preparation
Request a detailed protocolSamples for single timepoint imaging were acquired exactly as described for live imaging: dissected in growth media and imaged in Mattek dishes under a porous filter.
Acquisition
Request a detailed protocolAll imaging of single timepoint data was performed with the same microscope described above for Movie 5. We did not acquire timelapse data for these genetic perturbations, and thus we chose to image the entire disc using tiling, spaced Z-stacks. Approximately discs were imaged per dish, within approximately .
Processing
Request a detailed protocolTiled images were stitched together as Z-stacks; then we obtained the apical surface projection and its corresponding height map as described above and in Dye et al., 2017.
Circular laser ablation
Request a detailed protocolAcquisition
Request a detailed protocolWing discs from mid-third-instar larvae ( AEL) were dissected and mounted exactly as was done for the live imaging timelapses. Due to constraints on the speed of ablation, we only cut regions of the wing disc pouch that were completely flat, so that we could cut in a single plane at the apical surface (rather than having to cut in each plane of a Z stack). Where this flat region lies depends on how the disc happens to fall on the coverslip during mounting. Because the anterior compartment is larger and higher, most of our ablations are in this compartment. To access other regions, we also mounted some wing discs on an agarose shelf: stripes of 1% agarose (in water) were dried onto the surface of the coverslip and wing discs were arranged their anterior half propped up on the agarose shelf prior to adding the porous filter cover.
Ablations were performed using ultraviolet laser microdissection as described in Grill et al., 2001 using a Zeiss 63X water objective. First, we took a full Z-stack of the sample prior to the cut. Then, we selected a radius circle that would ablate in the flattest region of the tissue. No imaging is possible during ablation, but we acquired a timelapse immediately after the cut in a single Z-plane. This timelapse data was not used except to estimate whether or not the sample was fully cut. After, once the sample has finished expanding but not started to heal (), we took another full Z-stack to image the endpoint. We excluded a small number of data points if any of the following were true: (1) the inner piece remaining after the ablation was no longer visible (sometimes it floats or is destroyed); (2) the cut appeared to expand highly asymmetrically (rare); (3) the wing discs were clearly too young to be considered AEL (poor staging).
Immunofluorescence after ablation
Request a detailed protocolDue to a limited field of view on the microscope used for ablation, we performed immunofluorescence after the ablation in order to better estimate the position of the cut in the wing pouch. After all the discs in the dish were ablated ( discs/dish), the entire dish was fixed through the filter by adding 4% PFA and incubating at room temperature. After, the dish was rinsed and kept in PBST (PBS + 0.5% Triton X-100) until all discs from that image acquisition day were completed (). All samples were then blocked using 1% BSA in PBST + NaCl for , and then incubated in primary antibody overnight (diluted in BBX: 1% BSA in PBST). Initially, we labeled samples with SRF (Active Motif, 1:100 dilution), but later we switched to Patched and Wingless (DSHB, 1:100) to identify the compartment boundaries. Primary antibody toward GFP (recognizing the Ecadherin-GFP) was also included at 1:1000. After overnight incubation at , we washed with BBX, followed by BBX+ 4% normal goat serum (NGS), for at least . Secondary antibodies were added for at RT in BBX+NGS. Finally, samples were washed in PBST and imaged in this media. Imaging was performed on one of the two spinning disc microscopes described above for live imaging. We matched the stained samples with the ablation images by (1) keeping track of the position of each disc on the dish and where it was ablated during acquisition and (2) morphology of the disc before/after ablation. All antibody information are listed in the Key Resources Table.
Validation of Fat/D PCP RNAi
Request a detailed protocolSample preparation
Request a detailed protocolTo immunostain larval wing discs, late third instar larvae were partially dissected in PBS, inverting the heads to expose the wing discs, and then fixed in 4% PFA for . Samples were then washed twice for each with PBX2 ( PBS + 0.05% Triton X-100), and blocked in BBX250 (PBX2 + BSA + NaCl) for . Incubation in primary antibody diluted 1:50 in BBX (PBX2 + BSA) was performed overnight at . We used a mouse monoclonal antibody against Dachsous (Merkel et al., 2014). After, the sample was washed twice with BBX, for each wash, followed by one wash with BBX + normal goat serum (NGS). Incubation with secondary antibody was performed either for at room temperature or overnight at . We used a goat anti-mouse Alexa-Fluor 647 as the secondary antibody (ThermoFisher,Cat No. A28181) at 1:500 dilution in (BBX+4% NGS). The secondary antibody was removed by washing twice with PBX2, followed by washing twice with PBS. The tissues were stored in PBS and dissected from the body wall just before mounting. Wing discs were mounted on a slide within a thin channel created by two strips of double-sided tape (Tesa 5338). Wings were transferred to the middle of the channel with a p200 pipette, the excess PBS was removed, and the tissues were arranged with apical sides up. Next, the sample was covered with a #1 coverslip, adhering the coverslip to the double-sided tape. Vectashield mounting medium (Vector laboratories) was added to one side of the coverslip and allowed to seep in by capillary action. Excess Vectashield was removed, and the sides of the coverslip were sealed using transparent nail polish. The slides were then stored at until imaging.
To image the adult wing morphology, male flies of the desired genotype were collected and stored in isopropanol. Wings were dissected and collected in isopropanol. To mount, wings were transferred to a glass slide using a p200 pipet. In a centrifuge tube, of isopropanol was mixed with of Euperal, and about of this mix was added to the slide containing the wings. The wings were then aligned on the slide, adding more isopropanol/Euperal mix when necessary to avoid drying. Once the wings were all aligned, more isopropanol-Euperal mix was added and allowed to partially dry. Once the mix was almost dry, a clean coverslip containing about of Euperal was inverted on top of the wings. The slides was allowed to cure for before imaging.
Acquisition
Request a detailed protocolImmunostained samples were imaged on an Olympus IX81 microscope equipped with a spinning disk module (Yokogawa) and back illuminated EMCCD (Andor Technology, iXon Ultra 888) with an exposure time of and EM gain of 250. A confocal z-stack of immunostained samples was acquired using a silicon immersion objective with z-spacing.
Adult wings were imaged on an inverted wide field microscope (Zeiss) using a objective. The images were analyzed using a custom code written in MATLAB (Mathworks, USA) to measure wing area.
Quantification and statistical analysis
Cell segmentation, tracking, and alignment
Request a detailed protocolSegmentation and tracking
Request a detailed protocolUsing the 2D projections, we performed cell segmentation and tracking using the FIJI plugin, TissueAnalyzer (Aigouy et al., 2016). We manually corrected errors in the automated segmentation and tracking as much as possible and then generated a relational database using TissueMiner (Etournay et al., 2016). Images were rotated to a common orientation (Anterior left, dorsal up). We then manually identified three regions of interest at the last point in the timelapse, using the FIJI macros included with TissueMiner: the ‘blade’ was defined roughly as an elliptical region surrounded by the most distal folds; and the Anterior-Posterior and Dorsal-Ventral boundaries were estimated using the brightness of Ecadherin-GFP and apical cell size (Jaiswal et al., 2006).
For the timelapse data, we used only the cells within these manually defined regions that were trackable during the entire course of the movies. An example of this region is presented for Movie 1 in Figure 1A. Furthermore, we also excluded from all of our analysis the first of imaging, the so-called adaption phase, where cells uniformly shrink in response to culture (Dye et al., 2017).
For the RNAi data (Figures 5 and 7), we did not acquire timelapses, rather a single timepoint at late stages of development. We nevertheless chose to analyze these data using the same TissueMiner workflow for simplification. Because TissueMiner was developed for timelapse data, however, it requires at least two timepoints. Thus, we duplicated the data and labeled the two images as if they were timepoints 0 and 1 of a timelapse. We then only analyzed timepoint 0. The regions of interest in these data, therefore, are manually defined and not simply the region that is trackable (since it is static data).
Alignment on compartment boundaries
Request a detailed protocolTo average across all movies of timelapse data, or all discs of the same genotype of the static RNAi data, we generated a disc-coordinate system by normalizing the position of each cell to the AP and DV boundaries for that disc. To do so, we averaged the absolute xy positions of all the cells in the AP or DV boundaries over all time after the adaption phase. We then calculated the distance of each cell to the new X axis (average position of the AP boundary) and Y axis (average position of the DV boundary).
Analysis of cell size and shape
Request a detailed protocolDefinition of the cell elongation tensor
Request a detailed protocolThe cell elongation tensor is a traceless symmetric tensor that quantifies the anisotropy of cell shapes in a region of the tissue. We define cell elongation using a triangulation of the tissue obtained by connecting centroids of connected cells (Etournay et al., 2015; Merkel et al., 2017). The state of each triangle is described by a tensor , defined by a mapping an equilateral reference triangle to the triangles in the tissue. The state tensor contains information about the triangle elongation tensor , orientation angle , and area :
Here, is the area of the reference triangle, and is a 2-dimensional rotation matrix. Cell elongation in the tissue region is defined as an area weighted average of triangle elongation. For details about the method see Merkel et al., 2017.
Adjustment of cell area and elongation to account for tissue curvature
Request a detailed protocolThe wing disc pouch has a slightly domed shape (Figure 1—figure supplement 1A). After projection onto a 2D surface, the cell shapes and areas will be distorted. To ensure that the radial profiles we measure in a 2D projection (Figure 1B,C,E and F and Figure 3F,G) are not a result of tissue curvature, we account for the distortion caused by projection. We use the height maps generated by our projection algorithm, which identifies the apical surface of the pouch within the 3D Z-stacks. We smooth this height map with a gaussian filter of width to find the height field , and then we calculate the height gradient field. We then smooth the result again with the same gaussian filter to find . The deprojected cell or triangle area is obtained from the measured area as , where is evaluated at the cell center. To find the deprojected cell elongation of a triangle we evaluate at the triangle center. Then, we find the angle of steepest ascent in the projected plane and define the tilt transformation:
We apply this transformation to the triangle shape tensor , as defined in Merkel et al., 2017, to determine the deprojected triangle shape tensor :
The cell elongation tensor is obtained from the triangle shape tensor as the corresponding area weighted average using the deprojected cell area, as described in Merkel et al., 2017.
Spatial maps of cell size and shape
Request a detailed protocolTo generate the color-coded smoothed plots of area and elongation (Figure 1B–C), we divided the aligned wings into a grid with boxsize pixels (). For each position, we averaged cell area or performed an area-weighted averaged of triangle elongation in a neighborhood box = 20 pixels (). To create the nematic elongation pattern, we similarly averaged elongation of triangles whose centers lie within each grid box, with grid box size pixels (). Tables containing the corrected cell area and triangle elongation values for each movie at each timepoint have be uploaded to Dryad: doi:10.5061/dryad.jsxksn06b.
Calculation of radial elongation center point
Request a detailed protocolWe define the center of the wing pouch to be on the DV boundary. To determine the location of the center along the DV boundary we divide the pouch into four regions defined by the DV boundary and a line perpendicular to it located at some position , as shown in Figure 1—figure supplement 1E. Then, we define a function:
for all values of along the DV boundary. Here, are the areas of the four regions and are the area weighted averages of the cell elongation component calculated in the four regions. Finally, is the value that minimizes (Figure 1—figure supplement 1F). We find that the center point lies just anterior to the intersection with the anterior-posterior (AP) boundary (Figure 1, Figure 1—figure supplement 1, consistent with Legoff et al., 2013).
In time-lapse experiments, cell centers were determined using the last 20 timepoints. In single timepoint image experiments (Figures 5 and 7), all images of a single genotype were used together after alignment on the DV and AP boundaries.
Exclusion of the band of cells near the DV boundary
Request a detailed protocolIn contrast to the other regions (blade, AP and DV boundaries), the definition of the band around the DV boundary (Figure 1—figure supplement 2) and the region of the blade that excludes this region (Figure 1D) were not defined manually using the FIJI scripts of TissueMiner. Rather, we defined them after the TissueMiner databases were generated, using the position of cells relative to the DV boundary in the last frame. Using the plots of the radial elongation on each disc, we estimated the width of the stripe region and then filtered for cells included in and excluded from this region. Once these cells in the last frame were identified, we used the UserRoiTracking.R of TissueMiner to backtrack these two regions, producing a list of all cells belonging to this lineage that are traceable forward and backward in the movie.
For the static RNAi data, we used the average width of the band around the DV boundary from timelapse data as an estimate and excluded this region from the static images to analyze the radial gradients in cell area and elongation (Figures 5 and 7).
Plotting area and elongation versus distance to the center of elongation
Request a detailed protocolFor the timelapse data, we defined the radial elongation center (described above) for each movie and then calculated the distance away from this center for each cell. After excluding the band of cells around the DV boundary, we binned cells by radius by rounding to the nearest . We also binned the movie across five equal time windows (), excluding the adaption phase. Average cell area and area-weighted average of cell elongation were calculated for each of these bins in each time group for each movie. Including data from all five movies, we report the global average and its standard deviation (Figure 1E–F). For Figure 3, we averaged over the last . For the band of cells around the DV boundary (Figure 1—figure supplement 2), we binned cells in , rather than in radius, and report the gradient along the axis, defined as the DV boundary.
The static RNAi data were analyzed similarly. We report the global average of all discs in the genotype and the standard deviation (Figures 5D,I and 7D). The number of discs analyzed in each genotype is listed in the figure legend.
Regional analysis of isotropic tissue deformation
Request a detailed protocolSpatial maps of cellular contributions to isotropic tissue deformation
Request a detailed protocolTo analyze the pattern of isotropic deformation, we locally averaged cell behavior by dividing the tissue into a grid centered upon the crossing of the AP and DV boundaries. First, we determined the average position of cells in the AP and DV boundaries to generate a common frame of reference across all movies. Second, we divided the tissue visible in the last frame into a grid centered on these compartment boundary positions, with grid size = . Grid boxes that were incompletely filled (less than 33% of the area of the box filled by cells) were discarded to eliminate noise along the tissue border. Third, each grid box was considered an 'ROI' and then tracked through the entire movie using TissueMiner’s ROI tracking code. Fourth, the rate of deformation by each type of cellular contribution was calculated as a moving average (kernel ) for each grid position for each timestep. Last, we averaged these rates over all time points post-adaption period and plotted in space (Figure 1H,I,K,L).
Plotting the radial profile
Request a detailed protocolTo show tissue deformation as a function of distance to the center, we first calculated the distance to the center of symmetry for each cell. Second, we divided the tissue visible in the last frame into radial bins, rather than a grid, by rounding the radial position of each cell to its nearest . Third, as we did for the grid, we defined each radial bin to be an 'ROI' and then tracked the region forward and backward using TissueMiner’s ROI tracking code. Fourth, the rate of deformation by each type of cellular contribution was calculated as a moving average (kernel ) for each radial bin ROI at each timestep. Last, we averaged these rates over all time points post-adaption period and plotted this rate as a function of distance to the center. Note that we also show the spatial profiles of tissue growth in the band around the DV boundary in Figure 1—figure supplement 1, but here we binned along , not .
Regional analysis of anisotropic tissue deformation
Request a detailed protocolSpatial maps of cellular contributions to anisotropic tissue deformation
Request a detailed protocolWe previously published patterns of cellular contributions to anisotropic tissue shear from movies 1 to 3 (Dye et al., 2017); however here, we calculated these patterns in a more accurate way and report averages over multiple movies (Figure 2A). Previously, we assigned a grid at each timepoint as in Etournay et al., 2015; Etournay et al., 2016. While this method provides a simple first approximation of the pattern, it is not completely accurate because we are not tracking the box in time but reassigning it at each timepoint; thus, cell movement in/out of the box is not counted. Here, we perform grid box tracking, as described above for the calculation of cellular contributions to isotropic tissue deformation but with grid box size . Further, we present a global average of not just movies 1 to 3, but also the two new movies. We calculated for each movie the rate of tissue deformation by each cellular contribution as a moving average (kernel ) in time. Then, we averaged across all five movies at each timepoint and presented the pattern as a cumulative sum, starting from the end of the adaption phase (first ) and continuing through the end of the timelapse. All correlation terms were added together (see Merkel et al., 2017). The contribution to tissue shear from cell extrusion is very small and thus not shown.
Accumulating cellular contributions to radial tissue shear
Request a detailed protocolWe calculated a moving average (kernel size ) total value for each cellular contribution to radial tissue shear, averaged across the blade (excluding the band around the DV boundary). Then, we accumulated the contributions after the adaption phase (first ) and plotted over time (Figure 2B). In addition to the previously described types of correlations (Merkel et al., 2017), the radial decomposition also involves a correlation between cell area and shear (see Appendix 3). We added all correlation terms together in Figure 2B.
Quantification of circular laser ablation
Request a detailed protocolMeasurement of stress, cell elongation, and cell area
Request a detailed protocolFor the ablation data, we first projected the Z-stack images of the wing discs before and after ablation using our custom surface projection (Dye et al., 2017).
To quantify cell elongation and area, we used the images taken before the cut. Cells were segmented using the FIJI plugin TissueAnalyzer and then processed with TissueMiner to rotate the disc to a common orientation (anterior to the left; dorsal up) and to create a triangle network and a database structure. The area-weighted average of triangle elongation and the average cell area was calculated for those cells included in the center of a circle of a radius () centered at the center of the cut. Varying the size of this region only slightly affects the result: , the distribution becomes more noisy because we are averaging a smaller region with less cells, but regions that are too big () may begin to include cells that span different regions of the wing (i.e. the band around the DV boundary). In Figure 4C and Figure 4—figure supplement 1B,E, we plot cell elongation projected onto the stress axis.
To measure the shape of the tissue after ablation, we fit ellipses to the inner and outer piece left by the cut using the projected after-cut images. These images were first rotated using the transformation performed on the before-cut images to orient the tissue. We used Ilastik (v. 1.2.2) (Berg et al., 2019) to segment the cut region: we trained the classifier using all the data from a single day’s acquisition, delineating three regions: cells, membranes, and dark regions. Using the trained classifier, we then generated a thresholded binary image of the cut tissue’s shape and cropped the image around the cut. We then fit inner and outer ellipses to these images.
Grouping ablations into regions of the wing pouch
Request a detailed protocolTo quantify the position of the cut in the wing pouch, the immunofluorescence post-cut images were projected using maximum intensity. In FIJI, we manually measured (in ) the distance from the center of the cut to the DV boundary and AP boundary. We noted that fixation causes the tissue to shrink. We estimated the extent of shrinkage using a subset of the samples in which the compartment boundaries were visible in the Z-stack taken of the live sample before the cut. We used FIJI to measure distances from the center of the cut to the compartment boundaries in both the pre-fix live images and in the post-fix immunofluorescence images for this subset. We measure a discrepancy indicating that the tissue shrinks by ∼15%. Thus, we compensate for this shrinkage when measuring distances to the compartment boundaries in our analysis. We also noted the compartment (dorsal/ventral/anterior/posterior) in post-fix images and added a sign to the distances to the boundaries (ventral and posterior getting negative 'distances') to create the spatial map shown in Figure 4B.
To classify the cuts by position, we first estimated the width of the band around the DV boundary from the timelapse data to be (centered on the DV boundary). We then classified a cut as in the DV boundary if its cut boundary extended no more than (20% of the radius) outside this horizontal stripe region. Likewise, cuts were classified as outside this region if it did not extend more than into the horizontal stripe. We excluded all cuts that appear to straddle the border between these regions (gray in Figure 4B). We also excluded those cuts with centers lying within of the AP boundary, in case material properties at this boundary region are also different.
ESCA: determination of stress from response to circular ablation
Request a detailed protocolTo obtain the normalized shear stress , normalized isotropic stress and the ratio of elastic constant , we fit ellipses to the inner and outer tissue outlines simultaneously. These three parameters determine the large and small semi-axes of the two ellipses. Other fitting parameters are the center positions of the two ellipses and the angles of major axes. Stress is calculated by considering how the ellipses were generated by a circular laser ablation, as described in detail in the Appendix 2. When calculating the deformation from initial to final positions, we considered that the cells that were directly hit by the laser would die and thus not contribute to the evolution of the tissue after the cut. To account for this loss of cells, we first calculated average cell area for those cells that would be hit by the laser. Then, we adjusted the radius of the initial position by this amount.
Fits were performed on all of the laser ablations in the same region (either within or outside of the band around the DV boundary) for a range of fixed values of the ratio . The optimal value was considered to be the one that minimizes the sum of fit residuals (Figure 4—figure supplement 1C,F). We defined a threshold of fit residual, as shown in Figure 4—figure supplement 1A, to eliminate cuts with non-elliptical outlines. The uncertainty of was estimated by bootstraping a subset of 7 cuts in each region 100 times; we report the standard deviation of the results.
Appendix 1
Continuum model of the wing disc epithelium
1 Model definition
We use a previously developed continuum model to describe tissue mechanics (Etournay et al., 2015; Popović et al., 2017). Tissue shear flow consists of convected co-rotational change of cell elongation and a contribution from cell rearrangements
Tissue shear stress consists of an elastic part and a contribution from nematic cell polarity
Shear flow due to cell rearrangements consists of term accounting for cell shape relaxation and a contribution due to nematic cell polarity
Force balance couples shear stress and area pressure
which will be reflected in cell area. We use a constitutive relation
where represents cell area and corresponds to value of at .
We consider a simple model in which the tissue is a radially symmetric disc, see Figure 3C–E of the main text. We use polar coordinate system , where the radial force balance reads
Furthermore, in the model we consider a constrained tissue with negligible radial shear flow , consistent with experimental observations, see Figure 2 of the main text.
2 Cell area profiles follow from cell elongation profiles
Empirically, the radial cell elongation profile can be described by a power law
with values for and obtained by fitting the data on the last of the experiments, see red line in Appendix 1—figure 1A. We find and by fitting the data from the first . Data is averaged over five experiments, and errorbars represent a standard deviation among the mean values in each experiments.
At this point, it is useful to estimate the spatial profile of shear due to cell rearrangements . We calculate the difference between late and early cell elongation profiles, as shown in Appendix 1—figure 1B. Since we find with no significant variation in .
Now, we find an expression for cell area using Equations A2, A3, A5, A6 and A7
where
and is an integration constant. We fit the Equation A8 to the experimental data as shown in Figure 3G of the main text, and we find parameter values
where the reported uncertainties represent square roots of the fit covariance matrix and the parameter in Equation A12 was constrained to be positive, motivated by the fact that we do not find (see Equation A18 below).
3 Laser ablation experiments provide an estimate of model parameters
We write Equation A2 in the form
We fit this equation to the experimentally obtained shear stress as a function of cell elongation in the laser ablation experiments. We find
Using the estimate for we find
This value is consistent with the one found for the cell elongation relaxation time-scale in the pupal wing (Etournay et al., 2015).
In Equation A3, we can now show that cell elongation profile in the wing pouch reflects the profile of nematic cell polarity. Namely, in most of the tissue, except in the vicinity of the DV boundary, and thus
Comparison of laser ablation and area profile fits
It is interesting to notice that from Equation A11 and the value for obtained from the laser ablation experiments (main text and Materials and methods), we can estimate
This value is significantly higher than the one obtained by a direct fit to the laser ablation data (Equation A14). This apparent discrepancy stems in part from the fact that the elastic constant in Equation A5 is a non-linear elastic modulus. It can only be related to the one used in the laser ablation method by linearizing it around a typical value of in the data. When this value is not equal to the linearized elastic constant will differ from the original in Equation A5. Furthermore, Equation A5 does not take into account any adaptation of cell area to pressure or nematic cell polarity, which could influence the measured value when fitting the cell area profile. However interesting, these effects are beyond the scope of our current work.
4 Mechanosensitivity leads to spontaneous cell polarity by self-organization
a Homogeneous tissue
We propose a dynamical equation of the polarity tensor
where the parameter represents the mechanosensitive response of nematic cell polarity to shear stress. We use Equations A1, A2, A3 and A19 to obtain a dynamic equation for shear stress
The system of linear Equations A19 and A20 becomes unstable when
This instability will result in a finite value of nematic cell polarity in the tissue. To account for this polarized state we have to include a stabilising higher order term in Equation A19
In the steady state, the magnitude squared of the nematic cell polarity is
Here, we have introduced the strength of polarity , corresponding to the magnitude of polarity in a homogeneous system. Note that the definition of allows for negative values, which would correspond to no spontaneous polarization in the tissue.
b Polarization of a radially symmetric tissue
To account for spatial variations of nematic cell polarity, we introduce a Laplacian term in Equation A22
To fit the radial profile of observed in time-lapse experiments we account for non-vanishing as well. Consistent with the observation that during the experiments does not vary significantly in time (Figure 2B main text) we use approximation . Then satisfies
where the left hand side is the radial component of the Laplacian of in the polar coordinate system. From the solution of this equation, the profile is obtained as . This quantity is fitted to the experimentally measured tangential cell elongation profile by optimizing four tissue parameters, as well as two boundary conditions and required to solve Equation A25. Here is the smallest radius at which we measure the tangential cell elongation profile in Figure 6C of the main text. Parameters obtained by the fit are reported in Table 1 of the main text, and the fitted profile is shown in Figure 6 of the main text. Parameter uncertainties were estimated as the 10th and 90th percentile intervals obtained by fitting uniformly sampled profiles of from the interval and , see Figure 3 of the main text, with 101 realisations. Obtained uncertainty intervals are reported in the Table 1 of the main text.
c Predictions for experiments
Spontaneous cell polarity in our model depends on mechanosensitive feedback characterized by the feedback strength . The characteristic polarity strength grows monotonically with . A prediction from our model is therefore that when is reduced, cell elongation will be reduced because it is mediated by cell polarity. In order to test this prediction, we performed MyoVI knockdown experiments, which indeed show a decrease in cell elongation (Figure 7).
Our qualitative prediction is rather robust and does not depend on details of the system, such as the influence of compartment boundaries that disrupt the radial symmetry. However, boundary conditions, which we do not know, could have a significant influence on the system. Therefore, a full quantitative prediction of experiments is difficult. To show the self-consistency of our prediction, we test if the value of that we use to describe the data is consistent with the actual cell elongations observed. Indeed we have , see Table 1 which predicts a typical strength of cell elongation , consistent with experimentally measured values, see Figure 6C of the main text. This suggests that our qualitative prediction is robust given the uncertainties in the system.
Appendix 2
Theory of circular laser ablations
We developed a method to infer the stresses in an epithelial tissue from the tissue boundary shape after a circular ablation. Here, we first briefly review linear elasticity in polar coordinates. Then, we derive equations relating the shape of the inner and outer boundaries obtained by the ablation of an infinite elastic sheet in the small deformation regime. The solution to this problem is by no means new, see for example (Muskhelishvili, 1977). Here, we present our derivation of the solution. The application of the method and the obtained parameter estimates are described in the Materials and methods section.
1 Linear elasticity in polar coordinates
a Displacement gradient
The deformation of an elastic object is described by the displacement vector of each point on the sheet
where is the intial position of a point on the elastic sheet and is the position of that point after the deformation. Shape changes are described by the symmetric traceless gradient of displacement
which can be decomposed into the isotropic part and the shear part . Note that we use summation convention over repeated indices. In polar coordinates displacement gradient is given by
b Linear elastic constitutive relation
We are considering the linear elastic constitutive relation
where is the traceless symmetric component of shear stress, with shear elastic modulus and bulk elastic modulus . Here, we allow for the presence of active stresses due to a nematic field , corresponding to the nematic cell polarity in the wing disc. However, a constant active stress does not influence our results in the small deformation limit, since it simply redefines the reference configuration. We can therefore write this constitutive equation as
c Force balance and compatibilty: Airy stress function
Force balance imposes
In two dimensions, force balance is satisfied by all stress fields whose components are derived from the Airy stress function
which can be expressed in polar coordinates
Now, a compatibility condition, which follows from a requirement of continuity of the displacement field, can be written as
The general solution of this equation in polar coordinates, also called the Michell solution, is given by
Inferring the stresses from laser ablation experiments consists of solving two elastic problems using the Michell solution as we now show.
2 Circular ablation of an elastic sheet under stress
Laser ablation splits the elastic sheet into inner and outer pieces. We now show that the boundary shapes of the two pieces are ellipses by solving the corresponding elastic problems.
a Inner piece
We consider a piece of elastic sheet stretched by stress , , . After a circular laser ablation, the inner piece is under no stress. To infer the original stresses in the sheet, we determine the shape of the sheet boundary assuming the stresses are known. Then, given the boundary shape, we will be able to infer the stresses.
Stress in a sheet that is in mechanical equilibrium, under no external body force, is uniform throughout the sheet. Therefore, using Equation B3 we can integrate Equation B7 to find the components of the displacement vector
where we have introduced
To relate the displacement to the shape of the inner piece boundary, we have to solve Equation B1 for , in the configuration illustrated in Appendix 2—figure 1
Using the fact that
we obtain
To simplify the notation, we introduce
and we write the solution in the form
We denote semiaxes of this ellipse along the and axes, as defined in Appendix 2—figure 1, by and , respectively. We find that
b Outer piece
As in the case of the inner piece, we need to find the boundary shape as a function of the original stresses in the sheet. In this case, the normal stress at the ablation boundary is zero, but the sheet is not stress-free. To simplify the problem, we consider an initially stress-free sheet with a circular hole to which boundary stresses are applied. For small deformations, where the elastic problem is linear and solutions of the problems commute, the boundary shape is equivalent to the one created by the laser ablation of a sheet under stress .
Now, we have to solve the full elastic problem, since stresses are not uniform. We keep the terms of the Airy stress function , which obey nematic symmetry of the problem and relax sufficiently fast at infinity. We also assume that the polar axis is oriented along one of the principal axes of stress tensor, such that the off-diagonal component of stress vanishes. We obtain
and for stress components, using Equation B10
The stress boundary condition reads
allows us to determine coefficients in Equation B29
Now we can calculate and then integrate the deformation gradient components to obtain
Finally, as for the inner cut, we have to solve
for . We express the solution using the angle as a parameter
Note the ratio of elastic constants appears as a parameter of shape.
Now we show that the shape defined by Equations B40 and B41 is an ellipse. To this end, we first define
This allows us to write Equations B40 and B41 as
Now we show that is an ellpise. An ellipse equation can be written as
where and are semiaxes of the ellipse oriented along and axes, respectively. Starting from Equation B46 and expressing from Equation B44 we will show that we can recover Equation B45 for a particular choice of parameters and .
We first calculate
Inserting these identities into Equation B46 we find
Now, we only need to find and for which
for all . The solutions are
or in terms of the original parameters
If, in addition to Equations B23 and B24 we define
we can write
Equations B26, B27, B56 and B57 impose four constraints on three parameters: , , and . When determining these parameters from experimental data we are in general not able to satisfy all four equations simultaneously. Instead, we find the three parameters that best fit the segmented outlines of tissue boundaries, see Materials and methods.
Appendix 3
Correlation between cell elongation and cell movement contributes to the tissue radial shear flow
In this work, we use the previously developed method to define and calculate the cell elongation tensor and tissue shear due to T1 transitions, cell divisions, and extrusions (Etournay et al., 2015; Merkel et al., 2017). However, we apply the method in the polar coordinate system. To make it fully consistent, we have to introduce a new correlation term , which arises in the polar coordinate system when translation and elongation are correlated. To demonstrate the origin and derive an expression for , we consider a tissue that translates over the time interval but does not change otherwise. The tissue is triangulated by connecting the neighboring cell centers, see Merkel et al., 2017. Initially, the angle of a triangle with with respect to a center point is , blue triangle in Appendix 3—figure 1. During the translation by , the radial angle of the triangles changes by a small angle . Therefore, the components of the triangle elongation tensor in the polar coordinate system change by
However, since the triangles are only translating, there is no shear flow. Therefore, the average radial shear flow also vanishes and an additional term is necessary to account for the change in the average radial cell elongation in a translating tissue
Therefore, the correlation term is given by
In Appendix 3—figure 2, we show the new correlation term in black measured in the pouch of the wing disc, averaged over five time-lapse experiments. For comparison, we show the correlation term stemming from correlated fluctuations of triangle rotation and elongation in blue and the correlation term stemming from correlated fluctuations of triangle area and elongation in green. For definitions and discussion of the other two correlation terms see Etournay et al., 2015; Merkel et al., 2017. We find that the contribution of the new correlation term to the radial shear flow is much smaller that the other two correlation terms. Therefore, the contribution to radial tissue shear flow of the correlation between translation and triangle elongation is negligible in the wing disc.
Data availability
We have made all data analyzed during this study available. Data for Figures 1H–M, 2,4,5, and 7 are provided as source data files. The data on cell area and elongation in Figure 1A–F, 3F,G, and 6C are too large to be submitted here and are available on Dryad (https://doi.org/10.5061/dryad.jsxksn06b).
-
Dryad Digital RepositoryData from: Self-organized patterning of cell morphology via mechanosensitive feedback.https://doi.org/10.5061/dryad.jsxksn06b
References
-
Mutations in the cadherin superfamily member gene dachsous cause a tissue polarity phenotype by altering frizzled signalingDevelopment 125:959–968.
-
Segmentation and quantitative analysis of epithelial tissuesMethods in Molecular Biology 1478:227–239.https://doi.org/10.1007/978-1-4939-6371-3_13
-
ilastik: interactive machine learning for (bio)image analysisNature Methods 16:1226–1232.https://doi.org/10.1038/s41592-019-0582-9
-
Pattern formation in active fluidsPhysical Review Letters 106:028103.https://doi.org/10.1103/PhysRevLett.106.028103
-
Mechanical state, material properties and continuous description of an epithelial tissueJournal of The Royal Society Interface 9:2614–2623.https://doi.org/10.1098/rsif.2012.0263
-
Planar cell polarity in development and diseaseNature Reviews Molecular Cell Biology 18:375–388.https://doi.org/10.1038/nrm.2017.11
-
Apical cell shape changes during Drosophila imaginal leg disc elongation: a novel morphogenetic mechanismDevelopment 111:23–33.
-
The Frizzled Planar Cell Polarity signaling pathway controls Drosophila wing topographyDevelopmental Biology 317:354–367.https://doi.org/10.1016/j.ydbio.2008.02.041
-
Cell biology of planar polarity transmission in the Drosophila wingMechanisms of Development 120:1257–1264.https://doi.org/10.1016/j.mod.2003.07.002
-
A genetic analysis of the determination of cuticular polarity during development in Drosophila melanogasterJournal of Embryology and Experimental Morphology 68:37–57.
-
Turing's next steps: the mechanochemical basis of morphogenesisNature Reviews Molecular Cell Biology 12:392–398.https://doi.org/10.1038/nrm3120
-
Hydrodynamic theory of active matterReports on Progress in Physics 81:076601.https://doi.org/10.1088/1361-6633/aab6bb
-
Building an ommatidium one cell at a timeDevelopmental Dynamics 241:136–149.https://doi.org/10.1002/dvdy.23707
-
Front-Rear polarization by mechanical cues: from single cells to tissuesTrends in Cell Biology 26:420–433.https://doi.org/10.1016/j.tcb.2016.02.002
-
Hydrodynamics of soft active matterReviews of Modern Physics 85:1143–1189.https://doi.org/10.1103/RevModPhys.85.1143
-
Roles of the cytoskeleton, cell adhesion and rho signalling in mechanosensing and mechanotransductionJournal of Biochemistry 161:245–254.https://doi.org/10.1093/jb/mvw082
-
Active dynamics of tissue shear flowNew Journal of Physics 19:033006.https://doi.org/10.1088/1367-2630/aa5756
-
Active gel segment behaving as an active particlePhysical Review E 100:062403.https://doi.org/10.1103/PhysRevE.100.062403
-
Fiji: an open-source platform for biological-image analysisNature Methods 9:676–682.https://doi.org/10.1038/nmeth.2019
-
Laser ablation to probe the epithelial mechanics in DrosophilaMethods in Molecular Biology 1478:241–251.https://doi.org/10.1007/978-1-4939-6371-3_14
Article and author information
Author details
Funding
Max-Planck-Gesellschaft
- Natalie A Dye
- Marko Popovic
- K Venkatesan Iyer
- Jana Fuhrmann
- Romina Piscitello-Gómez
- Suzanne Eaton
- Frank Jülicher
Deutsche Forschungsgemeinschaft (EA4/10-1)
- Natalie A Dye
- K Venkatesan Iyer
- Suzanne Eaton
Swiss National Science Foundation (200021-165509)
- Marko Popovic
Simons Foundation (454953)
- Marko Popovic
Deutsche Forschungsgemeinschaft (EA4/10-2)
- Natalie A Dye
- K Venkatesan Iyer
- Suzanne Eaton
ELBE PhD program
- Romina Piscitello-Gómez
Deutsche Krebshilfe (MSNZ Dresden)
- Natalie A Dye
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Stephan Grill for the microscope used to perform circular laser ablations. We thank the Light Microscopy Facility of the MPI-CBG for assistance with all other imaging experiments. Benoit Lombardot, formerly of the Image Analysis Facility of the MPI-CBG, provided us with a FIJI script to perform the projection of the apical surface, originally designed by Dagmar Kainmueller. Franz Gruber generated the pk30,ecadGFP fly line. We thank Christian Dahmann, Charlie Duclut, Kinneret Keren, Yonit Maroudas-Sacks, Alphee Michelot, and Ioannis Nellas for critical review of the manuscript. Lastly, we dedicate this manuscript to our co-author, Suzanne Eaton, who passed away tragically toward the conclusion of this work. This work was supported by funding from the Max-Planck-Gesellschaft (NAD, KVI, JFF, RPG, SE, FJ), Deutsche Forschungsgemeinschaft (EA4/10-1, EA4/10-2; NAD, KVI, SE), Swiss National Science Foundation (Grant #200021–165509; MP), the Simons Foundation (Grant #454953 to Matthieu Wyart; MP), the Deutsche Krebshilfe/MSNZ Dresden (NAD), and the ELBE PhD program (RPG). Funding sources were not involved in the study design, data collection, data interpretation, or the decision to submit the work for publication.
Copyright
© 2021, Dye et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 5,059
- views
-
- 824
- downloads
-
- 40
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Developmental Biology
- Genetics and Genomics
O-GlcNAcylation is an essential intracellular protein modification mediated by O-GlcNAc transferase (OGT) and O-GlcNAcase (OGA). Recently, missense mutations in OGT have been linked to intellectual disability, indicating that this modification is important for the development and functioning of the nervous system. However, the processes that are most sensitive to perturbations in O-GlcNAcylation remain to be identified. Here, we uncover quantifiable phenotypes in the fruit fly Drosophila melanogaster carrying a patient-derived OGT mutation in the catalytic domain. Hypo-O-GlcNAcylation leads to defects in synaptogenesis and reduced sleep stability. Both these phenotypes can be partially rescued by genetically or chemically targeting OGA, suggesting that a balance of OGT/OGA activity is required for normal neuronal development and function.
-
- Developmental Biology
New research shows that the neural circuit responsible for stabilising gaze can develop in the absence of motor neurons, contrary to a long-standing model in the field.