Rigidity percolation uncovers a structural basis for embryonic tissue phase transitions

1Embryonic tissues exhibit structural properties of viscoelastic materials, the precise spatiotemporal regulation of which is essential for proper embryo development. Recent direct measurements of tissue material properties, such as viscosity, yield stress, and Young’s modulus, combined with genetic and mechanical perturbations have illustrated that spatial and/or temporal changes in material characteristics can affect morphogenetic processes, such as tissue spreading and body axis elongation. For instance, the dynamics of tissue spreading have been shown to be controlled by changes in tissue material properties without altering force application. This suggests that regulating tissue material properties constitutes an alternative mechanism of controlling tissue morphogenesis, in addition to other more direct mechanisms, such as the regulation of force generation. Yet, how tissue material properties change within the developing organism remains poorly understood.

2Intriguingly, changes in tissue material properties can be fast and drastic, resembling phase transitions (PTs), a phenomenon playing a fundamental role in the dynamics of many complex systems. In physics, PTs refer to a broad class of phenomena mainly characterized by the presence of abrupt changes in some of the macroscopic properties of the system, known as “order parameters,” as a consequence of smooth variations of a “control parameter,” when the latter reaches a critical value. At such “critical points,” universal core physical features can be observed, such as singularities in macroscopic observables or power-law distributions in the statistics of order parameters. These patterns are largely independent of the specific features of the system, providing a unifying, simple mechanism determining collective behaviors of many systems of disparate nature.

3Different classes of theoretical models have been proposed in the past decade to explore the possibility of tissues undergoing PTs. On the one hand, active particle models have been developed to explain density-dependent glassy dynamics observed during the collective migration of in vitro monolayers. On the other hand, vertex models have been proposed to describe the geometrical and material properties of confluent epithelial monolayers. In these models, specific geometric features of cells, such as their shape arising from cell-cell adhesion and cytoskeletal forces (control parameters), trigger, when reaching a critical value, large changes in tissue rigidity (order parameter) in a density-independent manner. Notably, in an unjammed state, cells can move at zero energy cost, thereby connecting the amplitude of relative cell movements to density-independent rigidity transitions. However, although such motility transitions can be observed in in vitro cell monolayers and some in vivo settings, direct rheological measurements demonstrating actual sharp changes in tissue material properties during such (un)jamming transitions remain challenging. This has started to be addressed in 3D embryonic tissues by methods such as micropipette aspirations (MPAs) or droplets, although the lack of dynamic measurements of a tissue undergoing a putative PT both in space and time impedes a proper characterization of the process. Moreover, from a statistical physics perspective, there are core universal features that are expected at critical points of tissue PTs, but they have not yet been traced in in vivo systems, in particular due to these experimental limitations. Finally, from a developmental biology perspective, the mechanisms underlying the robust physiological regulation of critical points in vivo remain largely unknown.

4We have tackled these outstanding questions by analyzing fluidization of the zebrafish blastoderm, a non-confluent tissue making up the blastula. By combining rigidity percolation theory with micropipette rheological measurements, quantitative image analysis, biophysical modeling of cell-cell contacts, and genetic perturbations, we show that the blastoderm undergoes a genuine rigidity PT, and that the cell-cell adhesion-dependent topology of the cell-cell contact network is an accurate predictor of its tissue material properties. Furthermore, we predict and observe key signatures of criticality in vivo and identify cell-cycle synchrony as a regulatory mechanism conferring robustness to the blastoderm when undergoing a PT.

5In order to address whether and how embryonic tissues undergo PTs, we turned to the zebrafish blastoderm (Figure 1A). We have previously found, using MPA, that the viscosity of the central part of the blastoderm abruptly drops at the onset of morphogenesis by more than an order of magnitude within a few minutes (Figures 1A–1C, yellow box; Figure S1A; Video S1; see STAR Methods). This tissue fluidization was only transient, with the blastoderm returning to its initial viscosity values before fluidization by gradually increasing its viscosity in a slow (∼1 h) “thickening” (defined in a rheological sense as increasing viscosity) process (Figure 1C, purple box; Video S1). We first asked which microscopic cellular process may trigger such stark macroscopic blastoderm viscosity changes. Given previous reports implicating changes of cell motion and shape in density-independent unjamming transitions of confluent epithelial tissues, we first asked whether blastoderm fluidization might be achieved by those processes. To this end, we analyzed the mean squared relative displacement (MSRD) of blastoderm cells and cell shape index (see STAR Methods) during blastoderm tissue fluidization and thickening, but found no clearly recognizable temporal correlation between these cellular parameters and changes in blastoderm viscosity (Figures S1B, S1B’, S1C, and S1C’). This argues against the blastoderm undergoing a motion-dependent/density-independent unjamming transition. Moreover, the increase in cell number in time via proliferation, which was proposed to give rise to density-dependent jamming transitions in vitro, could not explain the viscosity changes, as analysis of blastoderm cell nuclei density did not reveal any obvious correlation with blastoderm viscosity (Figures S1D and S1D’).

6However, we observed that, upon labeling the interstitial space within the blastoderm, cells were separated by interstitial fluid accumulations (Figures 1D, S2A, and S2B), suggesting that the blastoderm tissue is non-confluent. Importantly, when expressing cell density not as nuclear density, but as the percentage of space not occupied by interstitial fluid (cell fraction), we detected a clear nonlinear relationship between cell fraction and tissue viscosity, independent of variations in cell size (Figures S1E, S1E’, S1F, and S1F’). Such a nonlinear relationship between cell fraction and tissue viscosity is reminiscent of the physics of granular sphere packings, where increasing the packing fraction above a critical point triggers a jamming transition, characterized by the appearance of a non-zero shear modulus despite structural disorder. This occurs when contacting spheres are close enough to form a “rigid” graph of contacts, so that displacing any such spheres requires energy expenditure (see schematic in Figure S1G). Thus, we sought to explore this analogy in a biological setting by creating a temporal map of cellular connectivity, with cell-cell contacts being defined by the absence of interstitial fluid between neighboring cells (Figures 1D, 1E, and S2C). We found that blastoderm cell connectivity was directly related to cell fraction (Figures S1H and S1H’) and consequently displayed a nonlinear relationship with blastoderm viscosity (Figure 1F). Specifically, before blastoderm fluidization, when blastoderm tissue viscosity was largely unchanged, cell connectivity was smoothly and gradually reduced (Figures 1E, 1F, and 1F’, blue-green time points). This reduction in cell connectivity continued up to a threshold value of ∼3–4 contacts per cell, when blastoderm viscosity abruptly dropped at the onset of morphogenesis (Figures 1E, 1F, and 1F’, yellow time points), followed by a slow increase in cell connectivity during blastoderm thickening (Figures 1E, 1F, and 1F’, purple time points). Intriguingly, this transition point of connectivity is very close to the one predicted by rigidity percolation theory (two-thirds of maximal average number of contacts) separating floppy versus rigid connectivity graphs, suggesting that the blastoderm might undergo a rigidity percolation transition.

7To test whether the blastoderm indeed undergoes a rigidity percolation transition, we first generated in silico 2D triangular graphs of different sizes and connectivity levels typically used to characterize the network of contacts in a standard sphere packing (Figure S1I) and performed a rigidity analysis to determine which tissue regions are expected to be rigid (Figure 1G). A cluster is considered as rigid when all its degrees of freedom are absorbed by the structure and thus cannot be deformed without expense of energy (Figure 1G’, green cluster). As previously described, we found above the critical point of two-thirds of the maximum average connectivity the appearance of a rigid “giant cluster” (GC), which percolated over the size of almost the entire system (Figures 1G, 1G’, and S1G, red cluster). In this framework, the average connectivity of the network represents the control parameter, with two-thirds of the maximum connectivity (i.e., ∼4 contacts per cell), constituting the critical point for a rigidity PT. The size of the GC, in turn, corresponds to the order parameter of the system, being small below the critical point and close to the size of the system above it (Figures 1G and G’). Crucially, the presence of a GC implies that the network will resist any deformation, whereas in the absence of the GC, deformations and rearrangements of network nodes and regions are possible at no energy cost (Figure S1G). Therefore, cell connectivity, and thus the size of the GC, might represent a simple putative criterion for predicting the material property of a tissue, and corresponding PT, from purely topological metrics.

