Strain maps characterize the symmetry of convergence and extension patterns during zebrafish gastrulation

During gastrulation of the zebrafish embryo, the cap of blastoderm cells organizes into the axial body plan of the embryo with left–right symmetry and head–tail, dorsal–ventral polarities. Our labs have been interested in the mechanics of early development and have investigated whether these large-scale cell movements can be described as tissue-level mechanical strain by a tectonics-based approach. The first step is to image the positions of all nuclei from mid-epiboly to early segmentation by digital sheet light microscopy, organize the surface of the embryo into multi-cell spherical domains, construct velocity fields from the movements of these domains and extract strain rate maps from the change in density of the domains. During gastrulation, tensile/expansive and compressive strains in the axial and equatorial directions are detected as anterior and posterior expansion along the anterior–posterior axis and medial–lateral compression across the dorsal–ventral axis and corresponds to the well characterized morphological movements of convergence and extension. Following gastrulation strain is represented by localized medial expansion at the onset of segmentation and anterior expansion at the onset of neurulation. In addition to linear strain, symmetric patterns of rotation/curl are first detected in the animal hemispheres at mid-epiboly and then the vegetal hemispheres by the end of gastrulation. In embryos treated with C59, a Wnt inhibitor that inhibits head and tail extension, the axial extension and vegetal curl are absent. By analysing the temporal sequence of large-scale movements, deformations across the embryo can be attributed to a combination of epiboly and dorsal convergence-extension.


Scientific Reports
| (2021) 11:19357 | https://doi.org/10.1038/s41598-021-98233-z www.nature.com/scientificreports/ When these forces are coordinated across cells, tissue-level strains and tension are generated which contribute to the dynamics of embryonic tissue folding and dorsal closure as well as wound healing and other tissue-scale movements 14 . Thus, morphogenesis is the product of cell-to-cell signaling at a local scale regulated by gradients of morphogens at the whole embryo scale. While biochemistry, cell biology, and cell-level mechanics have explained key features of morphogenesis, the role of stress and strain in morphogenesis is now being explored. Recent advances in imaging and image processing approaches can now extract from whole embryo time-lapse images several key signatures of large-scale dynamics such as velocity fields, cell density changes, and 2D tissue deformations 15,16 . However, to detect and describe the mechanics of tissue-derived forces in the context of an entire embryo requires large-scale 3D in toto imaging with: (1) high speed to capture developmental dynamics, (2) high resolution to resolve individual cells, and (3) high penetration depth to image whole organisms. While confocal and multi-photon microscopy have contributed fundamental information on early development, light sheet configured microscopes have emerged as the instrument of choice in studies of Zebrafish or Drosophila embryogenesis 16,17 . With the temporal and spatial resolution and reduced photo-bleaching of light sheet microscopy, it is possible to track every cell during early development 16 . As a result, the technical imaging challenges are replaced by the computational challenge of extracting the tissue-level mechanics from very large image datasets.
To detect and characterize the large-scale tissue dynamics during gastrulation of the zebrafish embryo, we have adapted the 2D tectonics-based analysis of tissue-level strain to the 3D analysis of deformation over the whole embryo. The resulting strain rate maps reveal the time-course and location of areas undergoing compaction or expansion, i.e., linear strain, during well-known morphogenetic movements including convergence and extension, neural plate formation, and somite formation. In addition to axial and equatorial strains, new features such as rotational strain or curl appear during convergence and extension as well as the left-right laterally symmetric presence of strain across the embryonic axis. When embryos are treated with the Wnt inhibitor, C59, the resulting strain maps show that extension but not convergence is abolished.

