Epithelial colonies in vitro elongate through collective effects

Epithelial tissues of the developing embryos elongate by different mechanisms, such as neighbor exchange, cell elongation, and oriented cell division. Since autonomous tissue self-organization is influenced by external cues such as morphogen gradients or neighboring tissues, it is difficult to distinguish intrinsic from directed tissue behavior. The mesoscopic processes leading to the different mechanisms remain elusive. Here, we study the spontaneous elongation behavior of spreading circular epithelial colonies in vitro. By quantifying deformation kinematics at multiple scales, we report that global elongation happens primarily due to cell elongations, and its direction correlates with the anisotropy of the average cell elongation. By imposing an external time-periodic stretch, the axis of this global symmetry breaking can be modified and elongation occurs primarily due to orientated neighbor exchange. These different behaviors are confirmed using a vertex model for collective cell behavior, providing a framework for understanding autonomous tissue elongation and its origins.


Introduction
Tissue elongation is a central morphogenetic event occurring in many organisms during development (Lecuit and Lenne, 2007;Guillot and Lecuit, 2013), such as Drosophila or Caenorhabditis elegans. The tissue is transformed both in terms of area and shape. Such transformation takes place within typically hour timescale with or without cell division. During this process, symmetry of cells and tissues is broken by different mechanisms, such as neighbor exchange (Rauzi et al., 2008;Rauzi et al., 2010), cell elongation (Ciarletta et al., 2009;Vuong-Brender, 2016), and oriented cell division (Campinho et al., 2013). Rearrangement of neighboring cells or T1 transitions is essential in the germ band extension of Drosophila (Rauzi et al., 2008;Rauzi et al., 2010), allowing a group of cells to change their position by intercalation, eventually leading to tissue elongation. Cell deformation drives the threefold elongation process in C. elegans (Ciarletta et al., 2009;Vuong-Brender, 2016) while keeping the number of cells and their positions fixed. Finally, epithelial spreading during zebrafish epiboly is promoted by oriented cell divisions as a mechanism to limit tension (Campinho et al., 2013). Those mechanisms can act alone or in combination as in Drosophila pupal wing elongation (Etournay et al., 2015). While the phenomenon is known to involve remodeling of adherens junctions (Rauzi et al., 2010) and acto-myosin (He et al., 2010;Rauzi et al., 2010) at the molecular level, mesoscopic mechanisms leading to distinct morphogenesis processes are poorly understood. This is partly because inputs from morphogen gradients (Gilmour et al., 2017) or from neighboring tissues (Zhang et al., 2011;Etournay et al., 2015) can affect tissue autonomous selforganization in vivo. For example, changes in tissue shape can be influenced by neighboring tissues such as the cuticle and the blade in the Drosophila pupal wing elongation (Etournay et al., 2015), the coordination between amnioserosa and epidermis in dorsal closure (Hayes and Solon, 2017), and the muscle layer in gut development (Shyer et al., 2013) or in C. elegans morphogenesis (Zhang et al., 2011). Since in vivo, epithelial tissues are surrounded by other tissues and the cellular dynamics leading to elongation can result from interactions between tissues and boundary conditions, it is therefore difficult to disentangle cell intrinsic from externally mediated behaviors. In this context, it appears important to characterize elongation in an in vitro system where the epithelial tissue undergoes shape transition autonomously.
Here, we use in vitro experiments and numerical simulations to characterize the spontaneous behavior of a growing cell colony in vitro. We designed an assay to study the spontaneous elongation of a tissue that is not subjected to external orienting input, we studied the appearance of the symmetry breaking, and the effect that external forces have in this process. We show that in vitro tissue elongation arises from anisotropy in the average cell elongation. This anisotropy sets the direction along which boundary cells migrate radially outwards resulting in a non-isotropic elongation that arises primarily through cell elongation. For colonies submitted to a time periodic uniaxial stretch, the axis of global symmetry breaking can be imposed by external force, and tissue elongation arises through oriented neighbor exchange. Emergence of radially migrating cells and the interplay between cell elongation and cell rearrangements are confirmed by numerical simulations based on a vertex model. Our results suggest that spontaneous shape deformation is related to the mean orientation of the nematic cell elongation field in the absence of any external input. This provides a framework to explain autonomous tissue elongation and how contributions from different mesoscopic mechanisms can be modulated by external forces.

Isotropic colonies elongate in a non-isotropic manner
To study the spontaneous tissue deformation arising during epithelial growth, we designed an in vitro assay to track symmetry breaking, both spontaneous and driven by external force. We prepared isotropic colonies of Madin Darby Canine Kidney (MDCK) cells, which assume features of epithelial cells in vivo (Reinsch and Karsenti, 1994;Adams et al., 1998;Reffay et al., 2014), such as adherens junctions (Adams et al., 1998), cytoskeletal components, and the Rho signaling pathway regulating cell shapes and dynamics (Reffay et al., 2014;Fodor et al., 2018). The initial size and shape of the colonies were controlled by plating cells in microfabricated circular stencils (Ostuni et al., 2000). When cells reached confluency, the stencil was removed at time t 0 . Cell dynamics was followed over time by phase contrast (Video 1) or fluorescence microscopy with strains labeled with GFP cadherin (Figure 1b), that allowed to observe the behavior of individual cells. We observed that large colonies (750 mm in diameter) expanded isotropically ( Figure 1-figure supplement 1). In contrast, colonies of 250 mm in diameter (Figure 1b), the typical coherence length of such epithelial tissues (Doxzen et al., 2013), expanded in a non-isotropic manner ( Figure 1c). We then further characterized the process of symmetry breaking. In order to compare elongations in each experiment, we quantified the breaking of symmetry by ellipse-fitting the colony shape. Shape change analysis was quantified by a nematic shape elongation tensor Q. It has two independent components defined as Q xx = ½ ln(a/b)cos(2 . ( À a)) and Q xy = ½ ln(a/b)sin (2 . ( À a)), where a corresponds to the major axis, b to the minor axis, to the orientation of the major axis of the ellipse and a = q(t final ) ( Figure 1d). As can be seen in Figure 1e, MDCK colonies elongated persistently along the main axis of elongation (Q xx > 0 and Q xy » 0) for 6 hr ( Figure 1e). In addition, we explored if other epithelial cell lines would behave in a similar manner. Circular epithelial colonies of human epithelial colorectal adenocarcinoma cells (Caco2) and human mammary epithelial cells (MCF-10A) also elongated along the main axis of elongation and by the same magnitude that MDCK cells (Figure 1-figure supplement 2). We note that elongation observed during this time for the three epithelial cell lines was similar in magnitude to tissue elongation observed during in vivo morphogenesis, for instance in the wing blade in Drosophila (Etournay et al., 2015). Moreover, the elongation direction ( final = (t = 6 h)) converges to a constant value within 2 hr after t 0 ( Figure 1f). Altogether, large circular epithelial colonies (750 mm in diameter) expand isotropically, whereas small colonies (250 mm in diameter) expand in a anisotropic manner and shape symmetry breaking takes place within the first 2 hr. As a result, we focus here on these first 2 hr during which the elongation axis is established.
Cyclic uniaxial stretching rectifies symmetry breaking It has been previously described for C. elegans embryo elongation (Zhang et al., 2011) and in other organisms (Zhang and Labouesse, 2012) that time periodic stretch can play a role in morphogenesis. Motivated by these observations, we explored whether oscillatory external forces could have an impact on the direction of elongation. We designed an experimental setup where elongating colonies were submitted to cyclic uniaxial stretching ( Figure 2a and Video 2). Mechanical cycles of contraction-relaxation can range from 1 s in C. elegans epithelial elongation (Zhang et al., 2011) up to 200 s in dorsal closure (Solon et al., 2009). So, we explored frequencies and extensions around physiological values (Zhang and Labouesse, 2012). We selected three different cycle durations (20, 60, and 120 s period) and three different stretching conditions (5%, 10% and 15% strain). The stretch was applied to a silicon membrane and was transmitted to the colony. We then fitted the colonies shapes with ellipses at successive time and quantified Q xx and Q xy with respect to the angle of uniaxial stretching (set as x-axis, a = 0). Figure 2-figure supplement 1 shows the value of the components of the tensor Q along time for the different strains and periods tested. Among different conditions, we observed colony elongation along the direction imposed by the external strain when we stretched cyclically with 60 s timescale and 5% strain (Figure 2b and Video 3). The overall elongation of colonies under cyclic uniaxial stretching was similar to the spontaneous elongation in the absence of externally applied uniaxial stretching during the first 2 hr (Figure 2c). Also, the magnitude of the shape elongation tensor Q under cyclic uniaxial stretching was comparable to the  elongation is quantified by ellipse fitting and Q xx and Q xy measurement referred to the elongation axis (a = q(t final )). (e) Q xx (left y axis) and Q xy (right y axis) during 360 min of colony expansion. Mean value ± standard error of the mean, n = 4 colonies from N = 4 independent experiments. (f) Cosine of two times the angle difference between the instantaneous main direction of the colony (q(t)) and the main direction of the colony at 360 min (q(t final )). Colonies set the elongation direction within the first 120 min. Mean value ± standard error of the mean, n = 4 colonies from N = 4 independent experiments. The online version of this article includes the following source data and figure supplement(s) for figure 1: Source data 1. The Q xx (t)-Q xx (0) and Q xy (t)-Q xy (0) of individual colonies used in panel (e) and the raw values of cos(2Á(q(t)-q(t final ))) used in panel (f) of Figure 1.    spontaneous elongation of colony when stretch was not applied, but elongation was oriented in the direction of externally applied uniaxial cyclic stretching ( Figure 2d). Therefore, application of an external cyclic force can rectify symmetry breaking and set the direction of tissue elongation.