8Notably, classical rigidity percolation considers networks of elastic springs, giving rise to zero elastic modulus below the critical point and linearly increasing shear modulus as a function of distance above it. By contrast, the abrupt tissue material changes that we observed within the blastoderm are changes in viscosity. Such viscous response of tissues can arise from turnover or re-arrangements. To incorporate this behavior in a minimal way within our rigidity percolation network analysis, we analyzed 2D networks of elastic springs where the bonds can spontaneously break and heal (with a characteristic timescale τ; see STAR Methods for details of this modeling strategy), revealing that the critical point of rigidity percolation triggers a transition from low to high viscosity (Figures 1H, H’, and S1J; see STAR Methods). This suggests that rigidity percolation can be used to describe the changes in tissue viscosity observed within the blastoderm.

9To directly test the relationship between cell connectivity, GC size, and tissue viscosity in vivo, we sought to compare simulated random networks with experimental 2D cell connectivity networks within the blastoderm. We reasoned that a 2D analysis is applicable for the 3D tissue context of the blastoderm, given that no systematic changes in 2D network connectivity were found between 2D sections of the blastoderm at different heights within the tissue (Figures S1K and S1K’) and that measuring average connectivity changes in 3D showed a very similar pattern as observed in 2D (Figures S1L and S1L’), with a decrease of connectivity at fluidization up to a value close to that expected for the critical point of a rigidity PT in 3D networks. We then compared simulated random 2D networks of similar average size with those typically imaged within the blastoderm (∼100 cells per network) (Figures 2A and 2B; see STAR Methods). We found that not only the relationship between GC size and network connectivity was conserved between simulated and experimentally observed connectivity networks (Figures 2A and 2B), but also that the temporal change in the experimentally observed size of the GC in such networks (Figures 2A and S2D; Video S2) closely matched the temporal changes in tissue viscosity measured during the blastoderm fluidization/thickening process (Figure 2C): high tissue viscosity (before fluidization and after thickening) corresponded to rigid networks with connectivity above the critical point and large GCs, and low tissue viscosity (during fluidization) corresponded to floppy networks with connectivity below the critical point and small GCs (Figures 2C and 2D). Taken together, these findings show that the changes in blastoderm viscosity can be explained by changes in cell connectivity, which control overall tissue rigidity, as predicted by rigidity percolation theory.

10We next sought to experimentally perturb blastoderm cell connectivity to challenge the percolation model and address how reliably cell connectivity can predict tissue rheology (Table S1). To modulate cell connectivity, we interfered with cell-cell adhesion directly by reducing the levels of E-cadherin expression via injection of e-cadherin-morpholino (MO) (Figures S3B and S3C) and indirectly by analyzing different regions of the blastoderm or changing cell fate specification within the blastoderm. Directly interfering with cell connectivity by diminishing E-cadherin expression reduced global cell network connectivity within the blastoderm (Figures 3A and S3G) and, as predicted by the theory, GC size to values below the critical point during the period when control-MO-injected embryos were undergoing tissue fluidization/thickening (Figures 3A, 3A’’, 3D, and S3A–S3C). Importantly, measurement of tissue viscosity in e-cadherin morphant embryos during this period also revealed low values of viscosity, as expected for networks below the critical point within the floppy regime (Figures 3A’, 3E, S3F, and S3H).

11Furthermore, we had previously shown that the blastoderm margin, in contrast to the blastoderm center, does not undergo tissue fluidization and that this is due to the restricted expression of Wnt/PCP pathway components within the margin (Figure S3D). Consistent with our previous observations, we found that cell connectivity remained persistently high in the blastoderm margin during the period when the blastoderm center was undergoing fluidization (Figures 3B and S3G). As theoretically expected, we found that this high level of connectivity was accompanied by the presence of large GCs and high tissue viscosity (Figures 3B–3B’’, 3E, S3F, and S3H). By contrast, marginal tissues of mutant embryos for the Wnt/PCP pathway component wnt11f2 (slb), which undergo marginal tissue fluidization in a similar manner as the blastoderm center (Figures S3E, S3F, and S3H), displayed low connectivity values (Figure S3G) and small GCs (Figures 3E, S3E, and S3F).

12Finally, to indirectly interfere with cell connectivity within the blastoderm by changing cell fate specification, we ectopically induced mesodermal cell fate, which in wild-type (WT) embryos is confined to the blastoderm margin, in the blastoderm center (Figure S3D). Given that mesoderm specification induces the expression of Wnt/PCP components within the blastoderm margin (Figure S3D) and, consequently, high cell connectivity and tissue viscosity within this area, we predicted and found that ectopic mesoderm induction in the blastoderm center not only abrogated tissue fluidization but also raised the level of cell connectivity consistently above the critical point during the period when control embryos underwent tissue fluidization and thickening (Figures 3C–3C’’, 3E and S3F–S3H). Collectively, our perturbation experiments show that all of the tested blastoderm cell connectivity networks collapse on the same master curves describing the relationship between cell connectivity, GC size (Figure 3D), and tissue rheology (Figure 3E), suggesting that rigidity percolation analysis reliably predicted changes in tissue viscosity on the basis of associated changes in cell connectivity within the blastoderm. This central role of cell connectivity in determining tissue viscosity was further confirmed by our observation that interfering with blastoderm cell motility (Figures S3I and S3I’), but not connectivity (Figure S3J), by mildly reducing myosin II activity in those cells, had no recognizable effect on GC size and blastoderm viscosity (Figures 3E and S3F–S3H).

13A defining feature of high-order PTs is the divergence of specific observables near critical points as well as the appearance of power-law distributions in observables related to the order parameter. Our simulations showed a clearly defined power-law pattern in the distribution of cluster sizes other than the GC close to the critical point, with an exponent around 2.5 (Figures 3G and S3K; see STAR Methods for fitting techniques). This critical exponent implies that the variance of rigid cluster sizes (outside the GC) should diverge at the critical point in infinite systems. In the case of finite systems, traces of this singularity remained as a sharp peak in the variance of rigid cluster sizes around the critical point, whose height, given the critical exponent found, should diverge with system size as ∼N1/2, where N is the number of nodes (Figures S3L and S3M; see STAR Methods). The average of rigid cluster sizes, on the contrary, is nearly constant across system sizes and connectivity values (Figure S3N; see STAR Methods). To test the existence of such universal features within the blastoderm, we compared the average and variance of rigid cluster sizes as a function of connectivity between experimental and theoretical networks, which revealed a close match in a parameter-free manner (Figures 3F and S3N). We then went one step further and computed the experimental rigid cluster size distribution outside the GC for 30 networks from all of the above-described experimental conditions lying near the critical point (see STAR Methods). As predicted by the numerical experiments, and despite the relatively small size of the system, the experimental cluster size distribution near criticality closely followed the theory, displaying a power law with exponent around −2.5 with a correction for small cluster sizes (Figure S3K; note that this exponent is well defined only for large systems; see STAR Methods for details) (Figure 3G; STAR Methods). Collectively, these findings reveal key hallmarks of criticality in the blastoderm, further supporting the hypothesis that the blastoderm indeed undergoes a rigidity percolation PT.

14Our results so far demonstrate that average cell-cell contact number is sufficient to predict global contact topology (such as GC size), which itself serves as a reliable predictor for tissue material properties across time and experimental perturbations. But how is contact topology regulated from a biomechanical and biophysical perspective? Our previous findings had shown that cell-cell adhesion decreased due to mitotic rounding at the onset of blastoderm fluidization. To evaluate the role of cell-cell adhesion in determining global contact topology, we developed a mechanical toy model of a 2D cluster of 4 cells organized in a rhombus, one of the simplest non-rigid topologies (Figures 1G’ and 4A), based on previous studies of cell-cell adhesion and embryo compaction. We considered constant cell volume, so that the ratio of cell-cell tension to cell-fluid tension, which defines cell-cell adhesion strength, is the single parameter specifying cluster configuration (Figure 4A). For, cells are non-adhesive as they round up to minimize their total surface energy, while for, cells become adhesive by forming cell-cell contacts that result in surface energy minimization. Consequently, increasing cell-cell adhesion (decreasing ) in our 2D cluster of 4 cells organized in an initially floppy rhombus (Figure 4A, gray cluster, floppy regime) led to cell-cell contact expansion up to the point where tricellular junctions appeared, forming an additional contact in the cluster that collapsed the central interstitial space and transformed the configuration into a rigid topological structure (Figure 4A, green cluster, rigid regime).

