The principles underlying the biomechanics of morphogenesis are largely unknown. Epiboly is an essential embryonic event in which three tissues coordinate to direct the expansion of the blastoderm. How and where forces are generated during epiboly, and how these are globally coupled remains elusive. Here we developed a method, hydrodynamic regression (HR), to infer 3D pressure fields, mechanical power, and cortical surface tension profiles. HR is based on velocity measurements retrieved from 2D+T microscopy and their hydrodynamic modeling. We applied HR to identify biomechanically active structures and changes in cortex local tension during epiboly in zebrafish. Based on our results, we propose a novel physical description for epiboly, where tissue movements are directed by a polarized gradient of cortical tension. We found that this gradient relies on local contractile forces at the cortex, differences in elastic properties between cortex components and the passive transmission of forces within the yolk cell. All in all, our work identifies a novel way to physically regulate concerted cellular movements that might be instrumental for the mechanical control of many morphogenetic processes.
Hydrodynamic regression is used to infer mechanical parameters from microscopy data obtained during epiboly in zebrafish. This novel approach identifies a polarized gradient of cortical tension as the driver of tissue movements.
Hydrodynamic regression infers 3D pressure fields, mechanical poser, and cortical surface tension profiles during epiboly in zebrafish.
Epiboly tissue movements are directed by a polarized gradient of cortical tension.
This gradient relies on local contractile forces at the cortex, differences in elastic properties between cortex components and passive force transmission within the yolk cell.
The construction of complex functional structures during morphogenesis results from the unfolding of the developmental program of individual cells. Importantly, morphogenesis does not only rely on a tight regulation of gene expression and cell communication but also depends on fundamental physical principles (Kiehart et al, 2000; Desprat et al, 2008; Pouille et al, 2009; Rauzi & Lecuit, 2009; Aigouy et al, 2010; Lepage & Bruce, 2010; Martin, 2010; Grill, 2011; Lecuit et al, 2011). The stereotyped tissue movements involved and the embryo optical accessibility makes zebrafish epiboly ideal for exploring mechanical coordination during morphogenesis (Concha & Adams, 1998; Keller et al, 2003; Behrndt et al, 2012).
Epiboly is a conserved early embryonic morphogenetic movement that involves a striking physical reorganization characterized by the thinning and spreading of cell layers. Epiboly initiates when the enveloping layer (EVL), an external epithelial sheet, the deep cells—DCs, a spherical cap of blastomers centered on the animal pole of the embryo over the yolk cell, and the external yolk syncytial layer (E‐YSL), a specialized region at the yolk surface positioned ahead at the EVL, undergo cortical vegetalward movements. These are associated with the shortening of the rest of the yolk surface (yolk cytoplasmic layer—YCL) and the doming up of the internal yolk. Epiboly ends by the closure of the blastoderm at the vegetal pole covering the yolk cell (Kimmel et al, 1995; Solnica‐Krezel, 2006; Rohde & Heisenberg, 2007; Fig 1A).
Biomechanical studies experience a major limitation, that is the difficulty to directly measure forces in vivo. Several methods have been developed to measure the mechanical properties of living cells. So far, however, they cannot be readily applied to whole developing organisms. Exceptions are the use of laser microsurgery (Toyama et al, 2008) and the use of genetically engineered mechanical biosensors (Grashoff et al, 2010). Alas, such techniques are intrusive, fall short of quantitative characterization [e.g., laser cuts cannot distinguish active and passive responses (Rodriguez‐Diaz et al, 2008)] and cannot provide a global biomechanical representation of morphogenesis.
We have developed a novel non‐intrusive method that we refer to as hydrodynamic regression (HR), which can accurately and globally estimate the spatio‐temporal profile of mechanical power and cortical tension in living organisms. It relies on considering that in deforming tissues during morphogenesis forces are balanced and stresses are continuous at the cortex/fluid interface. This physical principle (boundary condition) has been shown to apply at the single‐cell embryo level (Niwayama et al, 2011) and on multicellular epithelia in gastrulating embryos (He et al, 2014). HR is mesoscopic and, although it cannot infer the mechanical properties of single cells, it helps identifying the regions where mechanical stress and power concentrate during tissue rearrangements. Importantly, it does not demand restrictive assumptions either on the mechanical nature of the cortex under stretch or on the physical properties of the embedding fluids. Thus, it can be applied to any morphogenetic process providing two conditions are met: (i) The geometry/shape of the tissues involved can be extracted from images and modeled and (ii) the flows around the cortex occur at low Reynolds number, where inertia is negligible, which is the case for most biological processes (Niwayama et al, 2011; He et al, 2014).
We applied HR to study epiboly global biomechanics and estimated the topography and temporal evolution of mechanical power and surface stresses of the embryo. These estimations were validated and confirmed by laser microsurgery, microrheology, and atomic force microscopy (AFM). In first place, we concluded that HR is suitable to explore the biomechanics of morphogenetic processes driven by the reorganization of cells or sheets of cells of finite width surrounded by fluids. In countless models, this may most probably be the case. Secondly, the HR analyses and in silico modeling let us suggest that epiboly mechanics and kinematics could be simply explained as a result of the integration of the contractile activity of the E‐YSL and the differential elastic properties of the EVL and the yolk surface. A tension imbalance will build up across the opposite sides of the E‐YSL during epiboly, enabling it to proceed (at least from 50% onwards). These novel predictions were tested in silico and experimentally confirmed by interfering in Misshapen expression [a regulator of myosin contractility active during epiboly progression (Köppen et al, 2006)] and reducing E‐YSL contractility. Our new biomechanical model, in contrast to previously suggested models (Behrndt et al, 2012), accounts for all known features of epiboly kinematics.
Global kinematics of zebrafish epiboly
To perform an analysis of epiboly biomechanics, it is essential to define the kinematic properties of the cortical vegetalward movements of the EVL, DCs, and E‐YSL leading to embryo closure (see Fig 1A). Experimental kymographs of the embryo surface along the animal to vegetal axis (φ axis) were extracted from 2D time‐lapse meridional sections, perpendicular to the dorso‐ventral axis (Fig 1B and Video EV1). They showed a sustained displacement of the EVL/yolk margin position φ*, whose speed v* continuously increased and that could be fitted, after crossing the equator (v0), to the function vφ = v0/sin φ* (Fig 1C and see Appendix Note S1). Kymographs also showed that the local surface speed decreased exponentially from the EVL/yolk margin toward the vegetal pole (Fig 1D). This decay outlined two domains, a proximal one ahead of the EVL (~E‐YSL) undergoing progressive constriction (see Appendix Fig S1) and a mostly static distal area.
The proximal E‐YSL constitutes a morphologically distinguishable region with a highly convoluted membrane (Betchaku & Trinkaus, 1978; Cheng et al, 2004; Fig 1E). In the E‐YSL, actin was conscripted within and beneath wrinkles (Fig 1F) underlined by nestled mitochondria (see also Popgeorgiev et al, 2011; Appendix Fig S2). Importantly, the width of the E‐YSL (δφ) and the velocity of the EVL/yolk margin v* were inversely paired (v* ~1/δφ) at all times. This pairing was just qualitative up to reaching the equator (v0), while quantitatively match theoretical kymographs after 50% (Fig 1C).
In parallel to these cortical (EVL, DCs, and E‐YSL) movements, we observed that the internal yolk granules undergo stereotyped motions. They spin symmetrically around a transverse annular axis creating a torus of vortices (Fig 2A and Video EV2). As epiboly progressed, the centroids of these vortices move gradually toward the cortex and the vegetal pole (Fig 2B–D). Importantly, they parallel the displacement of the edge of the EVL along the φ axis toward the vegetal pole, which after 50% epiboly showed a constant velocity of 200 μm/h (z‐axis) (Fig 2E).
In summary, during epiboly, the EVL leading edge movements, the shrinking and displacement of the E‐YSL, and the stereotyped vortices of the yolk granules are remarkably coordinated and seem mechanically coupled.
Development of a non‐invasive method, hydrodynamic regression, to determine mechanical stress and power patterns
Although some information on the power sources for epiboly is available (Cheng et al, 2004; Köppen et al, 2006), how the different elements involved in tissue spreading are mechanically coordinated remains unexplored. To study this coupling, we developed a mesoscopic method able to dynamically infer cortex stress patterns and mechanical power maps of general applicability (hydrodynamic regression—HR).
The cortex in a biological system can be mesoscopically described as a continuous quasi‐2D structure under tension embedded in viscous fluids. As such, we reasoned that the local differences in surface tension when the cortex expands or contracts could be estimated from the shear stress at its interface with the fluids. These stress values could themselves be indirectly inferred by analyzing the flows in the surrounding media (Fig 3A). Therefore, by analyzing the velocity fields of embedding fluids it should be in principle possible to infer the overall spatio‐temporal distribution of the 2D out‐of‐equilibrium cortical tension, τ (see below). Theoretical considerations and derivations are fully detailed in Appendix Note S2.
To implement HR, we have first to determine the velocity fields of the media surrounding the tissue/cortex. This can be done by employing bright field transmission or two‐photon confocal microscopies live and estimating tissue kinematics by particle image velocimetry (PIV) (Supatto et al, 2005). In a second step, we model the fluid motion by distributing Ne force dipoles (Stokeslet pairs) of variable strength βe perpendicular to the cortex according to its particular geometry at a distance δ away from each other (Fig 3B). These models assume that in developmental processes, fluids at small strain rates behave as incompressible Newtonian fluids at low Reynolds number (Cartwright et al, 2009). Biological fluid dynamics can then be described by a linear continuity equation, ensuring volume conservation, and the linear Stokes equation (1) (2)with ν being the shear viscosity, V the velocity field, Δ the Laplace operator, grad the gradient, and P the mechanical pressure. The linear combination of the fluid flows generated by the force dipoles provides a realistic model of the 3D velocity fields (strain rates). Third, we used regression analysis to fit the analytical models to the velocity fields. At each time point, a Monte Carlo scheme was used to determine the set of βe minimizing the mean square error (MSE) between experimental and theoretical fields. The mechanical pressure was determined from these same values. Last, considering the fluid strain rates and the mechanical pressure, we inferred the fluid shear stress, the local values of τ, the average surface tension, and the cortical mechanical power Π, that is, the rate of mechanical energy per unit of time (Fig 3C). Performing HR over time gave access to the spatio‐temporal evolution of both, cortical stresses and mechanical power maps.
The pipeline used to infer the mechanics of the cortex and its surrounding fluids by HR is shown in Fig 3D. The workflow and protocols involved are described in Appendix Note S2, and a flow simulation software suite with step by step documentation is provided in Appendix Note S3.
In principle, HR can be applied to different scenarios, materials, and geometries. To verify this versatility, we applied the method to diverse structures with known or easily predictable mechanical behaviors. Notably, we examined the reaction of the yolk surface of the zebrafish embryo to laser injury employing planar cortex (PC) models with Stokeslet pairs distributed on a 2D square grid. Further, we modeled a spherical fluid droplet sedimenting in an immiscible liquid building spherical cortex (SC) models with evenly distributed Stokeslet pairs around one main axis on a spherical shell with cylindrical symmetry (Fig 3E and Appendix Note S4). In both cases, HR faithfully inferred mechanical pressure fields, power, and surface tension maps from experimental velocity fields (see Appendix Figs S3 and S4 and Video EV3).
To understand the biomechanical behavior of the early zebrafish embryo, we studied the mechanical properties of the yolk and its cortex and employed HR to identify the power sources guiding and/or contributing to epiboly and to evaluate their dynamics.
Yolk mechanical properties.
We assessed the mechanical properties of the yolk by studying its rheological behavior recording the thermal fluctuations of microinjected fluorescent nanoparticles (Daniels et al, 2006). Images of yolk‐injected nanoparticles were obtained, and the position of their centroids, their trajectories, and the two‐dimensional mean square displacement (MSD) of each particle were computed. We found that the MSD exhibited a proportional dependence on time lag consistent with Newtonian fluid behavior and a viscosity of 129 mPa·s (Fig 4A). This viscosity of the yolk is thus on the same order to that recently determined for Drosophila (286 mPa·s; He et al, 2014). Considering this viscosity, the speed of the yolk flows (epiboly progression rates—200 μm/h), the yolk density (similar to water ~1.074 g/cm3; Fujimura et al, 2007), and the size of the embryo (~350 μm of radius), we calculated a Reynolds number for the yolk in the order of 10−7 (as in Drosophila; He et al, 2014), a very low value.
We then performed a direct quantification of the yolk cortex viscoelastic properties. Atomic force microscopy (AFM) was employed to acquire local surface tension measurements on living whole embryos (see Materials and Methods and Appendix Note S5). The viscoelasticity of the cortex was assessed applying low‐amplitude (100 nm) multifrequency oscillations composed of sinusoidal waves of different frequencies (0.35–11.45 Hz). Cortex rheology was dominated by an elastic‐like behavior with the viscous modulus ~5 times lower than the elastic modulus (Fig 4B). Moreover, no significant frequency dependence was found for any of these moduli in the experimental time window. As the ratio between the viscous and elastic moduli is expected to drop with decreasing frequency, our rheological data indicated that at a timescale of minutes the cortex, with an apparent elastic modulus G′ of ~100 Pa (see Equation (6), Materials and Methods), behaves as a very soft solid in the range of the softest biological materials (Butcher et al, 2009; Discher et al, 2009).
Last, from the measured relaxation times determined by AFM, we estimated the effective cortex viscosity at long timescales (fluid‐like response) in the order of ηc ~10 Pa·s. This implies that the ratio between the typical traction in the cortex at long times and the stresses associated with the measured flows of the fluid is just about ηcw/ηR ~0.1, where w stands for the cortex width (1–4 μm; Rawson et al, 2000), R the characteristic embryo size, and η the yolk viscosity. In conclusion, the yolk content can be considered a viscous Newtonian fluid with low Reynolds number and the yolk surface behaves at the timescale of epiboly movements as an extremely soft solid. Importantly, the magnitudes of the cortical and yolk fluid tractions were truly comparable.
Epiboly mechanical power sources.
We applied HR to unveil the spatiotemporal pattern of mechanical power during epiboly. To do so, experimental 2D velocity fields from time‐lapse videos of embryos meridional sections employing yolk granules movements as a reference were determined by PIV. Then, simulated 3D velocity fields were generated from a SC model with Stokeslet pairs distributed on five concentric spherical shells adjusted to the embryo cortex geometry. These theoretical fields were fitted by regression analysis to the experimental fields (Fig 4C and Video EV4).
We found that at the epiboly onset, the active power mainly maps to the blastoderm. Later, a distinct high power domain arises under the cortex between 50 and 60% epiboly at the position and time of the ingression of DCs. Last, once the EVL front crosses the equator, the largest mechanical power concentrates to the actomyosin‐rich E‐YSL, while the EVL adjacent cells display negative power densities opposing deformation. The EVL cells and DCs near the animal pole remain a weak power source (Fig 4D and Video EV5). At all times, the yolk always displays very low mechanical energy dissipation.
Cortical stress patterns.
HR also enabled to analyze the spatio‐temporal progression of local surface tension by mapping animal to vegetal (AV) and circumferential (CC) cortical stresses as a function of time (Fig 4E and Video EV6). We found that at the initiation of epiboly both, AV and CC stresses show an even distribution at the surface of the embryo. As epiboly progresses, however, a vegetalward‐oriented gradient of CC stress is built up throughout the EVL and the yolk, while the AV stress displays a complex distribution. From animal to vegetal, it diminishes with time at the animal pole, dramatically rises at the EVL front, drops at the E‐YSL, and increases again toward the vegetal pole. As a result, at the E‐YSL after 50% epiboly, the difference between AV and CC stresses sustained minimum negative values. In contrast, the stress difference in the adjacent EVL cells was positive and maximum. Summing up, a positive gradient of tension progressively develops at the yolk surface toward the vegetal pole from 50% up to the end of epiboly. This gradient is more pronounced for the CC component.
To verify the predicted axial tensional gradient, we employed targeted laser microsurgery (Colombelli et al, 2009; Appendix Figs S5 and S6 and Appendix Note S6). We found that the experimental disruption in the AV or CC directions of the yolk actomyosin cortex at the E‐YSL resulted in an immediate recoil, with exponentially decaying speed (see also Behrndt et al, 2012). This let measuring the relative magnitude and directionality of local surface tensions (Appendix Fig S7 and Video EV7). Significantly, the differences between AV and CC tension at each time point precisely recapitulated the stress differences estimated by HR (Fig 5A). Thus, the inferred temporal and directional pattern of tension at the E‐YSL was confirmed. Equally, the predicted vegetalward AV tension gradient at the yolk surface was also confirmed at different stages by laser microsurgery at different distances from the EVL/yolk margin (Fig 5B).
Stress‐graded profiles were further verified by AFM. We measured force‐indentation curves at the yolk surface and found that they showed a perfect match to a “liquid balloon” model. This is a viscoelastic mechanical model representing a viscous liquid surrounded by an elastic cortex under mechanical stress (Evans & Yeung, 1989; Lomakina et al, 2004; Krieg et al, 2008). By fitting this model to the data, we calculated absolute values of cortical tension (Tc) for the yolk membrane in the order of 100 pN/μm (Fig 5C). Performing this analysis at two different positions, 50–100 μm ahead of the EVL and at the vegetal pole, detected after 50% epiboly a significant long‐range positive gradient of tension. Remarkably, it also pointed, unequivocally, to the presence of tension at the vegetal pole (Fig 5D). These results are consistent with the short‐range gradient inferred by laser cuts in agreement with the HR analysis.
In short, the predicted yolk surface stress gradient and dynamics were accurately supported by AFM and laser microsurgery.
Mechanically relevant elements leading epiboly movements
Overall, epiboly kinematics suggest that the EVL behaves mesoscopically as a homogeneous deformable sheet, that the E‐YSL shrinks specially after 50% epiboly, and that the rest of the yolk cell membrane is mostly un‐stretchable and kept under increasing tension. Indeed, the observed spatial and temporal differences in surface speed can be consequence of the non‐homogeneous stiffness of the cortex. We evaluated this option analyzing the cortex displacement in response to laser cuts on the yolk membrane. After recoil, the geometrical center of cuts parallel to the EVL consistently shifted toward the vegetal pole (Fig 6A). These data emphasize the difference of stiffness between the EVL (softer) and the yolk (stiffer) surfaces. This difference correlates with previous analysis of micro‐elasticity using Brillouin Scattering Spectroscopy, which also revealed the relative higher stiffness of the yolk versus the blastoderm at early cleavage stages (Fujimura et al, 2007).
Differences in cortex stiffness may be relevant for tissue expansion, and we hypothesize that epiboly, at least after 50%, will be just reliant on the active contraction of the E‐YSL and directed toward the vegetal pole as a reaction to the difference of stiffness between the animal (EVL) and vegetal (yolk) domains, which would establish a cortical tension gradient (Fig 6B). To assess this hypothesis, we simulated epiboly in silico. This simulation took into account the experimental kinematics and assumed both, an exponential decay of the cortical contractility along the φ axis (in accord with the observed exponential decay of speed of the EVL/yolk margin) and that the EVL and the yolk behave as materials of different stiffness (see Appendix Note S7). Simulations precisely recapitulated, at least from 50% epiboly onwards, most of the kinematic and mechanical parameters associated, the yolk internal vortices, and the mechanical power and patterns of surface tension at the cortex (Fig 6C and D, and Videos EV8 and EV9). Simulations also rendered a remarkably accurate profile for the speed of the EVL front displacement along the φ axis, slow first, up to reaching the equator, and much faster afterward.
The in silico mechanical simulations predict that interfering with the contractile capabilities of the E‐YSL will result in a dramatic reduction in the speed of the advancing front and, in practical terms, in epiboly completion failure. Importantly, if we abolished the E‐YSL shrinkage, simulations provided an aberrant regime for the speed of the EVL front with a progressive decline from 50% epiboly onwards (Fig 7A).
We experimentally tested these inferences by genetically interfering in Msn1 expression. Msn1, an ortholog of the Drosophila Ste20‐like kinase Misshapen, is a well‐characterized mediator of E‐YSL contractility. Actin and phospho‐myosin 2 failed to be recruited at the YSL in morphants and the EVL cells at the margin do not stretch longitudinally, while ZO‐1 accumulation (septate junctions) was not affected (Köppen et al, 2006). We addressed the role of msn1 in the E‐YSL injecting morpholinos in the yolk at the sphere stage. This resulted, as previously described, in epiboly's slowness and arrest with 50% penetrance. We performed a kinematic and mechanical analysis of msn1 yolk morphants and found that the failure on E‐YSL contractility was associated with a strong reduction in the speed of progression of the EVL front mostly after 50% epiboly (Fig 7B and Video EV10). We further found by HR that power at the E‐YSL was critically compromised and that stress gradients along the cortex were essentially abolished (Fig 7C and D).
In short, simulations and experimental interferences in E‐YSL contractility show a remarkable correlation strongly supporting that epiboly kinematics and mechanical behavior can be minimally described by (i) the contractile capabilities and gradual change of dimensions of the E‐YSL and (ii) the dissimilar elastic properties of the EVL and the yolk surface (Fig 6B).
We have performed a detailed study to determine the, up to now elusive, major force‐generating elements driving zebrafish epiboly (Solnica‐Krezel & Driever, 1994; Cheng et al, 2004). We revealed several significant evidences: (i) both, the yolk surface (this manuscript and Behrndt et al, 2012) and the EVL tissue (Campinho et al, 2013), display elastic behavior at short timescales. This elastic response is sustained all throughout epiboly; (ii) the yolk cell membrane behaves as an extremely soft solid; (iii) the yolk content is fairly viscous and incompressible and flows at very low Reynolds number with negligible inertia, and (iv) the magnitudes of the cortical and yolk fluid tractions are comparable and of equivalent order. Further, we generated a theoretical approximation to the embryo as a tensile cortex of spherical shape (SC) surrounding a viscous fluid. Applying hydrodynamic regression (HR) to this SC model, we fitted a unique velocity and mechanical pressure field to experimental 2D+T time‐lapses and from here we built a minimal mechanical model of epiboly.
HR is a non‐invasive method able to quantitatively infer dynamic mechanical parameters (mechanical power and surface/cortical tension) from in vivo imaging. HR, in principle, could be applied to images acquired from any microscope including light‐sheet fluorescence (Keller & Stelzer, 2008) or non‐fluorescence nonlinear microscopy (second‐harmonic generation—SHG and third‐harmonic generation—THG) (Supatto et al, 2011) providing that velocity fields could be estimated. It relies on a well‐established hydrodynamic model, which exploits the analytical expressions of elementary solutions of Stokes flow. These elementary solutions (pairs of Stokeslets) are distributed along the cortex geometry and oriented normal to it. Stokeslet pairs model local cortical contractions or expansions and are used to fit simulated velocity fields to the experimental velocity fields estimated by PIV. Due to the force balance between the cortex and adjacent fluids (see Appendix Note S2), the cortical stresses are determined from the fluid velocity fields. Significantly, HR is largely over‐determined, that is, manage to fit a large number of variables with a small number of unknown parameters (the coefficients βe associated with the elementary solutions). A major advantage is that the analytical dynamic pressure field is always unique for a given velocity field. Thus, any change in dynamic pressure distribution leads to different velocity fields, and vice versa. Yet, the model is mesoscopic and potential microscopic inequalities, for example, heterogeneities in the conformation of actin and/or myosin filaments in a network, cannot be directly estimated. However, HR opens up the possibility to indirectly address spatial heterogeneities through the inhomogeneous strengths of the Stokeslet pairs.
HR demands that the fluids surrounding the cortex should behave as Newtonian. However, even if surrounding fluids are non‐Newtonian per se, this does not represent a problem at small strain rates at which non‐Newtonian materials behave as Newtonian. Indeed, the Newtonian behavior of a fluid at low Reynolds number is commonly accepted for processes such as the establishment of the body plan or organogenesis in Caenorhabditis elegans, Drosophila, Xenopus, or mice (reviewed in Cartwright et al, 2009); for example, during gastrulation in the early Drosophila embryo, stresses generated at the surface integrate with the hydrodynamic properties of the interior fluid, which behaves as Newtonian with a low Reynolds number, to transmit force (He et al, 2014).
Although HR can infer power densities and stress tensors regardless of the cortex material properties, it is important to point out that the macroscopic response of many tissues is viscoelastic (where the tissue behaves as an elastic solid over short timescales and a viscous fluid over long timescales). Indeed, zebrafish explants of embryonic ectodermal tissues are viscoelastic (Schotz et al, 2008). Yet, during epiboly in Xenopus (Luu et al, 2011), the ectodermal epithelium mimics an elastic solid when stretched isometrically at a timescale of hours. Considering that developmentally and geometrically amphibian and teleost epibolies are equivalent, and that not only do the relative values of their germ layer surface tensions fall in the same range but even their absolute values are essentially the same (Schotz et al, 2008), we would expect they will show equivalent mechanical properties.
In comparing to HR, current image‐based analytical methods aiming at inferring pressure and tensional parameters seem severely limited. Most of them rely on the application of vertex models to planar epithelia. In these models, the geometry of cell junctions has to be extracted from image segmentation and it is assumed that the cell edge linear tensions and internal pressures are at quasi‐equilibrium. To estimate the internal pressure of each cell and the tension of each cell‐to‐cell edge, two different approaches have so far been undertaken: (i) a statistical model of tension distribution (Ishihara & Sugimura, 2012) and (ii) the solution, in the least square sense, of an over‐determined system of equations imposing a common value for the internal pressure and cortical tensions of each cell (Chiou et al, 2012). None of these methods account for deforming tissues as they model an epithelium at quasi‐equilibrium and hence can only probe static mechanical properties. On the other hand, video force microscopy (VFM) (Brodland et al, 2010) stands apart and allows to derive 2D forces based on tissue deformations estimated from time‐lapse videos using a finite element model. The space is discretized, and partial differential equations are constrained at the boundaries of each partition. The system of linear equations is exactly determined and solved through matrix inversion. This method hence faces the typical error sensitivity inherent to inverse problems. So far, VFM has only been applied in 2D settings and both passive viscous dissipation and long‐range pressure propagation in the third dimension were disregarded. Thus, the pattern of tension simulated cannot model any resisting or active part outside the measurement plane, whereas HR can. We anticipate that as the method stands today, VFM would be unable to properly describe the patterns of flows inside the yolk during epiboly or accurately map cortical power densities and stresses in 3D. Finally, quantitative measurements of cellular stresses within tissues have been recently achieved by shape reconstruction of injected microdroplets using confocal microscopy and image analysis (Campas et al, 2014). Local spatial inhomogeneities in cellular stresses can be measured with this technique, but it fails to measure the isotropic tissue pressure consequence of large‐scale tissue flows. In some sense, this method may complement the information gathered by HR, and combining both methodologies could constitute a sound approach to analyze the overall biomechanics of developing tissues.
Cytoplasmic streaming has been observed in large plant cells like Chara (Yamamoto et al, 2006; Goldstein et al, 2008), in C. elegans (Wolke et al, 2007), and in the mice oocyte (Yi et al, 2011). Velocity fields have been studied and compared to kinematic models, and flows have been linked either to specific molecular motors or to stress anisotropies at the surface. However, in these works, the distributions of power and/or stress have not been inferred. HR provides relative values for mechanical parameters such as stress distributions, and these can be extremely relevant in these models. The ability of accurately inferring mechanical power patterns and surface tension profiles in vivo will enable the exploration of the genetic and cellular mechanisms that account for morphogenesis biomechanics.
A global mechanical analysis of epiboly
In its application to epiboly, HR enabled not only to estimate mechanical power and provide insight in the effective surface tension variations but also to reach a global concerted representation of epiboly in which cell movements can be fully explained by the sole mechanical properties of the system. Before crossing the equator, an active role for the EVL and DCs on epiboly is supported by the mechanical power at the blastoderm (see Fig 4D) and the stresses at the embryo surface (see Fig 4E). EVL and DCs early activity is instrumental for epiboly progression (Trinkaus, 1963; Warga & Kimmel, 1990; Concha & Adams, 1998; Köppen et al, 2006). From 50% onwards, however, we showed through simulations and experimental interference that polarized tensional properties of the yolk surface and the EVL altogether with the contractile capability of the E‐YSL [latitudinal (CC) and longitudinal (AV)] are sufficient to explain both surface and inner movements (see model in Fig 6B). The gradient of tension is validated by our own laser cuts (Figs 5 and 6A), together with those described (Campinho et al, 2013) for the yolk and EVL marginal cells. No differences in the recoiled velocity between parallel (AV) and perpendicular cuts (CC) are observed in the marginal EVL cells at the beginning of epiboly. However, as epiboly progressed, AV tension increased (parallel cuts), while CC tension decreased in the EVL cells. Thus, a difference between AV and CC tension, being AV higher than CC, was generated over time. In the yolk, however, the differences between AV and CC tension were opposite to those of the EVL cells. This model is further backed by the documented preferential orientation of EVL cells under anisotropic tension in the AV direction (Köppen et al, 2006). It is also in agreement with recent observations pointing to the EVL as a continuous elastic medium in which tensions and elastic deformations occur at appropriate length scales (see Campinho et al, 2013). This last issue is consistent with equivalent observations in Fundulus (Weliky & Oster, 1990).
It has been suggested, alternatively, that epiboly may be described as the result of the combined activities of cable‐constriction and flow‐friction motors (Behrndt et al, 2012). In this opposing scenario, tissue deformation would be driven by the ability of the cortical actomyosin ring at the E‐YSL to contract and generate active tension, while an additional pulling longitudinal force in the yolk will be generated by a regime of intermediate friction acting against an actomyosin retrograde flow. This cortical‐friction motor hypothesis derives from a theoretical description of YSL actomyosin network dynamics, and it has not been experimentally tested. It implies that the actomyosin ring, which is mechanically connected to the EVL on its animal side, has a free boundary on its vegetal side. This assumption is essential for the model and demands that the surface tension from the contractile ring toward the vegetal pole should decrease toward zero. We found, however, that laser microsurgery in the E‐YSL unequivocally pointed to an oriented gradient of AV tension from the EVL/yolk margin toward the vegetal pole (Fig 5C). Further, our measurements of surface tension (Tc) by AFM in vivo explicitly revealed a non‐null surface/cortical tension at the vegetal pole, which increased with time (Fig 5E). Therefore, the vegetal boundary of the E‐YSL is not mechanically free, as it was claimed (Behrndt et al, 2012), but linked to the tense yolk membrane. Under these circumstances, the tension originated by the constriction of the actomyosin ring (at the E‐YSL) will be conveyed by the yolk cell membrane toward the vegetal pole. Hence, for epiboly progression, a flow‐friction motor would be, in principle, needless.
In our view, the failure to detect tension at the vegetal pole by laser microsurgery (Behrndt et al, 2012) while we readily detect a significant surface tension with AFM strongly suggests that the yolk surface/cortical tension is indeed sustained at the pole by the elastic membrane itself and not by the underlying cortex. Remarkably, considering the notable fitting of our simulations, any potential anisotropy in material properties of the yolk membrane seems not to be mechanically relevant. In this context, the role of the observed yolk cortical retrograde flows of actin and myosin might be not to provide a friction force, but to replenish actin and myosin, which undergo rapid turnover at the contractile E‐YSL.
In summary, by applying HR, we have characterized the spatio‐temporal evolution of local cortical tensional stresses and power sources during zebrafish epiboly. We also found that the incompressible yolk displayed stereotyped laminar flows and low but measurable mechanical energy dissipation, suggesting a novel morphogenetic mechanism by which internal flows convey a passive transmission of cortical forces leading to tissue deformation. Our simulations and experimental interferences into E‐YSL contractility (msn1 yolk morphants) (Fig 7) further support that the coordinated action of cortical elements (local contractility at the E‐YSL and differential tension between the EVL and the yolk surface) and passive viscous flows in the yolk is instrumental in the mechanical implementation of epiboly movements (Fig 6B). Local surgical detachment of the EVL and the yolk contact edge in Fundulus (Betchaku & Trinkaus, 1978) results in E‐YSL faster progression toward the vegetal pole (as no resistance by the EVL was placed to its forward movement) and EVL cells retraction indicating that the EVL/yolk contact is under tension. This last outcome has also been observed for zebrafish (Behrndt et al, 2012). Thus, the EVL, after 50% epiboly, appears to limit the speed of epiboly progression resisting to deformation. Yet, an active role of the EVL before reaching the equator and a mechanical input of DCs ingression (Trinkaus, 1963; Warga & Kimmel, 1990; Concha & Adams, 1998; Köppen et al, 2006) may be necessary to refine the process.
How the increased stiffness at the yolk cortex is achieved? Why is the yolk cortex stiffer than the EVL? Why does its stiffness increase as epiboly progresses? We could speculate that the EVL cells could easily respond to nearby pulling forces flattening passively, while the yolk, on the contrary, has no way to increase its surface area. As a result, tension would accumulate at the yolk membrane as epiboly progresses, which indirectly would lead to an overall progressive reduction in stiffness of the EVL. Indeed, when we artificially increased contractility in the YSL by microsurgery, we observed that nearby EVL cells deformed more than faraway ones, while the vegetal part of the myosin ring did not deform (see Fig 6A). On the other hand, we have not yet explored in quantitative terms the molecular motors involved or the roles that microtubules or actomyosin flows may have in promoting epiboly movements. This would require developing specific analytical protocols as the hydrodynamic method, as such, does not provide any information about the physical mechanism that generates the stresses. Modeling making use of the measured stresses may be employed as input to define the roles for microtubule bundles or biofilaments and to predict their functional roles.
Local tensional stresses and the passive transmission of forces with elastic and viscous components in an incompressible medium provide a simple physical description for the mechanical control of morphogenetic movements. This biological logic appears to be widespread and it applies, for instance, to the folding of the epithelia in the early Drosophila embryo during gastrulation, where stresses generated at the surface integrate with the hydrodynamic properties of the interior to transmit force (He et al, 2014). We propose that hydrodynamics may be key for the comprehension of the biomechanics of multiple morphogenetic processes.
Materials and Methods
Zebrafish handling and maintenance, zebrafish lines and mRNA injection
Adult fish were maintained under standard conditions. Collected embryos were grown at 28.5°C in E3 embryo medium (Westerfield, 2000) and staged according to morphology as described (Kimmel et al, 1995). AB and TL wild‐type backgrounds were used throughout this study. Membrane‐GFP transgenic Tg(β‐actin:m‐GFP) animals (Cooper et al, 2005) were provided by Lilianna Solnica‐Krezel, and Myosin‐EGFP transgenic Tg(β‐actin:myl 12.1‐e‐GFP) fish (Behrndt et al, 2012) were supplied by Carl‐Philipp Heisenberg. A DNA construct encoding for Lifeact‐GFP (Riedl et al, 2008) and cloned in a zebrafish expression vector was kindly provided by Erez Raz. mRNA was in vitro synthesized (mMessagem Machine kit, Ambion) and injected into the yolk at one‐ or 512‐cell stages (150 pg) as described (Westerfield, 2000).
To create a condition where E‐YSL actomyosin contraction will be blocked, we employed Misshapen 1 yolk morphants (Köppen et al, 2006). The msn1 MO‐splice (5′‐ACACACAACT TACccttaaagtgga‐3′) morpholino that targets the exon 1/intron 1 boundary of the pre‐mRNA (Heasman, 2002) was used (kind gift of Antonio Jacinto). We injected 8 ng of the morpholino in the E‐YSL at the 1,000‐cell stage on Tg(β‐actin:m‐GFP) embryos. As previously described (Köppen et al, 2006), half of the injected embryos showed a delay in the epiboly of the DCs and the EVL. To image the internal yolk granules, we obtained medial sections using two‐photon microscopy as described below. The velocity fields generated with MatPIV were subjected to analysis by HR.
Embryos were dechorionated manually and fixed overnight at 4°C with 2% PFA, 2.5% glutaraldehyde in PB. They were post‐fixed in 1% osmium tetroxide with 0.08% potassium ferrocyanide, gradually dehydrated through an acetone series, infiltrated in EPON overnight, and polymerized for 48 h at 60°C. Ultrathin sections (50 nm thick) were cut on an ultra microtome Ultracut UC6 (Leica microsystems, Vienna). Images were acquired in a JEOL 1010 electron microscope at 80 kV equipped with a Megaview III CCD camera.
Embryos were manually dechorionated and mounted in 0.5% low melting agarose in E3 embryo medium in glass bottom Petri dishes (MatTek, No. 1.5 cover glass) for imaging in both laser scanning and spinning disk confocal microscopy, and in lumox dishes (Sarstedt) for two‐photon microscopy.
To image the internal yolk granules, medial sections (350 μm depth from the yolk cell surface) were collected from Tg(β‐actin:mGFP) embryos using a Leica SP5 two‐photon microscope equipped with a mode‐locked near‐infrared MAITAI Laser (Spectra‐Physics) tuned at 900 nm, with non‐descanned detectors and with a 25×/0.95 water‐dipping objective. Images were scanned at 200 Hz, and frames were averaged three times. Stacks of 30 μm, 10 μm step size, were collected every 2 min.
To visualize the surface of the yolk, Tg(β‐actin:mGFP) embryos were imaged in a Zeiss LSM700 confocal microscope with a 63×/1.4 oil objective. Stacks of 25 μm, step size of 0.2 μm, were acquired.
To study the kinematics of the membrane, dechorionated embryos were incubated in 100 μg/ml lectin‐TRITC (Sigma L1261) for 5 min at the sphere stage, 50 or 85% epiboly, rinsed, and mounted for imaging. Images were acquired in a Zeiss LSM 700 confocal microscope with a 40×/1.3 oil immersion objective. Stacks of 20 μm, 0.39 μm step size, were acquired every 30 s.
A PIV analysis (see below) was applied to time‐lapse videos to estimate the velocity field at each time point.
To simultaneously label the E‐YSL cortical actin and the membrane, Lifeact‐GFP mRNA was injected at 1‐cell stage embryos and they were incubated in lectin‐TRITC as mention above. Images were collected using a Zeiss LSM 700 confocal microscope with a 63×/1.4 oil objective, stack of 11 μm, step size of 0.3 μm.
Most image analyses were performed using Fiji (http://pacific.mpi-cbg.de) and Matlab (Mathworks). To obtain velocity fields through PIV, we applied the MatPIV software package written by Johan Kristian Sveen in Matlab (Supatto et al, 2005).
Laser surgery experiments and retraction analysis
Laser surgery of the actomyosin cortex was performed with a pulsed UV laser (355 nm, 470 ps per pulse) by inducing plasma‐mediated ablation as described before (Colombelli et al, 2009). To compare the cortical tension in the AV and CC directions at the E‐YSL at different stages, a 20‐μm laser line containing 50 pulses was scanned five times at a frequency of 800 Hz, parallel and perpendicular to the EVL front, centering the cut at a distance of about 20 μm, through a 63×/1.2 W objective lens. Fluorescence imaging was performed through a custom spinning Nipkow disk unit equipped with a 488‐nm laser line and a Hamamatsu ORCA CCD camera, acquiring at 1.5 frames per second.
The same type of laser ablation regime (same laser line length, same pulsed UV laser settings) was applied to validate the AV surface gradient of tension revealed by the regression analysis. This time, only cuts parallel to the EVL were performed at 0, 10, 20, 40, and 60 μm of the EVL border, and at different epiboly stages (50, 60–70, and 80–90% epiboly).
In all cases, to estimate the tensions, we followed the accepted assumption that the tension of the actomyosin cortex before the laser cut is proportional to the outward velocity of the immediate recoil (Grill, 2011). Retraction analysis was performed as follows: First, from the microscopic video time‐lapse recordings of the laser cuts, a region of 160 × 160 pixels was cropped and rotated so the cuts were centered and positioned horizontally. 160 pixels correspond to the apparent length of the original cuts. We then applied a PIV analysis (see Image Analysis) to measure displacements from frame to frame. We applied interrogation windows of 64 square pixels with three‐fourth overlap generating a 7 × 7 grid. From these grids, we averaged the vertical components of the displacement vectors of the three upper and lower rows and calculated the retraction speed V from their difference.
Quantification and statistics
Yolk granules kinematics.
From the velocity fields measured by PIV in the medial plane of three sibling embryos, we estimated the speed as an average over 10 frames, every 10 frames. We also computed the vorticity of the coarse‐grained velocity field and used its maximum absolute value to position the vortex centroid.
Kinematics of the membrane.
The velocity component along the AV axis obtained from lectin‐label embryos was averaged from parallel rows and consecutive time points at different stages. To compare the different videos, velocities were normalized to the velocity of the margin of the EVL.
To compare the cortical tension in the AV and CC directions at the E‐YSL at different stages, we calculated the mean recoil velocities for each perpendicular direction in each stage and tested the statistical significance of their inequality using the Wilcoxon rank sum test. To quantify the AV surface gradient of tension, we normalized the recoil velocities at each distance from the EVL margin to the recoil velocity at 20 μm distance and tested the statistical significance of the gradient using the Wilcoxon rank sum test.
Cortical tension measurements with atomic force microscopy
Zebrafish epiboly embryos were manually dechorionated and mounted for AFM as follow: Small holes of approximately twice the embryo size in diameter and half the embryo size in depth were made with thin forceps in an embryo medium 2% agarose layer prepared in 35‐mm Petri dishes. Embryos were placed in the holes made in the agarose layer and tightly held by surrounding them with 0.5% low melting agarose. They were rotated to probe the region of interest by AFM (e.g., E‐YSL or the vegetal pole).
The cortical tension of the embryos was measured with a custom‐built atomic force microscope (AFM) attached to an inverted optical microscope (TE2000, Nikon, Tokyo, Japan) using a previously described method (Alcaraz et al, 2003). The embryos were probed with a spherical polystyrene bead 4.5 μm in diameter attached to a Si3N4 cantilever with a nominal spring constant of 0.01 N/m (Novascan, Ames, IA). Force‐indentation (F‐h) curves were recorded in two different regions of the yolk surface. Each region was probed in five different points located at the center and the corners of a 10 × 10 μm square.
The vertical displacement of the AFM cantilever (z) was measured with strain gauge sensors coupled to the piezo‐actuators. The cantilever deflection (d) was recorded by a quadrant photodiode using the optical lever method. The slope of a d‐z curve obtained from a bare region of the coverslip was used to calibrate cantilever deflection. The cantilever spring constant (k) was determined by the thermal fluctuations method. The force (F) on the cantilever was computed as F = k × d. In each measurement point, five force‐displacement (F‐z) curves were acquired by ramping the cantilever at 1 Hz with a peak‐to‐peak amplitude of 5 μm (velocity = 10 μm/s) up to a maximum indentation of ~2 μm. The indentation of the sample (h) was computed as h = (z – zc) – (d – doff), where zc is the position of the contact point and doff is the offset of the photodiode.
Force curves acquired on embryos were interpreted in terms of a “liquid balloon” model consisting of an elastic cortex with cortical tension Tc enclosing a viscous liquid. Assuming a small indentation compared to the size of the embryo, force increases proportionally with indentation (Lomakina et al, 2004; Krieg et al, 2008) as (3)where Rb is the radius of the bead and Re is the radius of the embryo. Since the radius of the embryos as estimated by microscopy images (~350 μm) is two orders of magnitude larger than that of the bead (2.25 μm), this equation can be approximated as (4)Tc was computed for each F‐h curve by nonlinear least‐squares fitting (Matlab, MathWorks). For statistical analysis, the cortical tension in each embryo region was taken as the average value of Tc computed from the F‐h curves recorded at the five different measurement points.
The “liquid balloon” model has been previously described (Ishihara & Sugimura, 2012) and has been employed to estimate cortical tension of leukocytes by micropipette aspiration. This model has been extended to spherical cells indented with a microbead by means of micropipettes (Lomakina et al, 2004) and more recently applied to compute the cortical tension of spherical progenitor cells from gastrulating zebrafish embryos indented with spherical AFM tips 5 μm in diameter (Krieg et al, 2008). Here we adapted this AFM technique to measure Tc of zebrafish whole embryos. Although the “liquid balloon” model equation assumes small sample deformations, as previously reported (Krieg et al, 2008), we found a proportional F‐h relationship up to indentations close to the radius of the microbead. Taking into account the linear behavior of the force curves and the sensitivity of parameter estimation, we fitted the “liquid balloon” model to the whole indentation range of our measurements to obtain a robust estimation of Tc.
The force‐indentation curves recorded in each of the two regions exhibited a proportional relationship (Fig 5C), supporting the liquid‐droplet model. The model fitted all the data very well (R2 > 0.999). To evaluate the dependence of Tc on the velocity of the cantilever, two additional sets of five F‐h curves were acquired at 3 and 30 μm/s in the last measurement point. Although Tc increased by 13% (P < 0.01, t‐test) when the cantilever velocity was reduced from 10 to 3 μm/s, no significant change was observed from 10 to 30 μm/s. This weak dependence of Tc on cantilever velocity provides further support to the interpretation of this parameter computed by AFM as the cortical tension in the embryo.
Rheology of the cortex with atomic force microscopy
The viscoelastic properties of the cortex were measured in five embryos with AFM by applying a low‐amplitude (100 nm) multifrequency oscillation composed of sinusoidal waves of different frequencies (Alcaraz et al, 2003). An effective complex modulus g*(f) was computed in the frequency domain as (5)where i is the imaginary unit and F(f) and h(f) are the frequency (f) spectra of force and indentation, respectively. The term “b” is the correction for the viscous drag on the cantilever computed from oscillations applied above the surface. g*(f) was separated into real and imaginary parts as g*(f) = g′(f) + ig″(f). g′(f) is the elastic modulus and is a measure of the elastic energy stored and recovered per cycle of oscillation. g″(f) is the viscous modulus that accounts for the dissipated energy. The loss tangent was also computed as g″(f)/g′(f), which is an index of the solid‐like (< 1) or liquid‐like (>1) behavior of the material.
To compare the stiffness of the embryo with that of other biological materials, we also computed the complex shear modulus (G*(f) = G′(f) + iG″(f)) assuming the contact model of an equivalent homogeneous half‐space. Transforming equation 43 (see Appendix) to the frequency domain, G* can be computed as (6)where ho is the operating indentation (~0.5 μm).
Rheology of the yolk measured by particle tracking
The viscoelastic behavior of the yolk was assessed by recording thermal fluctuations of microinjected fluorescent nanoparticles (radius a = 100 nm, FluoSpheres F8811, Life Technologies). Two hours after microinjection, the embryos were placed on the stage of a confocal inverted microscope and images of nanoparticles were recorded for 26 s at a sampling rate of 25 Hz. The position of the centroid of the particles and the trajectories of the particles were determined with ImageJ. The two‐dimensional mean square displacement (MSD) of each particle was computed with custom‐made software (Daniels et al, 2006; Wirtz, 2009) as (7)where t is the elapsed time and Δt, the time lag. In a pure viscous liquid, MSD increases inversely with viscosity (ν) according to the Stokes–Einstein relationship (8)where T is the absolute temperature. The viscosity of the yolk was calculated by fitting this equation to the MSD data.
AH‐V and MM performed all biological tests; P‐AP designed the modeling and the regression analysis; JC, AH‐V, MM, and P‐AP performed the laser cuts and data refinement; TL and DN provided the AFM and particle tracking rheology data; ST and IP overview the hydrodynamic analysis and code implementation, and EM‐B designed the study, analyzed the data, and wrote the manuscript. All authors discussed the results and commented on the manuscript.
Conflict of interest
The authors declare that they have no conflict of interest.
In first place, all authors like to acknowledge the essential conceptual and experimental contribution of Philippe‐Alexandre Pouille, who passed away during the process of revision and publication of this manuscript. We would like to dedicate this work to his memory. We thank Jacques Prost for a thorough analysis of the HR methodology from a theoretical point of view. We also thank the Confocal Microscopy Unit from IBMB‐PCB, the Advanced Digital Microscopy Core Facility from IRB Barcelona, Xavier Esteban and members of the EMB laboratory for continuous support. We are grateful to Javier Buceta, Carolina Minguillón, Emmanuel Farge, Damian Brunner, Katerina Karkali, and Carla Prat for reading earlier versions of this manuscript and John P. Trinkaus for inspiring work. The Consolidated Groups Program of the Generalitat de Catalunya, and DGI and Consolider Grants from the Ministry of Economy and Competitivity of Spain to EMB, IP and DN supported this work.
FundingGeneralitat de Catalunyahttp://dx.doi.org/10.13039/501100002809
- © 2016 The Authors