Results
To detect and characterize the tissue dynamics during gastrulation, first we imaged all cells continuously from 80% epiboly (8 hpf) to 3 somite stage (11 hpf) (Fig. 1A, Movie-1), counted the cell population (supplementary Figure 2), calculated their velocity fields, and color-coded the direction of the velocity fields (Fig. 1). The velocity fields cover the entire surface of the embryo and capture the major morphogenetic movements in three principal directions: (1) epiboly as posterior translocation (red), (2) convergence as lateral translocation (blue), and (3) extension as anterior translocation (green). These movements persist through developmental stages 80% epiboly (8 hpf), 100% epiboly (10 hpf) and 1-somite (10.5 hpf) (Fig. 1B,C). In the ventral view of the 10 hpf bud stage embryo, the extending tail bud region is seen extending anteriorly from the vegetal pole.
From the velocity fields we then calculated 3D tissue deformation or strain by determining the 3D strain rate tensor in three principal directions: the anterior-posterior (AP), media-lateral (ML), and radial (r) directions (supplementary Figures 3-7). In our analysis, we display only the AP (axial) and ML (equatorial) as these directions align with the major changes during morphogenesis (Fig. 2B). Several patterns of strain are present in the maps of the dorsal and ventral hemispheres of the embryo (Fig. 2C, Movie-2). First from 8 to 10 hpf, the dorsal hemisphere is undergoing M-L compaction while simultaneous A-P expansion is restricted to an equatorial band around the dorsal but not ventral half of the embryo. Second, the dorsal M-L compaction is accompanied by extensive M-L expansion in the ventral hemisphere. Third, by 10:30 hpf the dorsal compaction narrows to the A-P body axis but is now flanked by a prominent a pair of M-L expansions (asterisk). This expansion is transient, lasting until 11 hpf and diminishes afterwards (n = 8) (Fig. 2C, Movie-2). At this time the ventral hemisphere compacts along the A-P axis as the body axis extends into the ventral side to form future head and tail of the embryo. We observed variabilities in terms of the magnitude of deformation ( Supplementary  Figures 3 and 4). However, major events such as the side expansion patterns during these stages are preserved over different samples.
To represent the time course in the evolution of strain, we constructed kymographs of the three strain components in two important regions of the body over different developmental stages (  (Fig. 3C). The M-L and A-P strains are bilaterally symmetric across the 0° dorsal point. Finally, radial strain ( Fig. 3D) is initially detected as compaction in areas midway between the dorsal and ventral sides of the embryo but by 10:30 hpf the dorsal point of the embryo shows prominent expansion (Fig. 3D).
We then analyzed strain along the body axis of the embryo (Fig. 4). Figure 4A indicates the continuous compaction along M-L direction at dorsal site. We observe that the dorsal compaction becomes stronger, extending from the equator plane towards 2 polar regions. Again, the kymograph shows this M-L compaction comes along with A-P expansion during epiboly (Fig. 4B), followed (10.5 + hpf) by an increase in radial strain at the same spatial position (Fig. 4C). This site also coincides with the formation of early somites as well as the neural tube.
A deformation or change of the domain, i.e., strain trace, can be calculated from the divergence of the strain rate tensor over a time-period (set at 30 min). The resulting strain trace map shows a large-scale change in the ventral hemisphere that extends medially and dorsally prior to 100% epiboly, followed by change along the dorsal body axis. These tissue deformations are inhomogeneous across the embryo and during morphogenesis (Fig. 5, Movie-3, Supplementary Figures 5 and 6).
Having identified the patterns of strain along the main embryo axes, we then visualized rotations in the velocity field (magnitude of the anti-symmetric part of the velocity gradient tensor) (Fig. 6, Movie-3  To test the sensitivity of our approach in the detection of abnormal developmental activities, embryos were exposed to the dorsalizing Wnt inhibitor C59. Because the epiboly was delayed by C59 treatment (Supplementary Movie 4), the embryos were staged on the progression through epiboly instead of time after fertilization (Fig. 7A,E). Compared with the untreated embryos (Fig. 7B), the M-L strain in the C59 treated embryos (Fig. 7F) showed similar pattern of dorsal strain with centralized compaction along dorsal midline and at later stages the lateral flanking regions becoming less compacted or even expansive). Similarly, in the C59 treated group (Fig. 7F), M-L compaction in the dorsal hemisphere and lateral expansion regions were prominent with the timings comparable to the control group.
However, strain in the A-P direction significantly differs between treated and control embryos (Fig. 7C). Most prominent is an equatorial band of compression at 50% epiboly that disappears by 80% epiboly and is replaced by a weak area of expansion centered in the vegetal quadrant by 90% epiboly. (Fig. 7G, supplementary Figure 7). Correlated with the effects of C59 in reducing A-P strain is an effect on rotational strain (Fig. 7D). In the control group, the four hemispheres of the embryo develop counter-rotating regions of curl by 60% epiboly with a region of zero curl at the dorsal and ventral poles. In contrast, the location, area, and magnitude of the rotational quadrants in the C59 treated group (Fig. 7H) showed clear differences with the control. By 60% epiboly, the two animal hemispheres displayed weak curl and were located away from the dorsal pole of the embryo. In comparison, the two quadrants of curl in the vegetal hemisphere developed at the margin of the blastdoderm by 50% epiboly but did not extend to cover the entire vegetal half of the treated embryo. Because the quadrants

Discussion
In this study we have detected and quantified three types of tissue deformation during epiboly and gastrulation: anterior-posterior (AP) and the medio-lateral (ML) strain, radial strain, and curl. Figure 2 presents a schematic representation of the deformation patterns. Notable features of these deformations are their long-range coupling and symmetry. An example of long-range coupling is first seen during epiboly (7-10 hpf) when the cells migrate in all directions to cover the yolk. This expansion is visible from the anterior view of the embryo as A-P strain (Fig. 2). However, before the cells reach the vegetal pole the strain maps also detect signatures of strain corresponding to convergence as ML compaction over the entire dorsal hemisphere. This compaction is accompanied by a broad equatorial band of expansion in the AP direction centered at the midpoint of the dorsal axis. The dorsal midpoint corresponds to the boundary between head and trunk structures. Interestingly, these lateral displacements of tissues on the dorsal side are accompanied on the ventral side by an opposite strain pattern (Fig. 4). Cells converging along the dorsal side must be balanced by cells expanding somewhere else, as observed on the ventral side during and after epiboly. Also, tissues tend to expand in AP on the dorsal side  www.nature.com/scientificreports/ (although this effect is diminished after the end of epiboly) and converge on the ventral side. Thus, the major morphogenetic movements of epiboly and convergence/extension appear to be mechanically coupled over the dorsal and ventral surfaces of the embryo. In addition to the expansion and compression of these movements is the bilateral left/right symmetry that is evident in the ML strain kymographs (Fig. 3B,C) during epiboly and continuing past gastrulation. This symmetry is defined by the animal/vegetal and dorsal/ventral axis from 8 to 10.5 hpf. Perhaps the most striking example of the symmetry in the strain maps is seen at the first appearance of a somite (10.5 hpf). There is a dorsal band of ML compaction which is flanked by a short-lived (10.5-11 hpf) and narrow region of ML expansion (asterisk). This "hot spot" region spatially coincides with the paraxial mesoderm, also known as presomitic or somitic mesoderm, which commits to the formation of somites 18 . However, this local strain precedes the time when morphogenesis corresponding to somitogenesis occurs and suggests that paraxial mesoderm expansion is a possible biomechanical signal for somitogenesis before apparent morphology formation. The correspondence between events on opposite side of the embryo illustrates that strain is bilaterally symmetrical and equally generated on the left and right sides of the embryo.
In contrast to linear strain that is expressed along the surface of the embryo, we are able to detect radial strain during morphogenesis (Figs. 3D, 4D). The most prominent deformations arise at 11 hpf when the neural structures are forming along the dorsal midline and in the ventral hemisphere. The ventral regions are undergoing significant expansion especially in areas corresponding to the head and tail structures. Thus, the radial strain corresponds to well characterized events in neurulation.
A third and particularly interesting type of deformation that we have characterized is rotational strain or curl. Curl appears in the dorsal hemisphere during epiboly and then in the ventral hemisphere when epiboly continues to the ventral pole. There are two striking aspects of the location and pattern of curl. First, the direction of curl is consistent with a large-scale integration of convergent movements toward the embryonic axis and extension along the embryonic axis in anterior and posterior directions into a unified rotational movement. Second, the direction of curl on the left and right sides of the embryo are opposite and complementary. When drawn on the surface of a sphere, these patterns of deformation clearly define the quadrants associated with the rotational patterns ( Figs. 6 and 7). These patterns reflect the bi-lateral symmetry of the cellular flows and coordination of morphogenetic movements over very large areas of the embryo.
Further, a curious feature of convergence and extension that emerges from our analysis of the velocity fields is a divergent point or saddle at the dorsal midpoint from where cells translocate either anteriorly or posteriorly. In the curl maps, this divergence point corresponds to the region where the domains of positive and negative curl intersect. This point of divergence may have a functional significance. It corresponds to the site of ingression during epiboly and the boundary between head and trunk structures in the embryo. Our maps do not reveal the mechanism that directs cells to move or expand in opposite directions. We are investigating possible structural and molecular explanations for divergent cell migration.
Notably, the sensitivity of our approach in detecting developmental abnormality was demonstrated by the drug perturbing experiment. To ensure reproducibility across imaging platforms, we measured the effects of C59 with the diSPIM which has an advantage of easier sample mounting and improved embryo survival rate compared with the DSLM. One significant drawback of this approach is that embryos could be imaged from only 2 orthogonal views in diSPIM instead of 4 in DSLM.
Wnt signaling is known as a ventralization factor during zebrafish development, the inhibition upon which causes dorsalized phenotypes [3][4][5][6]8 . Unlike a complete absence of convergence-extension (C&E) when blocking BMP gradient 1 , Wnt inhibition only reduced the magnitude and extent of C&E 6,19,20 , and therefore is a good candidate for testing the sensitivity of our approach. Previous studies 6,19,20 reporting the perturbation of Wnt signaling attributed the defects to an incomplete convergence and extension (C&E). It however remains unclear whether such perturbation affects both collective cell behaviors and whether other types of movement like epiboly also got affected. Here the strain maps not only detected the effects of Wnt inhibitor C59 on gastrulation but also revealed that extension but not convergence was predominantly affected. This weakening of extension easily explains the dorsalized phenotype.
At 50% epiboly, the A-P axis compacts at the blastomere margin from a combination of ingression/shielding and convergence/extension and then is soon replaced by a prominent A-P expansion in control group as the embryo axis extends. However, in the C59-treated embryos compaction remains possibly from a stalling of leading edge cells and the absence of extension from loss of AP patterning information from Wnt gradient [3][4][5]21,22 . Furthermore, according to the previous study, reduced friction forces between prechordal plate cells and neuroectoderms were observed in the silberblick morphants 23 . The reduced friction force is closely related to the slowdown of relative movement between these two populations, thereby contributing to a reduction in the extension 23 . Since C59 blocks both the canonical and non-canonical Wnt signaling pathways, the accumulated effects from the above would result in a predominant loss of the extension, and ultimately affect the linear and rotational strain map patterns. Without extension, cells collect at the dorsal hemisphere and further broadens the dorsal axis, a feature of dorsalization at later stages.
In conclusion, the strain maps and 3D tissue tectonics are the first step to quantify the global biomechanics of developing embryos and provide a sensitive way to detect dynamic spatial-temporal abnormalities much earlier than previous approach. To completely understand the force-extension pattern within the tissue, the lack of validated methods to measure tissue stiffness in real-time remains a challenge. However, it has been reported the cell density is correlated with local tissue stiffness 23 , and thus offers us an opportunity to gain an insight of relative tensions among tissues. A resulting tissue deformation-model can act as a mechanistic model to understand the coordinated activities of morphogen gradients, signaling pathways 24 and help identify local biophysical events like friction force mediated cell movements between different cell populations that drive embryo morphogenesis 25  Embryo microinjection. Wild-type Zebrafish were raised and obtained from Prof Christoph Winkler's lab and it is a part of NUS-DBS fish facility. The embryos were collected in the aquarium at 1-cell stage and immediately injected with 9.7 nL of H2B-GFP mRNA (150 ng/µL). Each embryo was kept at 28.5 °C in the incubator. At 8 h post-fertilization (hpf), the embryos were dechorionated, embedded in 0.4% low melting point agarose within a cleaned FEP tube, and then positioned in the DSLM for live imaging. The bottom of the FEP tube was sealed with 2% agarose solution. In diSPIM configuration, during sample mounting, the embryos were mounted into 0.8% low gelling agarose in a culture dish, which was filled with E3 medium or Wnt inhibitor C59 solution subsequently.  Figure 1). 3D datasets were collected with two independent light sheets for excitation and single detectors in the Multiview mode, at 2 min intervals from 7.5 to 12 hpf.
In diSPIM system, before transferring to the microscope, the sample were added with either Wnt inhibitor C59 pre-dissolved in DMSO at the final concentration of 1.5 µM, or the equal amount of DMSO for control group. Following this, the samples were transferred to the 3i Marianas Lightsheet microscope system (Intelligent Imaging Innovations, Inc) with the same environment settings as above. The imaging duration was around 12 h starting from 5 hpf, with 2 min interval. In this system, 2 orthogonal views were imaged, and the region covered by each view is extended beyond the midplane of the front hemisphere (Supplementary Figure 1). Image processing. Each experiment generates around 2 TB stack of TIFF image files which are processed with a MATLAB library (MathWorks, Inc) 10 to identify and track all nuclei in the embryo. The workflow includes: (1) 3D image reconstruction from 2D stacks, (2) image segmentation and cell identity registration using iterative thresh-holding method, (3) post-segmentation spatial filtering on intensity distribution, intensity value, anisotropy and connectivity, (4) temporal filtering that recovers false negatives and removes false negatives with reference to the nearby time points, and (5) 3D-cell tracking program adopting the nearest-neighbor algorithm.

Strain maps.
To calculate 3D strain, we modified a previously described tectonics-based method for measuring strain in 2D 9 . The general approach is to group all cells into overlapping spherical domains, calculate the velocity field gradients of the domains, and then extract the strain from the tensors of the velocity field gradient matrix. Because the early embryo is a sphere and to simplify the calculations, we start by mapping the embryo onto a unit sphere where the radius r 0 is set to 1. Then we discretise the spherical embryo surface into 2000 equally spaced points or "nodes" (Fig. 2A). Each node, − → r i , lies at the vertex of a triangular mesh and its location is described in polar coordinates ⇀ r (r, ϕ) where r is radial distance from the embryo centroid, is polar angle or elevation, and ϕ is azimuth. For each time point, t, we define local domains around each node by selecting all the cells within a certain distance from the node center. In this study, the radius of these domains is set to cover 110 microns which is small enough to capture relevant spatial patterns, while providing enough statistics to calculate the local strain rates. At 10 hpf, the domains would contain 10-250 cells. Once the domain is fixed, we obtain its local velocity field, v, from the mean of the cell velocities within a domain.
The velocity field is then analysed to extract the local strain rates, or the rate of tissue deformations. To calculate the strain rate, we differentiate the velocity field spatially, forming a 3 × 3 velocity gradient tensor L( ⇀ r , t): This tensor is then separated into its symmetric and anti-symmetric components. The symmetric part quantifies local deformations as the strain rate tensor: The tensor's eigenvalue components are the magnitude of strain rates along the corresponding eigenvectors (r,θ,ϕ ), and the sum of the diagonal components or trace of this matrix describes the total volume change of the domain. Because the major changes in development occur along the embryonic axis (A/P, Dorsal/Ventral) and equator (left/right), we decomposed all the vectors into anterior-posterior (A-P), medio-lateral (M-L), and radial (R) directions. The accumulated deformation or strain rates over time give rise to the magnitude of strains, in www.nature.com/scientificreports/ our case, we chose the time window of 30 min depending on the noise level to filter random noise and highlight significant signals. In addition, in the perturbation experiments the time window was fixed at 40 min for both control and drug treated groups, since the significant biomechanical signature only started to appear for such time intervals in drug treated group. The strains are colour coded where expansion is red, indicating diverging velocity fields or an overall increase in inter-nuclei distances while blue represents converging velocity fields or the cell cluster in the domain compacts and gets denser (Fig. 2B). In addition to linear strain, the anti-symmetric part of the tensor is the spin matrix and contains information about the local rotation or curl of the domain (Fig. 6): Positive or counter-clockwise rotation is red, negative or clockwise rotation is blue.
Kymographs. The time evolution of strain is captured in a moving 30 to 40 min window from an annular region with a span of 0.38 radians. The axial, equatorial, or radial strain is averaged over domains along the equatorial circumference which bisects the dorsal (0°) and ventral (± 180°) poles of the embryo or axial circumference in the AP/DV plane.