15To test whether this could be a mechanism to explain how cell-cell adhesion determines the experimentally observed transitions between floppy and rigid cell connectivity networks within the blastoderm, we determined relative cell-cell tension across different time points during blastoderm fluidization/thickening by measuring the contact angle between blastoderm cells at their interface to the interstitial fluid (Figures 4A and 4B; see STAR Methods for determining ). We found an increase in, and thus a decrease in cell-cell adhesion, during blastoderm fluidization (Figures 4B and B’) and qualitatively similar local cellular connectivity patterns as theoretically predicted (Figure 4A). Moreover, performing the same analysis across multiple time points and experimental conditions (e-cadherin-MO, marginal blastoderm, mesoderm induced) (Figure S4A), revealed a close correlation between and the average connectivity (Figure 4C), GC size (Figure S4B), cell fraction (Figure S4C) as well as tissue viscosity (Figure 4D), suggesting that is functionally linked to those properties. Remarkably, we also found that blastoderm fluidization occurred for values of that are close to the values theoretically predicted by our mechanical model for the transition between floppy and rigid local topologies (Figures 4A and S4A). Likewise, simulating cell-cell contact formation as a function of for large non-confluent 2D cell aggregates using our mechanical model not only produced cell connectivity maps qualitatively matching our experimental cell-cell connectivity graphs but also showed a floppy-to-rigid network transition for similar values as found in the 4-cell mechanical model (Figure S4E). Together, this suggests a simple explanation of how alterations in cell-cell adhesion can trigger global cell network topology and rigidity changes, in agreement with the ones observed during blastoderm fluidization/thickening process.

16Although the blastoderm displays hallmarks of criticality, suggesting that small connectivity changes can already trigger a rigidity PT, the actual changes in blastoderm viscosity are highly uniform, robust, and reproducible (Figure 1C). This implies that the cell network connectivity within the blastoderm must be tightly controlled in development. To obtain insight into these control mechanisms, we turned to the role of the last two rounds of meta-synchronous blastoderm cell cleavage cycles (∼1 h), which we have previously found to trigger a gradual reduction in cell connectivity. Given that cell divisions are randomly and thus homogeneously distributed throughout the blastoderm due to the synchrony of the cleavage cycles (Figure 5A, small variance in the fraction of dividing cells between the quadrants), it is conceivable that synchronous cell divisions are important for uniform tissue fluidization within the blastoderm center. To experimentally address this, we sought to interfere with cell-cycle synchrony in order to create a spatial heterogeneity in the distribution of dividing cells within the blastoderm and, consequently, break the spatial homogeneity of cell contact loss.

17To this end, we overexpressed Chk1, previously shown to slow down the cleavage cycle oscillator, and we found that it also reduced cell-cycle synchrony (Figures S5A–S5D), leading to an increase in the heterogeneity of the spatial distribution of dividing cells within the blastoderm (Figure 5A’). Consequently, Chk1-overexpressing embryos showed a highly heterogeneous spatial connectivity profile within the blastoderm (Figures 5A’ and 5B), altering its structural properties (Figures 5C–5E). The latter was evidenced by a broader rigid cluster size distribution than observed in WT embryos and theoretically expected for random networks of the same size (Figures 5C and S5E; see STAR Methods). However, by incorporating local correlations in links/adhesion in our simulations to account for the experimentally observed heterogeneous spatial connectivity profile (Figure S5F; STAR Methods), we were able to closely recapitulate the cluster size distribution in Chk1-overexpressing embryos (Figure 5C, orange-shaded line versus orange circles). Moreover, consistent with predictions from those simulations, we experimentally observed that the size of the 2nd biggest cluster after the GC (Figures S5G and S5G’; Figure 5D, orange cluster; Videos S2 and S3) was much bigger in Chk1-overexpressing embryos than in WT networks (Figures S5G and S5G’). This size increase partially scaled with the degree of heterogeneity in the spatial distribution of dividing cells (Figure S5H), thereby directly linking the randomness in cell divisions, and thus cell-cycle synchrony, to the structural properties of the blastoderm.

18Finally, we made use of Chk1-overexpressing embryos to assess how a heterogeneous spatial connectivity profile within the blastoderm affects the robustness of blastoderm fluidization. To this end, we measured tissue viscosity in Chk1-overexpressing embryos throughout the period of blastoderm fluidization and thickening in WT embryos. Strikingly, we found that the blastoderm in Chk1-overexpressing embryos did not show one distinct fluidization and thickening cycle, as seen in WT embryos, but instead went through several cycles of partial fluidization and thickening (Figures 5F and S5I), in agreement with the observation that the GC in these embryos appears and disappears several times (Figure 5D’; Video S3). Moreover, the variation in viscosity values obtained at distinct time points during the fluidization/thickening process was considerably higher in Chk1-overexpressing compared to WT embryos (Figure S5J). Strikingly, when comparing the variation in blastoderm viscosity to the heterogeneity in the connectivity of these embryos during the meta-synchronous cell cleavages, we detected a clear anti-correlation between these two features (Figure 5G). Collectively, these findings suggest that cell cleavage synchrony serves as an effective mechanism for uniform cell connectivity changes within the blastoderm that is required for uniform and robust blastoderm fluidization.

19Our findings demonstrate that the zebrafish blastoderm undergoes a genuine rigidity PT and provides mechanistic insight into the structural origin and developmental regulation of this PT within the organism.

20The major challenge in probing the structural basis of a PT—which we have addressed in this study—is to link microscopic control parameters to macroscopic rheological changes, an approach requiring combined measurements of local cell behavior with direct tissue-scale rheological measurements. In previous studies on confluent monolayers, material properties were indirectly inferred from a simple criterion derived from vertex or Voronoi models, the normalized cell perimeter (which coarse-grains details of the cell contact mechanics;;;;; ), as direct rheological measurements were difficult, hampering the characterization of a bona fide PT in those tissues. Recent in vivo studies in zebrafish have managed to probe tissue fluidization during development, specifically within the presomitic mesoderm (PSM) and early blastoderm, via direct rheological measurements, allowing to experimentally probe tissue PTs. However, a theoretical framework for modeling and inferring PTs in such non-confluent tissues was lacking. By applying rigidity percolation, a theoretical concept widely used in material science, to biological tissues, we provide a simple criterion readily testable in any non-confluent tissue, without requiring knowledge of the detailed microscopic properties (which are coarse-grained within the cell-cell contact topology). Importantly, the zebrafish blastoderm provides a system in which a direct rheological transition can be monitored in both space and time, a key prerequisite for testing our framework and pinpointing hallmarks of PTs. We show that rigidity percolation accurately predicts the material properties of the system throughout morphogenesis as well as across multiple experimental perturbations of cell fate, division, contractility, and adhesion. Furthermore, we demonstrate hallmarks of criticality in vivo, such as power-law exponents and its associated discontinuities of macroscopic observables at criticality. Given that those hallmarks have not yet been clearly demonstrated in biological settings—mostly due to difficulties in parameter fine-tuning, measurement precision, and finite size statistics—their identification provides additional strong support for biological tissues undergoing material PTs.

21Our findings suggest that analyzing cell-cell contact topology is sufficient to explain tissue PTs during morphogenesis. But how is this coarse-grained observable mechanistically regulated? For instance, our MPA experiments indicate that the blastoderm behaves largely as a simple viscous fluid, while classical rigidity percolation of spring networks provides information about elasticity. We have addressed this apparent discrepancy by showing that viscous networks, where contacts can break and heal, qualitatively change their response to external deformations at the same critical point of connectivity as elastic ones. This suggests that using rigidity percolation theory is appropriate to describe the viscosity changes observed within the blastoderm, although quantitative understanding of the absolute value of viscosity will require more detailed modeling of the nonlinear and time-dependent properties of the cell contacts. Indeed, for the jamming of sphere packings, spring network models, such as the ones we used, can only predict the critical point and linear response upon shear, while other responses such as bulk modulus or nonlinear rheology would require further extension of the theory to include mechanical elements such as re-arrangements, cell-fluid interaction, and cell-cell friction. The influence of adhesion, for instance, has been probed experimentally using a biomimetic system of attractive emulsion droplets. Adhesion is theoretically predicted to markedly change the nature of the jamming transition, leading to a slow growth of the rigid cluster as a function of sphere fraction around criticality, closely matching our findings (Figure S4D).