Collective effects are essential for rectification
To get further insight into the collective nature of the rectification of tissue elongation, we probed the roles of adhesion between cells. First, we stretched single MDCK cells, individually plated. We observed that cells oriented perpendicularly to the externally applied uniaxial cyclic stretching (Figure 3a and b) as previously reported for fibroblasts (Faust et al.,  (c) Colony elongation (Q xx ) of control colonies along the elongation axis, control colonies in the laboratory framework (non-aligned, a = 0) and colonies under cyclic uniaxial stretching in the laboratory framework (uniaxial stretching, a = 0). Box Plot between 25th and 75th percentile, being the line in the box the median value, whiskers and outliers (dots) are obtained following Tukey's method, N control = 11 independent experiments, n = 25 colonies and N stretching = 9, n = 20 colonies. Mann-Whitney test control vs control aligned p=0.0003, control vs stretching p=0.0281 and control aligned vs stretching p=0.3319. (d) Q xx (left y axis) and Q xy (right y axis) during 120 min of colony expansion for control colonies and colonies under cyclic uniaxial stretching. Mean value ± standard error of the mean, N control = 8, n > 14 colonies and N stretching = 9, n = 20 colonies. The online version of this article includes the following source data and figure supplement(s) for figure 2: Source data 1. Raw data used for the panel (c) of Figure 2 and Q xx (t)-Q xx (0) and Q xy (t)-Q xy (0) for the individual colonies used for panel (d) of Figure 2.  2011). Then we blocked cell-cell junction in circular colonies prior stretching. Briefly, we incubated confluent colonies in medium containing 5 mM EDTA and 10 mg/ml anti-E-cadherin blocking antibody that targeted the extracellular domain of E-cadherin for 30 min (Harris et al., 2014). Then, medium was replaced by normal medium and the evolution of colonies with and without stretch was followed ( Figure 3c). In the absence of externally applied uniaxial cyclic stretching, colonies treated with anti E-cadherin antibody expanded more than control colonies. Moreover, this expansion was still along one preferential direction (Figure 3d). Under cyclic uniaxial stretching, elongation was also non-isotropic and along the direction perpendicular to the cyclic stretching direction, in contrast to the parallel elongation observed when cell-cell contacts were intact (Figure 2b-d). This supports the collective nature of colony elongation. It is worth noting that cells inside the colony exhibited a decrease in their mean velocity ( Figure 3e) and a large recruitment of myosin within cells similar to reinforcements (Riveline et al., 2001) in stretching conditions, as shown by the appearance of stress fibers ( Figure 3f). However, this effect did not appear to affect the overall elongation rate. Altogether these data suggest that the asymmetric expansion of colonies in the direction imposed by cyclic uniaxial stretching is generally associated to a collective effect.

Fingers and symmetry breaking
We next sought to identify the source of symmetry breaking in both conditions, with and without application of cyclic uniaxial stretching. It has previously been reported that in MDCK monolayers, cells can migrate tangentially to the monolayer boundary when confined (Doxzen et al., 2013), or perpendicular to the boundary in the form of multicellular groups or fingers during monolayer expansion (Reffay et al., 2011;Reffay et al., 2014). During spontaneous elongation of MDCK colonies in the absence of externally applied cyclic stretching (Video 4), we observed that boundary cells tend to move either perpendicularly or tangentially to the colony boundary ( Figure 4a). In most of the experiments, an acto-myosin cable, similar to compartmentalization boundaries in vivo (Monier et al., 2010;Calzolari et al., 2014), was observed in the outer boundary of the colonies at stencil removal (see Figure 4-figure supplement 1). When this supra-cellular structure is intact, cells at the periphery are reported to undergo clockwise and counter clockwise rotations (Doxzen et al., 2013). In contrast, when a local interruption of this cable appeared, the cell at the vicinity could extend a lamellipodia and move away and radially from the center of the colony (Figure 4-figure supplement 2a). Apparently, a local defect in the cable could promote the local protrusion of a cell in the direction normal to the edge as shown in laser ablation experiments previously (Reffay et al., 2014). Several local defects could appear within the same colony, thus providing the opportunity for cells in the vicinity to protrude outwards. This cell has been termed leader cell (Reffay et al., 2011) and the collection of cells protruding from the circular colonies along this cell can be identified as the finger-like structures already reported for MDCK monolayers (Reffay et al., 2011;Reffay et al., 2014).
We performed cell tracking and observed that, on average, these protruding cells are faster than other boundary cells ( Figure 4a and Figure 4-figure supplement 2b and c). They are characterized also by radial and directional migrations, in contrast to tangential motion observed in the other cells of the outer region of the colony (Figure 4-figure supplement 2d and e). In general, the motion of these so-called leader cells was directionally persistent and on average the shape of the whole colony followed the same overall direction (Figure 4b and c). To correlate colony elongation with leader cell orientation, we analyzed the evolution of a larger number of colonies for 2 hr after stencil removal (Video 4). We quantified the breaking of symmetry by fitting an ellipse to the shape of each colony. We then tracked the positions where finger-like structures were appearing, as well as the direction and distance performed by each of them. We could observe that the elongation direction of the whole colony correlated on average with the direction of the leader cell migration and associated finger (Figure 4c).
We then measured the position and displacement for each finger in control colonies and colonies under cyclic uniaxial stretching ( Figure 4d). We observed that, when growing perpendicular to the direction of force application, finger cells performed shorter displacements than when growing parallel to it. In the absence of externally applied cyclic uniaxial stretching, fingers grew a similar amount as when growing parallel the direction of applied uniaxial cyclic stretching and no bias was observed vis-à -vis the nucleation position ( Figure 4d and Figure 4-figure supplement 2f). To further explore this effect, we grew MDCK monolayers with straight boundaries either parallel or perpendicular to the external force. Then, we let fingers appear and grow for 2 hr before applying cyclic uniaxial stretching ( Figure 4e). When fingers were growing perpendicular to the stretching direction, they shrank upon application of cyclic uniaxial stretching; in contrast, fingers further elongated when parallel to the direction of uniaxial cyclic stretching. Altogether, this suggests that direction of finger-like structures correlates with elongation direction, and that external stretching affects the dynamics of finger growth.

Collective effects and symmetry breaking
Finger growth correlates with colony elongation. However whether it is a cause or consequence of the symmetry breaking of the shape of the colony remains elusive. We therefore explored the possibility of inducing the growth of fingers and therefore set the direction of elongation of the colonies. Breakage of the acto-myosin cable by laser ablation induces the appearance of leader cells (Reffay et al., 2014). Hence we attempted to trigger the growth of fingers by locally injecting cytochalasin D using a micropipette. The transient injection of this actin polymerization inhibitor was followed by the disruption of the acto-myosin cable (Video 5 and Figure 5-figure supplement 1). However, the cable reformed, and fingers did not appear. This result shows that breakage of the cable alone does not trigger the growth of fingers in our colonies, and suggests that other mechanisms may be involved.
We observed that when a finger moves outward from the colony, cells in the immediate vicinity elongate and seem to reorient their elongation axis toward the finger (Figure 5a). Recent studies have shown that the nematic field of cell elongation and its topological defects could be involved in the growth of bacterial colonies (Doostmohammadi et al., 2016) and in controlling dynamics, death and extrusion of epithelial cells (Kawaguchi et al., 2017;Mueller et al., 2019;Saw et al., 2017). We wondered if the spontaneous elongation of colonies would also be related to the average cell elongation. We followed the evolution of the cell elongation nematic field in different MDCK and MCF 10A colonies during expansion. We first obtained the spatio-temporal cell elongation nematic orientation field f x; y; t ð Þ (see Materials and methods) on the experimental time-lapse images (see Figure 5b, Figure 5-figure supplement 2a-c, Video 6 and Appendix 1C). We could then obtain the angle nematic of the average cell-shape nematic field at t 0 which we compared with final colony orientation colony obtained using the ellipse fitting analysis ( The cell elongation orientation field f x; y; t ð Þ was not homogeneous at t 0 but exhibited a complex pattern with ±½ topological defects (Figure 5b and Appendix 1D). Interestingly, an expression that provides equilibrium orientation of liquid crystals with defects and having one constant Frank free energy, mimics the experimentally observed orientation field f x; y ð Þ very well with just one fitting parameter a (see Appendix 1D for details). Here k i = ±½ and (x i , y i ) are the strength and location, respectively, of the topological defects obtained from the experimental image. Thus, the defect position and  strength can be used to provide an approximate readout for the orientation of the cell-shape nematic field ( Figure 5d). Moreover, the location of finger nucleation seemed to be biased toward the position of topological defects. However, some defect locations were not stable in time and in some cases, the nematic field of cell shapes could only be interpreted in terms of virtual defects outside the colonies, thus suggesting that the mean nematic direction is a better readout for the cellshape nematic field.
On the one hand, finger nucleation seemed to be correlated with colony elongation direction ( Figure 4c). On the other hand, the orientation of tissue elongation correlates with orientation of average cell elongation at t 0 ( Figure 5 and Figure 5-figure supplement 2). This suggests that leader cells moving outward from the colony may not be the cause of symmetry breaking in colony shape, but rather follow from the initial cell shape elongation before stencil removal. Moreover, we found no correlation between breaks of the acto-myosin cable surrounding the colony and the mean nematic direction (Figure 5-figure supplement 3), which suggests that breaks are uniformly distributed along the colony border. We have also shown that breakage of acto-myosin cable after stencil removal, which is associated with leader cell formation (Reffay et al., 2014), did not necessarily induce the growth of fingers in our colonies. Altogether, our results could indicate that the orientation of the mean cell-shape nematic of the colony before stencil removal sets the direction of elongation by triggering the growth of fingers, which appear at the discontinuities of the outer actomyosin cable located along the nematic orientation, while preventing finger growth at discontinuities located in other directions.
Finally, when looking at the evolution of the mean cell elongation nematic field of colonies under uniaxial cyclic stretching, we observed that it did not change over time ( Figure 5-figure supplement 4). The initial mean direction of cell elongation, either parallel or perpendicular to the external stretching, was maintained throughout 2 hr of external stretching. This suggests that average cell elongation alone does not determine colony elongation direction when subjected to uniaxial cyclic stretching.

Contributions to symmetry breaking
We next sought to evaluate quantitatively the contribution of cellular processes to elongation. We quantified the contributions of each cellular event using image segmentation, cell tracking, and neighbor triangulation (Etournay et al., 2015;Etournay et al., 2016) (see Appendix 1A, Figure 6a and b, Figure 6-figure supplement 1 and Figure 6-figure supplement 2, and Video 7). This procedure decomposes overall tissue elongation, which is quantified in terms of cumulative total pure shear, into contributions axis) and cumulative Q xy (right y axis) during 120 min of colony expansion when E-cadherin are blocked by an E-cadherin antibody in control and under cyclic uniaxial stretching (a = 0). Mean value ± standard error of the mean, N control = 3, n = 8 colonies and N stretching = 4, n = 15 colonies. (e) Trajectories of cells (left) and single-cell velocity as a function of its distance to the center of the colony (right) in control colonies (top) and under cyclic uniaxial stretching (bottom). n control = 90 cells from 8 colonies of N control = 4 independent experiments and n stretching = 154 cells from 13 colonies of N stretching = 4 independent experiments. Individual cells in gray, red square and line corresponds mean (binned by distance to the center), shadowed area corresponds to SD. (f) Myosin distribution inside MDCK-GFP-myosin colonies at 0 min and at 120 min after expansion in control and under uniaxial stretching. Note the myosin structures appearing in the stretching case (yellow arrows). Scale bar 10 mm. The online version of this article includes the following source data for figure 3: Source data 1. The Q xx (t)-Q xx (0) and Q xy (t)-Q xy (0) (Figure 6f). Interestingly, we found however that in that case, shear created by T1 transitions is the main contributor for the total pure shear, while cell elongation does not contribute to tissue elongation ( Figure 6f and

Vertex model recapitulates symmetry breaking and shear decomposition
We then asked whether a model could reproduce experimental observations of shear decomposition and, in particular, in which conditions tissue elongation would arise from cell elongation or from topological transitions. We developed a vertex model which takes into account mechanical anisotropies such as active stresses and polarized cell bond tension (see Figure 7a and Appendix 2). We generated a confluent colony of circularly confined cells, in which a unit director p that represented the cell polarity was assigned to each cell. Based on orientation of the director, each cell generated an extensile active stress s a and bias l in the base value of its edge contractility to promote cell elongation and active T1 transitions. We assumed that the experimentally measured cell elongation nematic q is a readout of the underlying cell polarity p ( Figure 7b). Hence, the initial spatial distribution of p was based on the commonly observed pattern of q  direction (quantified by the average of the cosine of two times the angle of each finger trajectory, for example the angle corresponding to the vector between the position of the finger at t = 0 hr and t = 2 hr, of each finger weight by the finger's displacement). N = 5, n colonies = 12 colonies and n finger = 21 fingers. Blue line corresponds to the linear fitting of the data points and the shadowed area corresponds to the 95% confidence interval. Pearson's correlation coefficient r = 0.7724, p=0.0001. (d) Total displacement of finger growth as a function of its initial position in the colony (angular coordinate from the center of the colony) for control colonies and colonies under uniaxial stretching. N = 5 independent experiments, n colonies = 11 colonies and n finger = 21 fingers (control) and N = 6, n colonies = 10 colonies and n finger = 28 fingers (stretching). Averaged fingers in gray (both position and distance, Mean ± SD), blue line corresponds to the linear fitting of the data points and the shadowed area corresponds to the 95% confidence interval. (e) Distance between the tip of the finger and the edge of the monolayer along time, for monolayers in control conditions, stretched parallel and perpendicular to the finger growth direction. N control = 3, N parallel = 4 and N perpendicular = 3 independent experiments and n = 4, 6, and 3 fingers, respectively. Individual fingers in gray, blue line corresponds to the linear fitting of the data points and the shadowed area corresponds to the 95% confidence interval. The online version of this article includes the following source data and figure supplement(s) for figure 4: Source data 1. Raw data corresponding to panels (c), (d) and (e).   Video 5. Disruption of the acto-myosin cable. Movie showing the disruption of the acto-myosin cable. The colony and the micropipette used are shown at the left and the myosin signal is shown at the right. Cytochalasin D was mixed with Cy5 to allow its visualization. Time in hh:mm. Scale bar 50 mm (right).
https://elifesciences.org/articles/57730#video5 Figure 5. Collective effects and symmetry breaking. (a) A leader cell at the boundary of the colony pulls the colony outwards while inner cells deform and elongate. Scale bar 50 mm. (b) Cell shape is quantified using a nematic field (red segments). First, the mean cell shape nematic is quantified at the moment of stencil removal (0 hr) and the orientation q nematic of its mean for the entire colony is obtained. Then, the overall shape of the colony after 2 hr is obtained by fitting an ellipse, whose major axis makes an angle q colony . The yellow directors correspond to fits for the cell shape nematic field obtained with respect to the +1/2 and À1/2 topological defects of the experimentally obtained (red) nematic field (also see Appendix 1D for details.) Scale bar 50 mm. (c) The cumulative distribution function (CDF) for the difference D between q nematic (0 hr) and q colony (2 hr) is obtained from n = 19 colonies of N = 5 independent experiments. Red line corresponds to the CDF of a random distribution of the difference D. This plot shows a strong correlation between the cell shape nematic and the overall shape symmetry breaking (also see Figure 5-figure supplement 2d, and Appendix 1E). (d) The experimentally measured angle of mean nematic orientation q nematic obtained for 19 colonies at t = 0 hr is compared with its counterpart q fit obtained by fitting the experimental data with Equation 1 of the main paper with respect to the orientation parameter a (see Appendix 1D). The size of the red circles in (b) is proportional to the magnitude of anisotropy of the colony shape after 2 hr. n = 19 colonies of N = 5 independent experiments for MDCK. The online version of this article includes the following source data and figure supplement(s) for figure 5: Source data 1. Raw data corresponding to panels (c) and (d).        Figure 7c). To evolve p with time, we imposed that p of the exterior cells tended to be parallel to the boundary, whereas the inner cells tended to align their p with those of their neighbors (Appendix 2D).
Upon removal of confinement, we found that the simulated colony spontaneously elongates, in a direction set by the orientation of the mean cell elongation nematic field, along with the formation of finger-like structures near the +½ defect, as observed experimentally (see Video 8). Our simulations therefore reproduce experimental observations indicating that colony deformation can be understood without forces exerted by a leader cell at the colony boundary ( Figure 7c, Video 8, and Appendix 2). To test whether we could also recapitulate different contributions of the total pure shear inside the tissue, we performed a shear decomposition analysis of in silico movies. We found a qualitatively similar cumulative shear-strain decomposition as was observed in experiments ( Figure 7d and e), where the main contribution came from cell elongation. Moreover, by changing the relative contributions of the cellular active stress magnitude (s a ) and the edge tension bias (l), we could modulate the relative contributions from cell elongations and T1 transitions to the total pure shear (Figure 7-figure supplement 1) as was also observed in experiments with colonies in the absence or presence of cyclic stretching ( Figure 6-figure supplement 3). When s a was dominant, the colony elongation was primarily due to cell elongation, whereas when l was the stronger term, T1 transitions were the main cause of colony elongation. These results reveal possible cellular mechanisms that can govern the process of tissue deformation and influence whether cell elongation, or cellular rearrangements, dominate tissue elongation.
Our vertex model assumed that the cell elongation was the main readout for cell polarity, and it did not explicitly account for the effect of substrate stretching. To incorporate uniaxial cyclic stretching, we further developed the model. Our results show that initial cell shape elongation does not have a preferential direction ( Figure 5-figure supplement 4), but colony elongation under uniaxial cyclic stretching is along x axis, the direction of stretching, and mainly achieved through T1 transitions ( Figure 2c and Figure 6f). Also, we report that the elongation happens along y axis, perpendicular to the direction of stretching, in single cells and cell colonies with lowered Ecadherin levels ( Figure 3a-d). These two experimental observations can be implemented in the model. First, we introduced an additional term m stretch that oriented cell polarization p along x axis, for any given initial condition, upon confinement removal at t 0 -this term is inactive in the absence of uniaxial stretching (Appendix 2D). Then, by using cell active stress s a > 0, we mimicked the tendency of single cells to elongate perpendicular to the orientation of polarity, i.e., perpendicular to the uniaxial stretching. In contrast, the bias in the edge tension l > 0 induced T1 transitions along the polarity of the cell, i.e., parallel to the uniaxial stretching. Thus, in the presence of uniaxial stretching, m stretch oriented cell polarities along x, while the relative magnitudes of single cell active stress s a and edge-contractility bias l dictated the orientation of the colony elongation. When keeping a lower magnitude of s a and a higher value of l, colonies elongated along x through T1 transitions (Figure 7f and Video 9), mimicking colony elongation under uniaxial cyclic stretching (Figure 2  Comparison between the mean total pure shear obtained from TM and the overall colony pure shear obtained from ellipse fitting (left). Total shear corresponds to n colonies = 4 colonies from N = 2 independent experiments and Q xx control was obtained from n colonies = 25 colonies of N = 11 independent experiments. Relative contribution of the different processes to the total pure shear (right). Total shear and contributions were obtained from n colonies = 4 colonies from N = 2 independent experiments. (e) Cumulative pure shear decomposition for stretched colony. (f) Comparison between the mean total pure shear obtained from TM and the overall colony pure shear obtained from ellipse fitting (left) and relative contribution of the different processes to the total pure shear (stretching case). Total shear corresponds to n colonies = 4 colonies from N = 4 independent experiments and Q xx stretching was obtained from n colonies = 20 colonies of N = 9 independent experiments. Relative contribution of the different processes to the total pure shear (right). Total shear and contributions were obtained from n colonies = 4 colonies from N = 4 independent experiments. The online version of this article includes the following source data and figure supplement(s) for figure 6:   Figure 6f). On the contrary, by increasing s a and lowering l (Ecadherin deficient colonies), colonies elongated perpendicular to x (Figure 7g and Video 9), mimicking colonies under uniaxial cyclic stretching treated with E-cadherin antibody (Figure 3c-d). Therefore, we propose that a competition between the strength of active T1 transitions parallel to the external stretching and active cell stress perpendicular to the external stretching dictate overall colony elongation under uniaxial cyclic stretching. When cell-cell junctions are intact, colony elongation is along the direction of stretching and through T1 transitions (Figure 2c and Figure 6f), suggesting that the tendency of single cells to orient along y (Figure 3a-b) is partially screened by cell-cell junctions via T1 transitions. When cellcell junctions are weakened, active cell stress dominates, and colonies, which could be thought of to be closer to a collection of single cells, elongate perpendicular to the uniaxial stretching direction (Figure 3c and d).

Stretching-dependent elongation is mediated by ROCK
We showed that upon stretching, cells reduced their speed and myosin structures appeared (Figure 3e and f). These type of cellular responses to external stretching involve the Rho-associated protein kinase (ROCK) (Hart, 2020), which is also involved in cell-cell contacts integrity (Ewald et al., 2012;Nishimura and Takeichi, 2008). We treated MDCK colonies with a ROCK inhibitor (Y-27632 50mM) and followed their behavior for 2 hr, both in the absence and the presence of cyclic uniaxial stretching (Figure 8-figure supplement 1). When cyclic stretching was applied, elongation along the direction of the uniaxial cyclic stretching was absent (Q xx » Q xy ) (Figure 8-figure supplement  1a and b). However, colonies still elongated anisotropically, similar to colonies in the absence of application of uniaxial cyclic stretching (Figure 8-figure supplement 1c). According to our model, when edge tension bias is sufficiently large, the dominant mechanism for tissue elongation switches from single-cell elongations to T1 transitions (Figure 7-figure supplement 1). We observed that uniaxial cyclic stretching triggers this type of elongation ( Figure 6f). Therefore, if the effect that application of uniaxial cyclic stretching has at the cellular level was reduced, colonies subjected to cyclic uniaxial stretching should preferentially elongate as if cellular active stress becomes dominant, that is through single-cell elongation. Strikingly, shear decomposition analysis shows that the elongation mechanism shifts from T1 transitiondominant to cell elongation-dominant when colonies under uniaxial cyclic stretch have ROCK inhibited (Figure 8a-b). In summary, colonies under cyclic uniaxial stretching elongate through T1 transitions rather than through cell elongation (red dot in Figure 8c) similar to what our model predicts for colonies with increased biased edge tension. In contrast, colonies elongate spontaneously largely through single-cell elongation, which the model predicts when the effect of the cellular active stress is more dominant (green dot in Figure 8c). Strikingly, the application of a ROCK inhibitor leads to single-cell elongation dominating over T1 (orange dot in Figure 8c), effectively suppressing the effect of uniaxial cyclic stretching on the mode of colony deformation (Figure 8-figure supplement 1).

Discussion
Tissue spreading is key during embryonic development (Bénazéraf et al., 2010;Campinho et al., 2013;Etournay et al., 2015). Epithelial cells migrate in a collective and cohesive manner. In many cases symmetry is broken leading to shape transformation. The resulting tissue kinematics and the underlying mechanisms for this symmetry breaking has been studied using original approaches, and is understood with a variety of cell-based and continuum models in vivo (Blanchard et al., 2009;Ciarletta et al., 2009;Aigouy et al., 2010;Etournay et al., 2015;Guirao et al., 2015;Alt et al., 2017). On the in vitro situation too, combination of theory and experiments have been applied on similar phenomena (Mark et al., 2010;Tarle et al., 2015;Kawaguchi et al., 2017;Saw et al., 2017).
In the present work, we sought to integrate knowledge from in vitro and in vivo studies to test new ideas for breaking symmetry through collective effects. First we showed that initially circular colonies of three different epithelial cell lines spontaneously expanded in a non-isotropic manner and the elongation observed was similar in magnitude to Drosophila wing blade elongation in vivo (Etournay et al., 2015). However, the undeformed circular colonies have a non-zero average cell elongation even before spreading starts, which determines orientation of the final colony shape. Our analysis also showed that the cell orientation patterns are not homogeneous but spatially organized and directed by the presence of ±½ topological defects. It was already shown that topological defects in the cell elongation nematic field have a key role on epithelial dynamics (Kawaguchi et al., 2017) and on cell death and extrusion (Saw et al., 2017). Our results reinforce the idea that cell elongation nematic field, which could only arise from collective interaction between cells, can have an impact on epithelial morphogenesis.
Since developing tissues are regularly subjected to internal oscillations (Solon et al., 2009;He et al., 2010;Rauzi et al., 2010) and external pulsatory forces (Zhang et al., 2011;Hayes and Solon, 2017) in a number of living organisms, we explored the effect of external forces in the in vitro circular colonies. We observed that the direction of elongation could be rectified by imposing an external uniaxial cyclic force. This particular behavior is of great interest from an in vivo point of view: cyclic contraction could direct elongation in specific directions to trigger the next steps in morphogenesis. The generic localizations of muscles connected to epithelial layer (Zhang and Labouesse, 2012) could have this essential function for tissues which would otherwise elongate in any direction like in our assay. force arising from edge or cortical contractility L. In our model, each cell also has a polarity p associated with it through which active forces can act on the cell vertices due to anisotropic cell active stress (extensile in our case), cell motility v 0 and polarized or biased edge tension that depends on the orientation of the edge e ab with respect to the polarities of the adjoining cells. When the edge connecting two cells becomes smaller than a critical value cls , the cells are made to modify their neighbors by forming a new edge of length opn > cls as shown. Scheme depicting the different possibilities from the model parameters. (b) Experimental coarse-grained cell shape nematic at t=0 and t=120 min. Scale bar 50 mm.  Video 8. In silico recreation of a colony elongation. Vertex model simulation of a colony elongation. Cellshape nematics, topological defects, and cell polarization are followed over time. https://elifesciences.org/articles/57730#video8 Finally, direction of this tissue rectification is along the external force, but perpendicular to the reorientation of single cell under external uniaxial stretching, and this further supports the collective nature of the phenomenon. We also systematically quantified the shear deformation kinematics of the colony and demonstrated that some of the colonies exhibit shear decomposition patterns similar to those observed during Drosophila pupal wing elongation (Etournay et al., 2015). Moreover, while the colony deformation for the control colonies was dominated by cell elongation, T1 transitions were the main drivers of the colony shape anisotropy under cyclic uniaxial stretching thus indicating different mechanisms at work. Thus, the current work thus makes a direct comparison and contrast tissue kinematics between in vivo and in vitro cases.
Finally, we developed a vertex model that takes into account mechanical anisotropies. Cell anisotropic active stress arising in the cell core, cortical contractility at the cell-cell junctions and cell motility are three of the important forces involved in morphogenesis of epithelial monolayers. To our knowledge, this is the first attempt to systematically show how each of these activities acts on tissue kinematics. Specifically, we showed that cell anisotropic active stress results mainly in cell elongation, whereas anisotropies in cortical contractility primarily effects cell intercalations or T1 transitions. Including cell motilities appears to enhance tissue shear generated by the other two modes of internal forcing. We perturbed active stress and line edge tension by inhibiting ROCK, which has been reported to be involved in cell-cell contact integrity in vivo (Nishimura and Takeichi, 2008;Ewald et al., 2012), and recently, in cell responses to stretching in vitro (Hart, 2020). This led to experimentally blocking the ability of colonies to respond to the externally applied uniaxial cyclic stretching. By doing so, colonies which primarily elongate through cell intercalations, shifted to a single cell elongation driven mechanism.
From our simulations, we could demonstrate that symmetry breaking and finger formation in colonies and the corresponding tissue kinematics observed in our experiments could be brought about by collective active behavior of the colony cells and does not require special action of leader cells (Theveneau and Linker, 2017). This result echoes experiments in which the emergence of the leader cells and the fingering behavior at the border were suggested to arise due to the internal stress pattern in the tissue (Vishwakarma et al., 2018). On the other hand, there are many excellent models in which the boundary cells are ascribed special motility properties that could replicate the experimental results on border fingering (Mark et al., 2010;Tarle et al., 2015). Thus, although leader cells at the boundary may play a role in the border fingering, our experimental findings and simulations clearly indicate that the cell-level internal activities and cell-cell interactions are sufficient to cause symmetry breaking in the colony shape and its overall kinematics via the collective cell-shape nematic field.

Conclusion
Our results show that cell elongation nematic field can have an impact on epithelia morphogenesis. It was already reported that topological defects in the cell elongation nematic field have a key role on epithelial dynamics (Kawaguchi et al., 2017) and on cell death and extrusion (Saw et al., 2017). Now, we showed that circular epithelial colonies when in confinement build up a mean nematic orientation. This symmetry breaking results from the inner activity of cells, and sets the direction for

Fabrication of PDMS membranes and stencils
Poly(dimethylsiloxane) (PDMS) (Sylgard) was prepared by mixing the pre-polymer and the crosslinker at a 10:1 ratio. To prepare stretchable membranes, uncured PDMS was first centrifuged (3000 rpm for 5 min) to remove air bubbles. Afterwards, the PDMS was spin-coated on a flat polystyrene (PS) surface (500 rpm for 1 min) and cured at 65˚C overnight. PDMS stencils were prepared as described previously (Ostuni et al., 2000). Briefly, SU-8 2050 molds containing circles of 250 mm in diameter were prepared by standard photolithography. Uncured PDMS was then spin-coated on molds to a thickness lower than the height of the microstructures (50 mm) and cured overnight at 65C . Stencils for the finger experiments were prepared by spin-coating uncured PDMS on a flat surface.

Cell seeding on stencils
The PDMS stencils were cut, peeled off the mold, and placed in a 2% Pluronic F-127 (Sigma-Aldrich) in PBS for 1 hr. The stencils were then kept in PBS for 2 hr and dried under the hood flow. PDMS stretchable membranes were cut and then activated using O 2 plasma. The membranes and the stencils were exposed to UV light in the cell culture hood for 10 min. Afterwards, stretchable membranes were incubated with fibronectin 20 mg/ml for 1 hr, rinsed with PBS and dried. PDMS membranes were placed on a PS holder, and the PDMS stencils were deposited on top of the membrane right after. A rectangular PDMS chamber was attached onto the membrane using vacuum grease, and cells were seeded at a density of 20,000 cells/mm 2 (Serra-Picamal et al., 2012) for 4 hr. When cells were attached, the medium was changed and the membrane with the cells was kept in the incubator. Local cell density could vary within each colony. We followed the dynamics of assembly of the colony prior removal of the stencil and we could see that cellular clusters size distribution and respective location within the cavity at plating could contribute to these variations. Once they formed confluent circular colonies, the stencils were removed with tweezers carefully before starting the experiment. Some of the colonies exhibited elongation in the short time window between stencil removal and the start of image acquisition.

Time-lapse microscopy
After stencil removal, the medium was replaced by L-15 (Leibovitz's L-15 medium, Invitrogen) supplemented with 10% FBS. Cells were then observed under a Nikon Ti inverted microscope using either a x10 or a x20 objective for 6 hr at 37˚C. Images were acquired every 5 min.

Stretching experiments
The stretching device was designed in the lab. Briefly, a Thorlabs motor (Thorlabs z812 motor, Thorlabs) was controlling the motion of a PDMS membrane, and everything was mounted on a custommade microscope stage. Circular colonies were plated on PDMS membranes, and after removal of the stencils, samples were placed on the microscope. Cyclic uniaxial stretches were applied and images were taken every 30 min typically shortly to prevent interfering with the time course of the experiments. We probed three times for cycles, 20, 60, 120 s, and three extensions, 5, 10, and 20%. The shape of the cycles was triangular. We checked that the PDMS membranes were elastic at all extension and frequency ranges.

Chemical treatments
To prevent the formation of E-cadherin-mediated adhesion, cells were incubated for 30 min with L-15 medium containing 5 mM EDTA (Sigma-Aldrich) and 10 mg/ml anti-E-cadherin blocking antibody that targeted the extracellular domain of E-cadherin (Gumbiner et al., 1988) (uvomorulin, monoclonal anti-E-cadherin antibody, Sigma); after incubation, the medium was replaced by normal L-15 and the experiment started. The inhibition of ROCK was done by incubating cells with Y-27632 50 mM solution (Sigma-Aldrich) from 30 min before the experiment started until the end of the experiment.

Finger dynamics experiments
For the finger test after growth, we let finger grow for 2 hr, and we subsequently applied the cyclic stretch.

Colony shape change analysis
Shape change analysis was performed using ImageJ (http://rsb.info.nih.gov/ij, NIH). The outline of the colony on phase contrast images was ellipse fitted at every time point. Major axis a, minor axis b, and ellipse orientation were obtained. We computed Q, defined as Q xx (t) = ½ ln(a(t)/b(t))Ácos2Á ((t) À a) and Q xy (t) = ½ln(a(t)/b(t))Ásin2Á((t) À a ) being a = (t final ) to quantify cell colony elongation. In uniaxial stretching experiments, the x axis corresponds to the direction of the external stretch and Q components are defined as Q xx (t) = ½ ln(a(t)/b(t))Ácos2Á((t) À a) and Q xy (t) = ½ln(a(t)/b(t))Ásin2Á((t) À a) being a = 0.

Velocity analysis
The centroid trajectories of cells were tracked using the manual tracking plug-in in ImageJ. Data analysis was performed using a custom-made code in MATLAB (The MathWorks). Cell positions were characterized by a vector r(t), with t denoting time and r position in space (bold letter refers to a vector). Every recorded cell position during the time-lapse experiment was defined as r i = r(t i ), where t i = iDt are the times of recording and Dt denotes the duration of time-lapses. The average velocity of each cell was then defined as v = (1/T)Á P i r i , being r i the module of the vector r i and T the total duration of the trajectory.

Colony segmentation and cell tracking
Movies acquired using an MDCK GFP-E-cadherin strain were first pre-processed with FIJI. The subtract background function was applied to remove noise. Images were then loaded to Tissue Analyzer (TA) v8.5 (Aigouy et al., 2010) for edge detection and cell tracking.

Orientation field of the cells and topological defects
First, the background noise of the time-lapse images of the elongating MDCK colonies was reduced with the subtract background function by using a rolling ball radius of 40 px. The resulting images were then subjected to band-pass filter with upper and lower limits of 40 px and 3 px, respectively. The background noise from this output was reduced by using the subtract background command again with a rolling ball radius of 40 px. The processed images from each experiment were analyzed with the OrientationJ plugin of FIJI to quantify their local spatial orientation that reflects the underlying cell elongation. Within this plugin, we used a local smoothing window of 20 px (approximately of the size of the cells) to obtain the structure tensor at discrete points on a grid of 20 px Â 20 px. The plugin provides the dominant direction f i of the structure tensor at each point x i ; y i ð Þ that represents the local orientation field q i ¼ cosf iê1 þ sin f iê2 for cell elongation. The OrientationJ plugin also provides the coherence C of the structure tensor to quantify the strength of the orientation; C » 0 and C » 1 would approximately correspond to rounded and elongated cells, respectively. The orientation or the director field q thus obtained was further quantified by studying the spatiotemporal evolution of ±1/2 topological defects that were obtained by calculating the winding number over unit-cells of the underlying grid. The local smoothing window of 20 px for obtaining the structure tensor, which is approximately of the size of cells, ensured that the most robust defects were observed. The validity of this procedure was cross-verified with the smoothed cell-shape nematic field and the corresponding ±½ defects from the segmented and triangulated data of the experiments processed in Tissue Miner (TM) (Figure 5-figure supplement 2 and Appendix 1). Finally, the orientation of mean cell-shape nematic calculated at 0 hr was compared with shape orientation of the colony at t = 2 hr. For obtaining the cell-shape nematic field for MCF 10A colonies, the procedure was the same as for MDCK control but the numerical parameters used were, rolling ball radius for subtract background 50 px, no band pass filter, and local smoothing window of 15 px and grid size of 30 px for OrientationJ. Similarly, for obtaining the cell-shape nematic field for stretched colonies, the procedure was the same as for MDCK control but the numerical parameters used were, rolling ball radius for subtract background 50 px, no band pass filter, and local smoothing window of 40 px and grid size of 30 px for OrientationJ. See Appendix 1 for more details (also see Video 6).

Acto-myosin cables
In order to image acto-myosin supracellular cables at the boundary of colonies either cells expressing mCherry-actin / GFP-myosin strain or cell immunostaining were used. To disrupt acto-myosin cables, local injection of 4 mM Cytochalasin D using a micropipette was performed. Glass micropipettes were connected to a microinjection system (CellTram vario, Eppendorf). The position of the pipette tip was controlled in x, y, z by using a micromanipulator. The system was mounted on an epifluorescence inverted microscope to record the process. Cytochalasin D was released locally for about 10 min. To detect discontinuities in the cable, fluorescent images of actin and myosin were treated using Fiji as follows. First, a median filter was applied and the background subtracted. Then, a contrast-limited adaptive histogram equalization (CLAHE) was used and the background was again removed using an exponential function. Images corresponding to actin and myosin were added. Finally, a Laplacian of Gaussian filter was applied to the resulting image. Once the cable structure was revealed, the positions of the defects were identified.

Quantification of cellular deformations, topological transitions, and their contribution to pure shear deformation
After obtaining the geometrical and topological information of the colonies from the series tracked images generated using TA, TM was used to extract, triangulate and store the data with the help of an automated workflow. The database obtained after this stage of analysis was used to quantify various state properties such as cell area, neighbor number, cell elongation and the contribution of different cellular processes to deformation using scripts written both in R and in Python. TM was further used to quantify the contributions of various cellular events such as cell elongation and topological transitions to the colony deformation. More details about this analysis can be found in Appendix 1 (also see Video 7).

Vertex model simulations
A vertex model was developed with an addition of a unit nematic director p to every cell. The orientation of the boundary cell p was maintained parallel to the boundary, whereas p for internal cells were modeled to tend to align with the p of its neighbors. In these simulations, an extensile active stress s a ðpp À 1 2 IÞ with s a <0 acts to increase cell elongation along p. In addition, a bias l was also applied on the basic edge tension with respect to the director p of its neighboring cells. For positive l, this bias reduces (increases) the tension of the edge along (perpendicular to) p. Consequently, the closure (formation) of edges is enhanced in the direction perpendicular (parallel) to p. Hence, the T1 transitions in the region around the cells are oriented to cause shear elongation (contraction) along (perpendicular) to p. The colony is provided with an initial condition for p that mimics the initial configuration of experimentally frequently observed cell shape nematic fields q with two +½ defects that are separated by a distance. The initial polar vector orientation along the nematic axis are chosen at random such that the total polarity is zero, and the dynamics of polarity reorientation is independent on a p ! -p flipping of the polarity axis. To begin with, the cell positions and director orientations were evolved under colony confinement until cell position and p do not change significantly. The confinement is then removed to see how the colony breaks symmetry in its shape. In another set of simulations, a small motility v 0 was also provided to the internal cells (Figure 7-figure supplement 1). Similar to the experiments, the output of these simulations was also processed in TM and analyzed for topological defects and pure shear decomposition. For colonies subjected to uniaxial cyclic stretching, p for any cell was provided with an additional tendency to align along the stretching orientation (x-axis). Moreover, the active stress in this case s a >0 was contractile, that is the cell tended to elongate perpendicular to the orientation of p. More details of the simulations are provided in Appendix 2 (also see Videos 8 and 9).

Statistics
No statistical methods were used to predetermine sample size. Measurements were performed on different colonies (  The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. Additional files

Author contributions
Supplementary files . Transparent reporting form

Data availability
All data generated or analysed during this study are included in the manuscript and supporting files.

A. Analysis of cell elongation
Cell-shape anisotropy is an important marker of stress in the cells in response to either internal activity or external forcing (Paluch and Heisenberg, 2009;Rauzi and Lenne, 2011). To understand and quantify the elongation of cells over time, TM computes an elongation parameter proposed in Aigouy et al., 2010. In their study to understand the role of flow in orienting the axis of polarity in Drosophila wings, Aigouy et al., 2010 use a nematic tensor with a magnitude and orientation to quantify the elongation of each cell. For a given cell with area A c , the elongation nematic tensor is proposed as This tensor E is traceless and symmetric with components given by where RðfÞ is the distance to cell boundary at an angle f (SI of Aigouy et al., 2010). For this nematic tensor, the magnitude is given by and the orientation is given by angle f cosð2fÞ ¼ xx ; sinð2fÞ ¼ xy : (4) Based on the above expressions, the spatial distribution of the magnitude of shape nematic for individual cells within the colony for three time points is shown in Figure 6-figure supplement 2ac. The fine-grained nematic segments for individual cells whose magnitude is proportional to and orientation is f (Figure 6-figure supplement 2d-f). The coarse-grained representation of this tensor ( Figure 6-figure supplement 2g-i) is obtained by spatially averaging the elongation nematic over several cells by using a Gaussian smoothing kernel of size comparable to the size of cells (30 pixels in this case). From the figures it is clear that the elongation of interior cells is strongly correlated with the orientation of the future 'finger' even at t ¼ 0 h. B. Analysis of tissue deformation from components of pure shear strain via triangulation As discussed earlier, during different stages of morphogenesis, epithelial tissues are reported to undergo a series of physical changes. The various internal and external forces that are being applied to the tissue during these processes result in deformation and rearrangement of constituent cells. The shape modification rate of the tissue is quantified in terms of the anisotropic component of shear strain rate. However, at the tissue level in epithelial sheets, shear strain can arise due to a combination of cell elongation, T1 transitions, cell division, cell extrusion, and correlation effects. By studying each both component separately and with respect to the other components, one can get hints regarding the mechano-chemical state of the colony. As is described in Etournay et al., 2016, for any general 2-D velocity vector v, the gradient tensor can be defined as The above velocity gradient tensor for a group of cells can be split into three parts: (i) isotropic component to obtain the rate of area change, (ii) symmetric traceless component to quantify the pure shear strain rate, and (iii) an antisymmetric part that represents the local vorticity (Etournay et al., 2015;Etournay et al., 2016)  Here, more specifically, v kk ¼ v xx þ v yy is the trace of the velocity gradient tensor, or the divergence of the velocity field and quantifies the rate of contraction or expansion in the cell area. The traceless symmetric partṽ ij is the pure shear strain rate and quantifies shape changes in the tissue. The antisymmetric component ! is the vorticity of the velocity field and provides values of rotational rates of the group of cells during deformation.
The total pure shear strain rate of the tissueṽ ij , now represented in a tensorial formṼ, can be further split into traceless, symmetric, components as follows where Q is the shear contribution coming from average cell elongation, T is the shear contribution from T1 transitions, C from cell divisions, E from cell extrusions and D due to the correlations effects in the tissue which arise due to coarse graining effects. Detailed derivation and discussion of this expression are available in Etournay et al., 2015;Etournay et al., 2016;Merkel et al., 2017. The contributions from above mentioned components are quantified using the triangulation approach . As shown in Figure 6-figure supplement 2j-l, the triangulation method considers the entire colony without the boundary cells as a network of triangles, connecting the centers of the neighboring cells. Any topological changes in the form of cell division or extrusion will lead to a new triangulation. The shear strain rates in Equation 7 can be integrated over time to provide cumulative shape deformation and the corresponding cumulative shear decomposition for the cell colony, and is obtained for the complete colony without including the external boundary cells (Figure 6b). Figure 4 of Etournay et al., 2015 schematically represents this decomposition of total tissue shear into individual cellular events.
In Drosophila pupal wing the x axis is naturally provided by the A-P axis. In our in vitro situation, however, the x axis is chosen as the dominant direction of colony elongation, such that at the end of the experiment (t ¼ 2 h) the cumulative value of the total shear strain rateṽ xy ¼ 0. C. Correlation between the orientation of mean cell-shape nematic and the orientation of the overall colony elongation As discussed earlier, the shape of a cell can be thought of to be represented by a nematic tensor Q cell . Upon averaging this quantity over multiple cells, we obtain nematic field Q nem ðx; y; tÞ over the colony. As described in the previous sub-section, this nematic tensor is traceless and symmetric with an orientation fðx; y; tÞ (Figure 5-figure supplement 2 and Figure 6-figure supplement 2g-i) and magnitude Q nem ðx; y; tÞ such that Q nem xx ðx; y; tÞ ¼ Q nem cos 2f and (8) Q nem xy ðx; y; tÞ ¼ Q nem sin 2f: This nematic field is obtained for our MDCK colonies by analyzing time-lapse images with the Ori-entationJ plugin of ImageJ (Püspö ki et al., 2016;Saw et al., 2017). We verified that this image based approach gives results that are comparable with those obtained by spatially smoothing (with Gaussian kernel) the individual cell-shape nematics extracted from the segmented images of the colonies (Appendix 1A, Figure 5-figure supplement 2 and Figure 6-figure supplement 2). At any particular time t, the mean nematic field for the colony is Q ncol xx ðtÞ ¼ hQ nem xx ðx; y; tÞi x;y and (10) Q ncol xy ðtÞ ¼ hQ nem xy ðx; y; tÞi x;y : The orientation of the mean cell-shape nematic is then simply where tan À1 is implemented in Matlab using the inbuilt function atan2 (MATLAB, 2018). In this analysis, the coherence of the image structure tensor is the counterpart of the magnitude of the mean cell-shape nematic (Püspö ki et al., 2016). We then compare the final colony orientation colony ðt ¼ 2 hrsÞ, which is obtained from ellipse fitting, with the orientation of the mean colony nematic nematic ðt ¼ 0 hrsÞ (Figure 5b and Figure 5-figure supplement 2). Also refer to Materials and Methods: Orientation field of the cells and topological defects in the Main Paper for additional details.

D. Topological defects as readouts of cell-shape nematic orientation field
Any nematic field can contain topological defects (Kleman and Laverntovich, 2003). A visual inspection of the cell-shape orientation field q ¼ cos fê x þ sin fê y for MDCK colonies indicates the possible presence of topological defects ( Figure 5 and Figure 5-figure supplement 2). The strength k of topological charge for this field (Appendix 1C) in a given region is obtained as where the line integral is calculated by traversing the curve that encloses this region in an anticlockwise sense. For the experimental colonies, the nematic field is not continuous but obtained on the vertices of the grid-cell as shown in Figure 5-figure supplement 5. In this case, the topological charge k enclosed within the grid-cell is where the notation is as shown in Figure 5-figure supplement 5. If no topological defect is present within the cell then k ¼ 0. Else, the net topological charge would evaluate to k ¼ AE1=2, corresponding to the effective strength of the topological defects within the cell.
The location and strength of topological defects thus obtained for the experimental data could then be used to approximately depict the entire nematic field as follows. We first model the energetics of the nematic orientation field with one constant Frank free energy (de Gennes and Prost, 1995). If we then assume that the orientation field f is in local mechanical equilibrium at every timestep, it will then satisfy the Laplace equation r 2 f ¼ 0 (Kleman and Laverntovich, 2003). The resulting orientation field fðx; yÞ for the nematic field in the presence of a single topological defect of strength k at the origin has a singular solution of the form where is the polar angle and a is a harmonic function that satisfies that Laplace equation (Kleman and Laverntovich, 2003). However, for the current paper we find that taking a to be a constant that determines the orientation of the topological defect is sufficient to approximate the orientation of nematic field (Figure 5d). Hence, as per Equations 13 and 15 the integral for a closed curved around the defect gives A ¼ k. If multiple defects are present in the domain, this solution for f can be generalised to give where fx i ; y i g and k i are the coordinates and the strength, respectively, of the i th nematic topological defect obtained from the experimental data, and a dictates the orientation of the defect. Here, the subscript a in f a ðx; yÞ is added to distinguish it from the experimentally obtained fðx; yÞ. As before, tan À1 is implemented in Matlab using the inbuilt function atan2 (MATLAB, 2018). We now use the fitting parameter a in f a ðx; yÞ to see how well it represents the orientation field fðx; yÞ by using the following procedure. The experimental field is obtained on a grid with N points fx k ; y k g, where k ¼ 1; . . . ; N. The difference between f and f a is quantified by the residue obtained over the grid points k. The minimization of the value of r a over a would provide us with a field where f » f a in a mean sense. Interestingly, by using this simple approach with just one fitting parameter a, we can get an excellent representation of the experimentally obtain nematic orientation fieldq (Figure 5d).
E. Statistical test for cumulative distribution of the orientation differences between nematic ð0 hrsÞ and colony ð2 hrsÞ To check if the orientation of the colony elongation at t ¼ 2 >hrs is correlated with the orientation of the mean cell-shape nematic at t ¼ 0 >hrs, we examine the difference D ¼ j colony ð2 >hrsÞ À nematic ð0 >hrsÞj. Specifically, as shown in Figure 5c, we obtain the cumulative distribution function (CDF) for D. If the quantity D is completely random, i.e., colony ð2 >hrsÞ is completely uncorrelated with nematic ð0 >hrsÞ, then the probability density for D 2 ½0 ; 90 is equal to 1=90 . In such a case the CDF(D) as a function of D should simply be a straight line with a slope of 1=90 . However, it can clearly be seen that the experimentally observed CDF is significantly above this line, especially at relatively smaller values of D 30 . This indicates that there indeed is correlation between colony ð2 >hrsÞ and nematic ð0 >hrsÞ.
We now note that the total number of experiments (colonies) that were analysed for Figure 5c is N ¼ 19. To check if the correlation between colony ð2 >hrsÞ and nematic ð0 >hrsÞ holds up for this finite sample size, we do the following. Let us take the probability density pðDÞ ¼ 1=90 ; >8D 2 ½0 ; 90 , i.e., this relative angle D is uniformly random between 0 and 90 . If we now sample N ¼ 19 values for D then we have created a random counterpart to the actual experimental result. We repeat this trial using a straightforward Matlab script and obtain M ¼ 10 6 samples, each with N ¼ 19 experiments. Of all the samples thus randomly generated, for a given D we select those samples i for which CDF random i ðDÞ ! CDF experiment . If for a particular D the total number of such selected samples is mðDÞ, the probability of observing equal or higher bias towards alignment of the two angles from a random experiment than that actually observed experimentally is PðDÞ ¼ mðDÞ=M ( Figure 5-figure supplement 2d). It can clearly be seen that the percentage probability of observing such high values of CDF for smaller D<30 from a completely random sampling of D is » 1, i.e., approximately one sample of N ¼ 19 experiments for every 100 samples, a very low number. The few colonies for which nematic ð0 hrÞ does not match well with colony ð2 hrÞ appear to align with the orientation colony ð0 hrÞ of tissue initial elongation -the associated mechanism is, however, still open (see Figure 5-figure supplement 2e). This entire analysis clearly indicates that the orientation of the mean cell-shape nematic at the initial time has correlation with the final anisotropy orientation of the colony with a high probability.
If we can find a relation between du ab and fdx i a g, then Equation 22 will provide us with expression for F i a . We first note that, assuming an affine virtual displacement field, the virtual displacement dx i a of any vertex i of a given cell k will be given as where the superscript c denotes the centroid of vertices, with position : It can very easily be seen from Equation 37 that for a given cell, the sum of all the active forces on the vertices of that cell will be zero. Moreover, the total torque of all these forces with respect to the vertex centroid as defined in Equation 24 can also be shown to be equal to zero. This is as expected from a symmetric stress tensor. In the current work, we take the cell active stress as extensile, and hence s a <0 (Saw et al., 2017).
Finally, the contribution of this active force to any vertex is a vector sum of the force contributions coming from all the cells that the vertex is a part of, and is given by, ðF i a Þ total ¼ X cell k2vertex i ðF i a Þ total ; a 2 f1; 2g: The active force on the vertex would be added to the forces from cell area deformation and cell edge contractility as discussed in the earlier section.

B. Polarised or biased edge tensions
It is observed that T1 transitions can be influenced by the polarity of the cells (Sato et al., 2015). Such bias in T1 transitions can be introduced in the vertex model by modifying the tension l ab of any edge ab based on its relative orientationê ab with respect to the polarityp a andp b of the shared cells (Figure 7a) as L ab ¼ L 0 ab þ lð1 À ½ðê ab Áp a Þ 2 þ ðê ab Áp b Þ 2 Þ þ ðtÞ; where the quantities L 0 ab , l, ðtÞ are, respectively, the base edge tension, angle dependent edge tension bias, and correlated noise that is calculated as Here, t and DL are the correlation time and strength, respectively, for the correlated noise. This noise promotes T1 transitions and introduces tissue fluidity in the vertex model (Curran et al., 2016). The angle dependent bias ensure that the edges that are parallel (perpendicular) to the polarisation of the cells have a tension that is lower (higher) by l than its base value L 0 ab . This bias enhances the propensity of tissue shear along cell polarisation due to T1 transitions.

C. Inclusion of cell motility
In addition to the forces described in the previous section, cells can also generate motility force (Mogilner and Oster, 1996;Pollard and Borisy, 2003). Following the idea of self-propelled cells (Bi et al., 2016), any cell a can also be modeled to have a speed v 0 , such that in the absence of any other forces it moves with a velocity v 0pa , wherep a is the polarisation direction of the cell as described earlier (Figure 7a). A simplest representation of motile force on a particular vertex i is where N i is the total number of cells that contain the particular vertex i, and b is the index of each of these cells (Sussman, 2017). The motile force, if included, will be in addition to the forces as described in the previous sections.

D. Dynamical evolution
Dynamical evolution in these simulations involves updating the position of vertex i as per the following equation separated by a distance of 1:6R colony with respect to the center of the circle (also see Appendix 1D). An orientation shift of p is randomly added to f p of individual cells such that the vector sum ofp over all the cells is zero. This step is irrelevant for the implementation of forces due to cell active stress and biased edge tensions due to their inherent nematic symmetry (Equations 21 and 39) but ensures that the total motile force from all the cells is zero. The cells are then provided with small active stress (s a <0) and evolved under confinement for some time using the above rules till frg and fpg are almost equilibrated. This configuration is then used as the initial condition for all the numerical experiments, in which the confinement is removed to allow the colonies to elongate (see Video 8). Different simulation conditions correspond to varying proportion of active cell stress ðs a <0Þ, biased cell edge tension ðl>0Þ and cell motility ðv 0 Þ. The results from the simulations are then analysed in a similar manner as for the experimental data (Figure 7-figure supplement 1).
This set of polarity evolution rules for the interior and the boundary cells ensure that in the steady-state the orientation field f p satisfies in a coarse-grained sense, i.e., on longer length-scales, the Laplace equation r 2 f p ¼ 0. The nematic field for the interior cells tends to have a uniform orientation, whereas the nematic field for the boundary cells tries to align along the colony boundary. Because of this boundary alignment tendency of the border cells, at least for the parameters used in our simulations (Appendix 2E), the total topological charge of the polarisation field is mostly conserved in our simulations, i.e., remains equal to the total initial charge of þ1 (for a relevant study in another context see van Bijnen et al., 2012.) Due to the action of active cell stress (Equation 21) the cells have a tendency to elongate predominantly alongp thus creating a cell-shape nematic orientation field q ¼ cos fê 1 þ sin fê 2 (also Appendix 1A-C.) Hence, the cell elongation anisotropy, as described earlier in Appendix 1, can loosely be thought of to be a readout of the internal polarisation of the cells. However, since the actual cell shape results from the interplay between cell mechanical properties and active cell stress,p and q are not necessarily always fully aligned. Thus, for the cell-shape nematic field, the total topological charge can be less than þ1 if the aligning tendency of internal cells dominates, e.g., when the cell active stress is large (see Video 8).
When the effects of external uniaxial stretching are included in the model, we take the m stretch term to dominate due to which the polarisation of the cells align along the orientation of uniaxial stretching AEê x . Hence, in this case, the role of initial and boundary conditions, as discussed above, is not quite important. As described in the section "Vertex model recapitulates symmetry breaking and shear decomposition" of the main paper, the competition between the cell active stress (s a >0) and biased cell edge tension l>0 determines if the colony elongates either along the orientation of stretching or perpendicular to it (see Video 9).

E. Parameters used and the non-dimensional groups present in the Vertex Model
Using reference scales for length (l 0 ), time (t 0 ) and energy (U 0 ), we get the following set of nondimensional groups for the various parameters that are involved in the simulations. The non-dimensional counterparts of the parameters used in the simulations have been assigned an additional 0 to make them distinct from the original.