22Interestingly, decreasing adhesion in 2D vertex models of fully confluent tissues has been shown to increase tissue rigidity. However, such an effect would only apply up to the point where loss of adhesion breaks tissue confluency and where cell-cell connectivity then becomes key. To gain insight into this process, we have analyzed the microscopic basis of this framework, via a model of cell-cell contacts arising from minimization of cell-cell and cell-fluid surface energy, suggesting that changing the ratio of both tensions (absolute values of adhesion) is enough to transit between rigid and floppy network configurations. Similarly, have extended the vertex model to incorporate fluid pockets, cell-cell adhesion, and active fluctuations in order to describe a fluidization transition within the PSM, as assessed via local micro-rheology and a caged-to-diffusive transition in the MSRD of cells. This fluidization transition in PSM has further been proposed to arise from changes in cell-cell tension/adhesion fluctuations, together with some changes in the absolute values of these tensions influencing cell-cell connectivity. Intriguingly, our functional experiments, where MSRD and cell-cell tension/adhesion fluctuations were reduced by partial inhibition of myosin II activity without affecting cell-cell connectivity changes (Figures S3G and S3J), showed no effects on rheological tissue properties (Figures S3I, I’, S3J, S5K, and S5K’’), consistent with the observation in WT embryos that contact fluctuations do not correlate well with changes in blastoderm tissue viscosity (Figures S1B and S5K). This suggests that the specific cellular properties controlling tissue PTs might be highly context dependent, while the framework of rigidity percolation can provide generic insights into material properties in a parameter-free manner based on a single coarse-grained measurable observable. More generally, applying the concept of rigidity percolation to biological tissues allows importing key concepts and the multi-faceted toolbox of statistical physics to understand the universal features of tissue material properties, although the specific properties of different tissues are likely, as in material science, to depend on different microscopic parameters within this framework.

23Finally, being close to a PT has been proposed to be advantageous, allowing a tissue to change its material properties rapidly and drastically. However, it also might pose a risk as the tissue, when placed close to criticality, could ectopically change its material properties due to noise. To avoid ectopic PT, cell-cell contact loss must therefore be tightly controlled in space and time. Our findings suggest that the meta-synchrony of the cell cycle provides such control mechanisms, where largely synchronous cell divisions lead to random and thus homogeneously distributed contact loss throughout the tissue. This is consistent with previous findings that have linked cell divisions to tissue fluidity. We also show that cellular fate modulates cell-cell connectivity and thus distance to the critical point, hinting at an intricate relationship between cell division, fate specification, and tissue material properties in development, with potential application in pathological processes such as tumor invasion.

24In our pipette aspiration assays (Figure 1B), we have restricted ourselves to the study of the regime where the deformation distance grows linearly in time (consistent with a simple Newtonian fluid), which allows us to extract a unique rheological quantity, i.e., viscosity (Figure 1C). In the future, it would be important to probe different regimes systematically with different pipette aspiration pressures, to understand potentially more complex rheological behaviors, such as yield stress or viscoelasticity. Moreover, while we have shown that the transition in rheological property of the blastoderm is well predicted by rigidity percolation (Figure 1H), a more quantitative understanding on what sets the absolute values of rheological quantities would require an extension of the model, taking into account cellular shape but also the mechanics of cell contacts, including frictional dissipation and re-arrangements. To complement this experimentally, it would be important to further qualify what constitutes a cell-cell contact. Here, we have relied on a binary definition (presence or absence of interstitial fluid between two membranes; see Figure 2A) and have assessed average contact strength at a given time point via junctional angle measurements (Figure 4), but studying the heterogeneities of contact strength both in time and space might help in better understanding the overall tissue material properties.

25Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Carl-Phillip Heisenberg (heisenberg@ist.ac.at).

26This study did not generate new unique reagents.

27The Fiji macro used to quantify the MPA experiments and the codes used for the rigidity simulations, correlation model and the deformation of elastic/viscous networks are available at https://github.com/bernatze/Lattice_Deformation.

28Zebrafish (Danio rerio) were maintained under a 14 h light/10 h dark cycle. The following zebrafish strains were used in this study: wild-type (WT) AB and slb/wnt11 f2. Zebrafish embryos were grown at 25-28.5°C in E3 embryo medium and staged as previously described. For precise staging before and during blastoderm spreading, the last rounds of meta-synchronous cleavages and resulting changes in cell size were used as temporal hallmarks defining developmental time relative to the onset of doming. Chk1-overexpressing (oe) embryos displayed approximately 1 h of developmental delay, thus we extended their analysis by 1 h. Onset of doming (t = 0 min) was defined as the time when during the viscosity measurements the first embryo within an egg-lay/batch underwent fluidization. All animal experiments were carried out according to the guidelines of the Ethics and Animal Welfare Committee (ETK) in Austria.

29Zebrafish embryos were injected using glass capillary needles (30-0020, Harvard Apparatus, MA, USA), which were pulled by a needle puller (P-97, Sutter Instrument) and attached to a microinjector system (PV820, World Precision Instruments). Microinjections of mRNAs and morpholinos (MOs) were performed at the one-cell stage. mRNAs were synthesized using mMACHINE SP6 kit (Ambion). The following mRNAs were injected: 70 pg membrane RFP, 70 pg H2B GFP, 200 pg CAMypt1, and 72 pg chk1. The following MOs were injected: 4 ng e-cadherin MO, 4 ng control MO for e-cadherin. To induce mesodermal progenitor cell fate, the following combination of mRNAs and MOs was injected: 100 pg ndr2l cyclops mRNA and 2 ng casanova MO. To label the interstitial fluid, 1 nL of 0.6 mg/ml dextran Alexa Fluor 647 (10,000 MW; D22914, Invitrogen) was injected in the blastoderm of 1k-stage embryos (∼3 hpf).

30Blastoderm viscosity was measured by micropipette aspiration based on previously established protocols. Briefly, embryos were placed on 3% methylcellulose coated glass coverslips in 1x Danieau’s solution on an inverted Leica SP5 microscope. A fire-polished, passivated (with heat inactivated FBS) micropipette of 35 μm inner diameter, 30° bent, with a spike end (Biomedical Instruments) was inserted into the blastoderm center or margin, just below the EVL. The micropipette movements were controlled by motorized micromanipulators (Eppendorf Transferman, Nk2). Upon insertion of the pipette in the blastoderm, an aspiration pressure of 150 Pa was immediately applied using a Microfluidic Flow Control System Pump (Fluigent, Fluiwell) (with negative pressure ranging from 7-750 Pa, a pressure accuracy of 7 Pa and change rate of 200 Pa.s-1) and the Dikeria micromanipulation software. The value of the applied pressure for the measurements at all stages during blastoderm fluidization/re-solidification (150 Pa) was set according to prior test aspiration experiments, where the applied pressure was continuously increased in a stepwise process (10 Pa / 20 s) until the aspirated tissue started flowing into the pipette. For creep and recovery experiments, pressure was applied until the tissue flew into the pipette at a constant velocity (for ∼3 min, except the cases where the tissue was fluidized and the deformation was too fast) and then pressure was immediately released. Constant tissue flow velocity was similar for aspiration times of 3, 4, 5, and 6 min, and experiments were thus performed with 3min aspiration time (Figure S1A). Images for monitoring the aspiration and relaxation of the tissue were acquired every 500 ms. Viscosity calculations were performed as previously described. Briefly, the tongue length for each time point was measured using a customized Fiji macro, and changes in tongue length during aspiration and relaxation were then plotted over time (Figure 1B). The slope of the aspiration curve at the point of constant flow depends on the viscosity, with R being the radius of the pipette, P the applied pressure and Pc the critical pressure. When the pressure is set to zero during the relaxation, the tissue retracts at a velocity. From the aspiration and retraction rates, viscosity can be calculated as. In case the retraction rates are very low, then the major determinant of viscosity is the aspiration rate.

31Embryos were fixed in 2% paraformaldehyde (PFA) for 4 h at RT. After fixation, embryos were washed in 0.5% Tween-20 (in 1x PBS) overnight at 4°C, dechorionated and washed in 0.5% Tween-20, 0.5% Triton X-100, 0.1 M glycine (in 1 x PBS) for 1 h at RT. Embryos were then incubated in blocking solution (0.5% Tween-20, 0.5% Triton X-100, 1% DMSO, 1% BSA in 1 x PBS) for 3-4 h at RT and then incubated with primary antibody (rabbit E-cadherin anti-zebrafish, 1:200, generated at MPI-CBG) diluted in blocking solution overnight at 4°C. Embryos were subsequently washed 4 × 20 min in 0.5% Tween-20 and incubated with secondary antibody (goat Alexa Fluor 546 anti-rabbit, 1:500, A11010 ThermoFisher Scientific) diluted in blocking solution overnight at 4°C. Finally, embryos were washed 4 × 20 min in 0.5% Tween-20, post-fixed in 4% PFA for 30 min at RT and imaged (Figure S3B).

32In situ hybridization was performed as previously described. Briefly, 4.3 hpf embryos were fixed in 4% PFA overnight and then dehydrated and stored for at least 1 day in 100% Methanol at – 20°C. Embryos were then rehydrated, permeabilized and incubated in hybridization buffer and with the RNA probes overnight. Embryos were washed in serial dilutions of SSC, incubated with an Alkaline phosphatase anti- digoxigenin (DIG) antibody overnight (A2237, Sigma), then washed and stained with a mixture of 4-Nitro blue tetrazolium chloride (NBT) (11383221001, Roche) and 5-brome-4-chloro-indolyl-phosphate (BCIP) (11383221001, Roche). Antisense RNA probes were synthesized from cDNA for ntl, gsc and wnt11, using SP6, T7 or T3 RNA polymerase from mMessage mMachine kits (AM1344, ThermoFisher) with DIG RNA-labeling kit (11 277 073 910, Sigma) (Figure S3D).

33Dechorionated embryos were mounted in 0.5% low melting point agarose (16,520-050, Invitrogen) on a glass bottom dish (P35G-1.5-14-C, MatTek Corporation). Mounted embryos were kept in an incubation chamber at 28.5°C throughout acquisition. Whole embryo single plane bright-field/fluorescence imaging was performed on a Nikon Eclipse inverted wide-field microscope equipped with CFI Plan Fluor 10x/0.3 objective (Nikon) and a fluorescent light source (Lumencor). For high magnification confocal imaging of deep cells to reconstruct connectivity maps, a Zeiss LSM880 inverted microscope, equipped with a Plan-Apochromat 40x/NA 1.2 water-immersion objective (Zeiss), GaAsP and multialkali PMT detectors was used. The Fast Airyscan super-resolution module was used for high magnification confocal imaging of cell-cell contact length fluctuations. For imaging the micropipette aspiration experiments, a Leica SP5 inverted microscope equipped with a resonant scanner and a HC Plan-Apochromat 10x/NA 0.4 objective (Leica) was used. For imaging the immunofluorescence experiments on E-cadherin levels, a Zeiss LSM880 upright microscope, equipped with a Plan-Apochromat 20x/1.0 water immersion DIC objective was used. For imaging the in situ hybridization experiments, an Olympus SZX 12 stereo-microscope, equipped with a QImaging Micropublisher 5.0 camera was used.

34Cell connectivity was defined on 2D confocal sections of the 1st-2nd or 3rd-4th deep-cell layer, in which the nuclei, cell–cell contacts and interstitial fluid accumulations were differentially labeled. The image from the interstitial fluid channel was converted to a binary image and processed by a median filter of 2-pixel width in Fiji, to emphasize the interstitial fluid gaps. Cell nuclei of cells within the same focal plane that had no interstitial fluid between them were considered as contacting cells. By using the multi-point tool in Fiji, the coordinates of each nuclei and their contacting-neighbors were marked and plotted using standard python package matplotlib. Due to the manual identification of connections, the connectivity of each image was corrected twice by overlaying the connectivity map on top of the blastoderm image. The rigidity analysis to identify floppy and rigid areas, Giant and 2nd Giant clusters, and cluster size distributions were performed on the connectivity maps using a python version of the pebble algorithm pebble.py (Key resources table).

35In this section we present the precise method used to generate the lattices that serve as a model of the cell contact networks analyzed in this study.

36We worked in general with random triangular lattices,, embedded in a flat space (see Figure S1H). Throughout the text, we will refer indistinctly to them as networks or lattices. is the set of nodes and is

37the set of links connecting nodes of. We call average connectivity the average number of edges connected to a given node, defined as:A fully connected triangular lattice contains all possible links between the nodes, and is typically used to model packing problems. For triangular lattices, in the limit,. For finite lattices, the low number of links connecting the nodes in the boundary implies a reduction of this maximum connectivity,, and the specific number will depend on the exact value of and the node arrangements. Considering that nodes are arranged in a square with sides (corresponding to the shape of experimental networks), can be very well approximated as:Square geometries are commonly used in the study of rigidity percolation, where finite size effects are well understood. Our data fits satisfactorily with the square-geometry assumption, although an exception is found in the networks belonging to the margins, since, in these networks, nodes are disposed along irregular geometries. We therefore focused our rigidity analysis to the networks located in the center, since for them we could construct model networks with a unified criterion, thereby reducing the potential effects of the boundaries and enabling fair comparison. As discussed in the main text, we restricted the rigidity analysis to 2D networks, which are more tractable both experimentally and computationally. However, the 2D connectivity of the embryo was identical at different z-slices (Figures S1K and S1K’), meaning that no anisotropy in the packing is expected, and that one expects to be able to predict 3D average connectivity from these 2D sections. To test this more directly, we have also measured the local average number of contacts per cell both in 2D and 3D for a subset of randomly selected cells (see Figure S1L’ for further details on the methods of measurement), and confirmed that we also observed a drastic decrease of average connectivity in 3D measurements at the fluidization point (Figure S1L). As expected for triangular lattices, the average number of contacts in 3D was around the double of the average number of contacts in 2D at all time points. Further extension to 3D would require a closer quantification of the exact structure of 3D connectivity, and in particular the type of packing of the blastoderm cells in 3D, which would impact on the threshold. In addition, in order to compare systems with different cellular numbers (which increases as a function of time due to division in the blastoderm), we worked instead with the normalized connectivity of a given finite lattice defined as: corresponds to the probability that a link exists in the triangular lattice of a given size. In Figure S1H we provide an example of a fully connected lattice and a lattice with average connectivity 3.6.

38The methods described here are implemented in the python package Lattice.py. In brief, in silico random triangular lattices are constructed based on two parameters and. The construction steps are:Start from a fully connected triangular lattice made of nodes arranged uniformly in a square with links of equal length Remove links at random until the average connectivity value is achieved.

39Start from a fully connected triangular lattice made of nodes arranged uniformly in a square with links of equal length

40Remove links at random until the average connectivity value is achieved.

41In the study of Chk1 overexpression (oe) tissues, we needed to generate spatially correlated lattices, implying that the presence/absence of links is not purely random throughout the tissue, but spatially correlated. In silico correlated triangular lattices are constructed following four parameters:,, and. is the correlation length, which we fixed to be a.u. is the parameter accounting for the strength of spatial correlations. Before starting the description of the algorithm, we defined the vicinity of a link as all the links whose middle point is located at a distance of the middle point of. The construction steps are:Start from a fully connected triangular lattice made of nodes uniformly arranged in a square with links of equal lengthChoose a link at random and remove it. Call this link.With probability, choose a link in the vicinity of the previously chosen link and repeat the step. With probability go to the previous step.If the links in the vicinity of have been all removed, randomly chose a link from those that have been removed in the vicinity of, call it and go to the previous step.keep the process running until the average connectivity value is achieved.

42Start from a fully connected triangular lattice made of nodes uniformly arranged in a square with links of equal length

43Choose a link at random and remove it. Call this link.

44With probability, choose a link in the vicinity of the previously chosen link and repeat the step. With probability go to the previous step.

45If the links in the vicinity of have been all removed, randomly chose a link from those that have been removed in the vicinity of, call it and go to the previous step.

46keep the process running until the average connectivity value is achieved.

47Effectively, this process removes links sequentially creating holes in some areas of the network while leaving other areas highly dense. Holes represent sets of links that have been sequentially removed, and its average size follows a Poisson statistics with average hole size of order. In Figure S5F the method is schematically outlined. As shown in Figures 5C and S5E, this reproduced well the network statistics observed in Chk1 oe embryos, while random networks reproduced well the network statistics of WT embryos. From a biological perspective, such correlations are expected from heterogeneities in cell divisions, which are linked to cell contact loss. Given that cells mix relatively little during early morphogenesis, this means that some regions of the embryo will have clusters of “delayed” cells sharing a common ancestor, and which will remove their bonds in a correlated manner, given the experimentally determined correlation between cell division and contact loss.

48The main theoretical concept underlying the present study is rigidity. Rigidity is a property related only to the connectivity pattern of the nodes (i.e., it is a topological property), but it has direct impact on the material response of the net. Rigidity is a property of a subgraph of a given network, and this is why we refer to rigid clusters. Following, a rigid cluster is defined as “A collection of nodes form a rigid cluster when no relative motion within that cluster can be achieved without a cost in energy.”

49In that context, energy costs arise from the stretching or compression of links after a node rearrangement. We will consider only maximal rigid clusters: Given a lattice, a subgraph is a maximal rigid cluster if i) it is rigid, according to the above definition, and ii) there is no subgraph of that is rigid and contains properly. In other words, there is no such that and is rigid. In addition, we consider that clusters are induced subgraphs of a given lattice. An induced subgraph is a subgraph if all the connections between the nodes present in are also present in. Therefore, we take this working definition for our purposes and, for the sake of clarity, adapt all the theory to the framework of 2D triangular lattices:

50Given a triangular lattice over a plane, a rigid cluster is a maximally induced subgraph by which any movement over the plane of any proper subset of nodes of it will unavoidably result into a stretching or compression of at least one link.

51Importantly, the above definition does not exclude i) the presence of more than one rigid cluster in a network – in general, there is a number of them– and ii) the fact that a node participates in more than one rigid cluster. Rigid clusters as defined above are identified using the Pebble game algorithm.

52As pointed out above, a triangular lattice may contain an arbitrary number of rigid clusters, constrained only by the fact that the minimum number of nodes for a rigid cluster is, as a triangle configuration is the minimal rigid structure. If we name a given rigid cluster, we define the set of all rigid clusters of as R = {R1(V1,E1), …,Rk(Vk,Ek)}. Given a cluster, its size is the number of nodes it contains

53A particularly relevant element of the set R is the largest or giant rigid cluster,, of the lattice, due its prominent role in the material response of the network. Since the size of the may diverge with the size of the lattice, the most relevant quantity is the relative size of the giant cluster relative to the total number of nodes of the whole lattice,:The second largest rigid cluster, which we refer to as the second giant cluster (also referred to simply as second cluster in the text), also plays an important role in our analysis (in particular to identify spatial correlations and heterogeneities in connectivity patterns, as discussed in Figures 5D, 5D’, 5E, S5G, S5G’, and S5H), but it is less relevant in order to understand the material response of the system.

54At low/mid values of connectivity, the expected fraction of nodes belonging to the largest rigid cluster is negligible. However, when the connectivity crosses a critical value, the size of the largest cluster experiences an abrupt jump from the previous negligible value to almost all the network. This phenomenon is a second order phase transition and has, as an order parameter, the size of the largest rigid cluster and, as a control parameter the connectivity. This phase transition has important consequences on the rheological properties of the network.

55Specifically, in the limit of infinite triangular lattices, the value of tends to for any value. For, experiences a sharp growth close to until it saturates at. We call the critical point of the rigidity phase transition for triangular networks, also known as the isostatic point for rigidity. At the critical point, we observed the emergence of a rigid cluster that spans a finite fraction of the lattice in a second order phase transition. Standard analysis of signatures of criticality are based on the study of the behavior of the correlation length close to the critical point, showing a scaling behavior, with.

56Finding traces of criticality within our data using this classical metric can be problematic: To grasp this kind of behavior, one needs to deal with very large systems arbitrarily close to the critical point. The number of networks available and their intrinsic small sizes makes this approach impracticable. However, there are other traces of criticality in percolation phenomena that revealed more robust when looking at real systems, such as the cluster size distribution. For example, in standard percolation in random graphs, cluster size distribution near the critical point displays a clearly defined power-law pattern with exponent. As a consequence, the average cluster size – excluding the giant cluster in the computations – diverges at the critical point. For the specific study of rigidity transitions, these markers of criticality were also explored, but in slightly different network processes. We took this approach in order to elucidate whether the experimental data had traces of criticality. Numerical explorations with square lattices with side size, containing nodes show that, close to the critical point (found around, for lattices of this size), the probability distribution of rigid cluster sizes other than the,, follows a well defined power law:with, fitted with the method of maximum likelihood estimation, as described in (see Figure S3K). Furthermore, computing the rigid cluster size distribution for smaller graphs (of the same average size as experimental graphs), predicted similar exponents although with a finite-size correction, both features providing excellent fit to the data (Figures 3G and S3N).

57This critical exponent implies, as opposed to standard percolation in random graphs, that the average cluster size converges in a value independent of the network size and connectivity around for the interval of connectivities under study, whereas at the critical point, the variance of the cluster size,, is expected to diverge with the size of the lattice (see Figure S3L). In particular, given the value of the exponent, one expects to grow as a function of the system size,, as:where, which agrees with our computational findings (Figure S3M). Importantly, turning to experimental cell-contact network revealed a behavior fully consistent with these results, with both the position and amplitude of the peak of well-fitted by the rigidity percolation model in the absence of adjustable parameters.

58It should be noted that previous works had also computed cluster size exponents, for slightly different phenomena, and found slightly different values of for non-deformable adhesive spheres and subcritical or supercritical regimes of self-organized, less stochastic network growth.

59An important insight of the rigidity percolation is its prediction on the rheological properties of the tissue. In particular, at the rigidity percolation critical point, networks acquire a non-zero elastic modulus upon deformation (matching the shear response of sphere packings). Below the critical point of rigidity, small deformations of the network may be performed at no energy cost. On the contrary, beyond the critical point of rigidity, the network as a whole, starts exerting non-negligible force against any deformation. Similarly, work on the influence of weak attractive forces on the rheology of colloidal suspension predict a divergence of viscosity at a critical particle fraction, after which elasticity should be observed. In our experimental system, however, the tissue behaves as a viscous fluid even for very high connectivity values (Figure 1F, or alternatively very high cellular fractions, Figure S1E). In fact, beyond the rigidity percolation critical point, viscosity varies relatively little with further increase of connectivity (Figure 1F), so that even fully connected graphs display viscous rheology under micropipette aspiration. Given that such transition from elasticity to viscosity is frequent for biological material undergoing turnover or undergoing passive or active topological re-arrangements (T1 transitions), and that re-arrangements and finite contact times were observed in our system, we explored the idea that incorporating in a toy model turnover would result in viscosity behaving similarly as elasticity (i.e., below the percolation transition, neither elastic or viscous network can resist deformation while above both can). We start by recapitulating the response of a purely elastic network as previously described, using a large triangular lattice with the bottom layer fixed in space with links as springs of elastic constant. Given that is the rest length of the spring, force balance on a given (overdamped) node at position leads to the following differential equation:where is the damping parameter of the spring, is an index summing over all neighbors of, the length of the link connecting nodes and and the unit vector associated to this link. The nodes therefore define a system of coupled differential equations determining the evolution of the network under potential external forces or deformations. To infer the linear response of the lattice against deformations, we complement this with two boundary conditions: i) we quasi-statically apply a small lateral displacement of the top layer by (see Figure 1H) small enough to grasp the linear response against this deformation performed by the network and ii) we fix the position of the bottom layer, to calculate the force exerted by this layer, against the imposed displacement at the top. This is similar to rheological assays with controlled strain.

60In rigidity percolation, elastic networks cannot transmit forces below the critical point. Numerical experiments testing against confirmed this, showing that displays negligible values until (Figure S1J). For connectivity values above the rigidity percolation threshold, the transmitted force (i.e., elastic modulus) increases linearly with, again as expected from rigidity percolation in networks, as well as sphere packing under shear.

61We now include the dynamic nature of cell-cell contacts to a permanent energy dissipation that will provide the network a viscous behavior. In general, viscosity can be defined as:where is the force exerted over the system to deform it and the velocity of the part of the network where this force is applied, as a response to it. In networks with turnover, energy is released spontaneously, enabling the network to be deformed irreversibly. The simplest way to implement such energy dissipation is to consider breaking and healing of bonds at a constant rate in our system, which will give rise to viscous networks as measured experimentally. In this framework, a bond with rest length with actual length may break and heal again with timescale, resetting the rest length to. In the case the network is subject to some external force/deformation, these dynamics will result in an irreversible deformation. In the case the network is subject to a constant force e.g., pushing up the top layer, the network will reach a regime where it will flow at constant speed depending on. Let us analyze this situation and how it projects into viscosity. Now assume that energy is dissipated due to a relaxation of the rest length of the springs occurring stochastically with an average time window of. In this case, there is a constant,, that determines the relation between velocity and, leading to Equation 2 to read Since, one can rewrite the above expression as, which relates the relaxation time and the viscosity, providing therefore a very simple framework to explain viscous behaviors observed in our experiments. A crucial point is that, at the linear order, we expect the response to be similar to elastic networks with respect to the critical point, which we proceeded to check. To that end, we followed the following modeling criteria:We consider random triangular lattices with increasing connectivity and we deformed them by displacing all the nodes of the upper layer a distance (imposed shear), keeping the coordinate constant for this layer.The total number of integration time steps was chosen such that a fully connected, elastic lattice of the same size than the ones used by the simulation reached a stationary state in terms of force propagation through the whole network under the displacement of of the nodes belonging to the upper layer.

62We consider random triangular lattices with increasing connectivity and we deformed them by displacing all the nodes of the upper layer a distance (imposed shear), keeping the coordinate constant for this layer.

63The total number of integration time steps was chosen such that a fully connected, elastic lattice of the same size than the ones used by the simulation reached a stationary state in terms of force propagation through the whole network under the displacement of of the nodes belonging to the upper layer.

64We performed numerical tests under the above considered conditions. We generated random triangular lattices of size,,,, and we performed a permanent displacement of the top layer of. ranged from to in 20 equidistant connectivity intervals, and we performed replicas per connectivity value. Assuming that the equilibration time of the elastic network is, we considered three cases: i) elastic network ii) A viscous network whose links have an average lifetime of and iii) a viscous network whose links have an average lifetime of (Figures 1H and S1J). Calibration was performed under the criteria that force distribution should reach a stationary pattern in the fully connected elastic network of this size. Under that criteria, we found that was enough.

65To quantify the adhesive properties of the cells, we compute the ratio between cell-cell and cell-fluid surface tensions ( and, respectively) following the methodology developed in. Let us suppose two cells, whose geometry in isolation is spherical with identical radius. When we put them in contact, the balance of surface tensions will create a circular region of contact, and the contact point will define an angle between the two cell membranes and the fluid (Figure 4A) of the main text for a schematic. Following Maître et al., 2015, the area of contact between cells is:In turn, the area of the cells in contact with the fluid is related to the contact angle as:According to these two equations, the energy of these two cells in contact is given by:From the above relation, it can be shown that the relation between surface tensions when the system reached an energy minimum is given by:Therefore, by defining the relative cell-cell tension (which decreases for increasing cell adhesion) as:one has a single parameter depicting the relation between surface tensions, which can be extracted from experimental data by computing the angle of contact between cells at the fluid interface. This analysis shows that relative cell-cell tension (resp. cell-cell adhesion) decreases (resp. increases) at the fluidization point (Figure 4B and B’), and displayed a strong anti-correlation with tissue viscosity across time and conditions (Figure 4D).

66To relate these changes to the changes of cell connectivity observed from blastoderm network reconstruction, we generalized this doublet model to configurations of arbitrary cell numbers. In particular, we used a two-dimensional toy model of 4 cells with equal volumes and tension, organized spatially into a rhombus with 4 links. This mirrored one of the simplest floppy topologies sketched in Figure 4A. We started from non-adhesive spheres (, and gradually decreased cell-cell tension (increased adhesion). This caused an expansion of cell-cell doublet contacts, which above a critical threshold in relative cell-cell tension merged to form tri-cellular vertices, triggering the spontaneous creation of a cell contact that makes the system to transit from floppy to rigid, in topological terms (Figure 4A). The trigger of rigidity is found at around –, although this value depends on the exact initial configuration used, and displays hysteresis upon increasing again above the threshold. Finally, we also considered disordered arrangement of cells, and found that the threshold value toward rigidity was located around the same ranges of (Figure S4E). Simulations to compute equilibrium arrangements of cells have been performed using the free available software ‘Surface Evolver’.

67All acquisition data were processed using Fiji (NIH) and/or Imaris 9.0.

68Cell connectivity (< C > ) was calculated in each confocal section as total number of contacts (defined as described in connectivity map reconstruction) divided by the total number of cells in the image. Normalized connectivity (< k > ) was calculated in each confocal section as connectivity C divided by the maximum potential connectivity (kmax) (computed as described in Supplementary Information Text). 3D number of contacts per cell was calculated in Imaris by randomly choosing blastoderm cells, followed by manual identification of the cell-cell contacts using the “Ortho-slice” and “measurement” tools (Figures S1L and S1L’).

69The mean squared relative displacement (MSRD) was calculated from individual 3D cell trajectories following the nuclei position in Imaris of two cells (CellA, CellB) that were in contact at the first time point. The relative xyz position between the CellA and CellB was defined as:where,,. The MSRD was defined as:with t0 being the first time point and t1 3 min later (Figures S3I and S3I’). MSRD for the different time periods was defined as t1 = −30 and t0 = −60 for t-60 min, t1 = 0 and t0 = −30 for t-30 min, t1 = 30 and t0 = 0 for t0 min, t1 = 60 and t0 = 30 for t30 min and t1 = 90 and t0 = 60 for t 60min (Figures S1B and S1B’).

70Cell fraction (cf.) was calculated as:ff was calculated from 2D confocal sections at the 1st-2nd deep cell layers, where the interstitial fluid channel (labeled with Alex Fluor 647 Dextran) was first converted into a binary image and then processed by a median filter of 2-pixel radius in Fiji. A signal intensity histogram analysis of the whole binary image was performed in Fiji (range of 256 different signal intensities, with 256 being the maximum), and the ff was obtained by the average gray value of the histogram normalized to the max signal intensity of the histogram (256) (Figure S1E).

71Cell shape index (q) analysis was calculated from 2D confocal sections at the 1st-2nd deep cell layers, where the cell perimeter P and cell area A were obtained using the freehand selection tool to draw the cell perimeter and q was defined as previously described with (Figure S1C). Cell area analysis shown in Figure S5C was also performed using the freehand selection tool to draw the cell perimeter in Fiji. Cell area variability shown in Figure S5D was expressed as the coefficient of variation of the cell area measurements shown in Figure S5C.

72E-cadherin expression was quantified by the fluorescent intensity of E-cadherin levels as judged by antibody staining. Control-morphant and e-cadherin-morphant embryos were processed for immunofluorescence within the same tube. The average fluorescence intensity was calculated in Fiji within an ROI covering the central deep-cell blastoderm (Figure S3C).

73The cell-cell contact angle, θe, was calculated in degrees from 2D confocal sections using the angle tool in Fiji, as indicated in Figure 4B, and then converted to radians. The cell-cell tension α was then calculated as:The spatial heterogeneity in cell connectivity (khet) and fraction of dividing cells (fdhet) (Figures 5A and 5A’) was analyzed by computing the variance of the two parameters between the quadrants of the 2D confocal sections used for the tissue rigidity analysis. For khet, in each quadrant (Area I, II, III, IV), the normalized connectivity was calculated as described above (kI, kII, kIII, kIV) and then the variance was calculated between kI, kII, kIII, kIV. For fdhet, in each quadrant the fraction of dividing cells was calculated as the number of mitotic nuclei (excluding cytokinesis stage) divided by the total number of nuclei in the plane (fdI, fdII, fdIII, fdIV). In order to compare the variance in fd between WT and Chk1 overexpressing (oe) embryos displaying different number of cells, the standard deviation between the quadrants was first normalized using a binomial test. The probability q of a cell dividing was approached as, with n, total number of mitotic nuclei and N, total number of nuclei. This defines the null model underlying the statistical test. The probability that a dividing cell is in a given quadrant is 1/4. The probability that Ni cells in quadrant i from which ni are dividing, p(Ni, ni) is then given by:Here, is a binomial distribution with parameter 1/4 and number of trials N. is a binomial distribution with parameter q and number of trials Ni. The expected average number of dividing cells per quadrant is then. The expected standard deviation of dividing cells per quadrant is:The absolute value of the z-score for each quadrant and the number of dividing cells is:

74, with Ki, the observed number of dividing cells in quadrant Qi. The spatial heterogeneity coefficient, SHC, is then, plotted in Figure 4B.

75The fraction of dividing cells shown in Figure S5A’ was calculated as the number of mitotic nuclei (except cytokinesis phase), divided by the total number of nuclei in the 2D confocal sections. The relationship between connectivity and fraction of dividing cells was analyzed during the meta-synchronous cell cleavage cycles (−60 to 0 min for WT embryos, −60 to 60 min for Chk1 oe embryos) in higher magnification 2D confocal sections by following every 5 min a randomly selected cell and its direct neighbors (presence of a cell-cell contact). The increase in the fraction of dividing cells and the reduction in the relative number of contacts every 10 min was plotted in Figure S5B.

76The quantification of differences in the viscosity fold-change between WT and Chk1 oe embryos over time was calculated from viscosity measurement experiments where all embryos originate from the same batch. The relative change in the viscosity value was calculated as, with η being the viscosity values and the time interval between t1 and t2 being < 10 min. From each batch the percentage of embryos exhibiting 0-1, 1-5, 5-10 and > 10 was calculated for each time point (Figure S4I).

77The variability in the viscosity was expressed as the coefficient of variation for each time point of viscosity measurements for the same number of embryos between WT and Chk1 oe embryos (Figure S5J).

78The robustness viscosity factor was expressed as the inverse of the coefficient of variation (Figure 5G).

79Cell-cell contact length fluctuations were calculated from 2D confocal time series acquired using the Fast Airyscan mode (10 s interval, 10 min time lapse duration). For each time point, the contact length was measured using the free-hand line tool in Fiji (Figure S5K). Then, contact lengths were normalized to the average of the contact lengths during the 10 min duration and the frequency distribution (in percentage) was calculated in GraphPad prism (contact lengths were binned in 0.1 intervals) (Figure S5K’). The contact length fluctuations in Figure S5K’’ was expressed as the coefficient of variation,, for each contact during the 10 min measurement, and the average of CV of 10 contacts per time point and condition was plotted. The kymographs shown in Figure S5K were generated in Fiji, using the KymoResliceWide plugin.

80The statistical analyses were performed with GraphPad Prism 6.0. Statistical details of experiments are reported in the figures and figure legends. Sample size is reported in the figure legends, and no statistical test was used to determine sample size. The biological replicate is defined as the number of embryos or independent batches of embryos, as stated in the figure legends. No inclusion/exclusion or randomization criteria were used and all analyzed samples are included. Unless differently stated in the figure legends, the graphs show mean ± sem and the error bars are calculated and shown based on the number of cells or embryos, as indicated. The statistical test used to assess significance is stated in the figure legends and was chosen after testing each group with the normality distribution test D’Agostino. For comparing two groups, a two- tailed Student’s t test for parametric distributions and a Mann–Whitney U-test for non-parametric distributions were used. To compare more than two groups, either an analysis of variance (ANOVA) or Kruskal–Wallis test for parametric and non-parametric distributions, respectively, was used. No blind allocations were used during the experiments or in the analysis. At least more than three independent experiments were performed except for the data shown in Figure S3C (embryos from two different embryo batches, in which an independent experiment was defined as the embryo batch). This information is also stated in the figure legends.

81To obtain the universal curve accounting for the expected size of the GC in random triangular lattices along connectivity values, we first computed the average size of all experimental networks (except the Chk1 overexpression networks, which we treat separately as discussed below), leading to. We then generated model lattices for equidistant connectivity values from to of size. For every connectivity value, we performed 50 replicas and computed the average and standard deviation of the GC.

82To compute the model curve appearing in Figure 1G’, we proceeded exactly as above with different lattice sizes. In particular, we considered the following network sizes respectively.

83To obtain the curve accounting for the expected sizes of the clusters other than the GC in real networks along connectivity values, we considered equidistant connectivity values from to. We then created, for each connectivity point, 50 lattices with the same criteria as above, computed the overall cluster size average and variance and compared them to the curve obtained for real networks.

84To produce the part corresponding to real data of the plot shown in Figures 3F and 3G, we first considered that, for networks of size, the peak in the – which we take as the indicator of the critical point – was located around. The cluster size analysis near criticality was performed by collecting the cluster sizes other than the GC for all the networks whose laid in the interval and then computing the statistics. This interval ensured certain degree of resolution as well as certain volume of data to perform statistics. According to that criteria, we found that of the studied experimental networks were in the critical regime. We then computed the same metrics with the same criteria in computational networks of the same average size.

85Cluster size distribution offers a detailed map of the potential heterogeneities and deviations from pure randomness in the connectivity patterns. To avoid confusion, it is worth to emphasize that here we are considering the whole range of connectivity values, not only those near the critical point. As we see in the results, WT cluster size distribution can be fairly approximated using random lattices as models. However, for Chk1 overexpressing (oe) embryos the case is different: The cluster size distribution departs clearly from the expected behavior if the underlying null model are random networks. We needed to use, instead, a model lattice that creates spatial correlations, as the one defined above. Below we detail the numerical experiments performed for the cluster size analysis.

86To perform reliable statistics in order to derive the cluster size distribution, we needed to be more detailed than in the above computations, given the small nature of the system and the fine-grained observable we are dealing with. To that end, we refined the generation of the null model for lattice structure.

87We created uncorrelated lattices, which is the exact number of WT real networks in the sample. Lattice sizes were picked at random following the lattice size distribution obtained from real data. The average connectivity of each generated random lattice was chosen at random uniformly within the interval and.

88With the above methodology, we create a null model for the topology of (random) lattices belonging to the WT set. We performed 30 replicas of the whole process and plotted the results in Figures 5C and S5E.

89For the modeling of cluster size distribution of the Chk1 oe data (Figure 5C) of the main text, we used the method for generating correlated lattices with data parameters, defined as above, and as the single free parameter. We created, which is the exact number of Chk1 oe real networks in the sample. Lattices sizes were picked at random following the lattice size distribution obtained from real data. The average connectivity of each generated random lattice was chosen at random uniformly within the interval and. With the above methodology, we create a null model for the topology of (correlated) lattices belonging to the WT set. We performed 30 replicas of the whole process and plotted the results in Figure S5E.

90To find the optimal fit of the parameter between real data and simulated results we used the method of minimization of the Jensen-Shannon divergence between real and numerical data. The global optimum was found in, implying a spatial correlation in the distribution of cluster sizes of a few cell length, which is quite large given the finite size of the experimental networks, providing an explanation for the much less uniform viscosities measured experimentally in Chk1 oe embryos during the studied period.

91This analysis was performed identically for the WT and Chk1 oe datasets. To check the statistical relevance of the potential deviations of the size of the second largest cluster, for each real network, we generated 100 replica networks with approximately the same size and average connectivity and, from this population, we computed the expected value for the second largest cluster. We then aggregated the real networks in their respective different time points and computed the average size and standard deviation of the second cluster. We did the same operation using the expected second cluster sizes we computed before. Consistent with the deviation from the expected cluster size distribution observed in Chk1 oe embryos, we observe that, in general, the size of the second cluster is bigger than the expected cluster size (Figures S5G and S5G’).