Quantitative analysis of tissue deformation dynamics reveals three characteristic growth modes and globally aligned anisotropic tissue deformation during chick limb development

Tissue-level characterization of deformation dynamics is crucial for understanding organ morphogenetic mechanisms, especially the interhierarchical links among molecular activities, cellular behaviors and tissue/organ morphogenetic processes. Limb development is a well-studied topic in vertebrate organogenesis. Nevertheless, there is still little understanding of tissue-level deformation relative to molecular and cellular dynamics. This is mainly because live recording of detailed cell behaviors in whole tissues is technically difficult. To overcome this limitation, by applying a recently developed Bayesian approach, we here constructed tissue deformation maps for chick limb development with high precision, based on snapshot lineage tracing using dye injection. The precision of the constructed maps was validated with a clear statistical criterion. From the geometrical analysis of the map, we identified three characteristic tissue growth modes in the limb and showed that they are consistent with local growth factor activity and cell cycle length. In particular, we report that SHH signaling activity changes dynamically with developmental stage and strongly correlates with the dynamic shift in the tissue growth mode. We also found anisotropic tissue deformation along the proximal-distal axis. Morphogenetic simulation and experimental studies suggested that this directional tissue elongation, and not local growth, has the greatest impact on limb shaping. This result was supported by the novel finding that anisotropic tissue elongation along the proximal-distal axis occurs independently of cell proliferation. Our study marks a pivotal point for multi-scale system understanding in vertebrate development.

: Implantation experiments using beads soaked with inhibitors of FGF (SU5402) and SHH (cyclopamine). Most beads were implanted in the posterior and/or distal sides of the right limb bud. In some embryos, beads were also implanted in the middle of the right limb bud. The left limb bud was used for the normalization of limb bud size to remove embryo-to-embryo variability. We performed the implantation experiment for two consecutive time intervals . (A) As in the estimation of the deformation map, a few tens of DiI and DiO fluorescent markers were injected into the mesenchyme of each limb bud. The arrow heads indicate implanted beads. (B) Each limb bud was divided into multiple triangles based on the marker positions. From the relative positional changes of the markers, we calculated the area growth rate and deformation anisotropy vector for each triangle in each limb bud. We then compared the spatial patterns of growth rate and deformation anisotropy with those Development | Supplementary Material of the wild type limb buds, which were calculated by using the distribution of fluorescent markers before deformation and the 2D estimated map 2D . (C) To quantify the effects of inhibiting signaling on deformation anisotropy, "difference in deformation anisotropy vectors at X for the i -th bud and for the wild type. and are eigenvalues and eigenvectors of deformation gradient tensor (see Materials and Methods). The reason we use this index for quantifying the effect on the anisotropy is that this value tends to be larger in the case where the differences in the orientations of both vectors and the magnitudes of the vectors are larger. In the region where the anisotropy is small, its orientation is sensitive to various randomness, such as embryo-to-embryo variability of tissue deformation dynamics and measurement noise, which should not be included in the evaluation of the effects of inhibiting SHH or FGF signal on deformation anisotropy. (D) Distributions of "change in tissue growth rate" and "difference in deformation anisotropy" calculated for different triangles for different zones (distal or posterior) and different time intervals . The graphs in Fig. 6B were made based on these distributions. The change in tissue growth rate is defined as the ratio of the growth rate with beads to that in the control case.

S1. Mathematical details for Fig. 1F
The 1D tissue growth dynamics is given by the following equation: is the position at time t of the cell that is initially located at X . ) , ( t X s is the tissue growth rate that the cell X experiences at time t . In this equation, the coordinate of the leftmost (or the proximal end) is set to be 0. In Fig. 1F, for simplicity of explanation, we assumed spatially and temporally uniform growth, i.e., constant ) , , so that the growth dynamics is simply given by . Next, let us consider the camera position. When the camera is fixed at the proximal end, since the position of each cell X relative to that of the camera (denoted by ) , , the velocity of the cell is observed as: Therefore, in this case, the velocity becomes higher for cells located farther to the right or more distally. It should be noted that this does not mean that the motility of cells located at the distal side is higher because uniform-growth is assumed. This example shows that tissue deformation results in apparent movement of cells. On the other hand, when the camera is fixed at the distal end (i.e., the rightmost), the relative position is given as: where 0 L is the initial length of the limb bud. Thus the velocity is: In this case, the velocity becomes higher for cells located farther to the left or more proximally. Again, this does not mean that the motility of cells located at the proximal side is higher. It should be noted that the limb bud does not grow inside the body, although the velocity of cells on the proximal side looks faster. In the case where the camera is fixed so as to focus on the cell located at the center of the tissue, the following holds: Development 142: doi:10.1242/dev.109728: Supplementary Material In this case, the velocity becomes maximal at the two ends of the tissue.
Finally, we note that the above argument holds not only for uniform growth, but also for arbitrary non-uniform growth dynamics.
S2. Data for quantifying tissue deformation dynamic (details) (i) DiI and DiO markers on the frontal plane (spanned by the P-D and A-P axes). We injected the fluorescent marker DiI/DiO at approximately 20-50 sites in limb mesenchyme at a depth of about 100 m from the dorsal ectoderm using a sharp tungsten needle that had been dipped in a 1 mg/ml DiI/DiO solution (Fig. 2B). We used only limb buds of the same length at each stage by measuring with an ocular micrometer. Fluorescent images were taken using a Leica M205FA microscope with a digital zoom focus to adjust the precise magnification every time. The diameter of each marker was about 50 m . To estimate the 2D deformation map for each time interval, we used marker positional data from 7 or 8 embryos (representing a few hundred markers in total). Since limb bud size was slightly different between embryos, data for each embryo were normalized along the PD and AP axes.
(ii) Measurement of dorso-ventral thickness with an OPT scanner The thicknesses of the limb buds along the dorso-ventral axes were measured with an OPT (optical projection tomography) scanner 3001 (Sharpe et al., 2002). The shape of each limb bud was captured by detecting autofluorescence of the limb bud embedded in 1% low melting point agarose diluted in PBS. To avoid shrinkage and deformation of the tissue due to dehydration, samples were not cleared with Benzyl benzoate/Benzyl alcohol. Data on DV thickness were extracted by Avizo software.

S3. Validation of the estimated deformation map for chick hindlimb development
We evaluated the accuracy of the estimated deformation maps by their predictive performances (cross-validation) (Figs. 2F, G). For each time interval, 90% of the positional data for injected fluorescent markers before and after deformation was used for estimating the deformation maps by the proposed Bayesian method, and the remaining 10% of data was used for calculating the predictability that was evaluated by the mean square error of predicted positions of markers after deformation. The error was calculated along the proximo-distal and antero-posterior axes separately as follows: Error along P-D axis: is the positional data of the i -th marker. p N is the number of data points used for calculating the predictability.

S4. Correlation analysis between area and DV growth rates
In Fig. 3B, for the different time intervals, we calculated correlation coefficients between area and DV growth rates, AD , as follows: where ) ( i A x and ) ( i D x are the area and DV growth rates at each lattice point i x , respectively (see also Fig. S2). A E and D E are their averages over the whole limb bud. A and D are their standard deviations. The correlation between area and DV growth is low ( 4 . 0 1 . 0 AD ) although there are a few stages in which the spatio-temporal patterns look similar.

S5. Measurement of cell cycle time (details)
Cell cycle time was measured by pulse-chase labeling with BrdU/IddU (Martynogaa et al., 2005). First, 100 l of IddU (100 mg/ml) was injected into the yolk sac, and embryos were incubated at 38.6 o C for 1.5 h. Next, 100 l of BrdU (100 mg/ml) was injected into the yolk sac, and embryos were incubated for another 30 min, followed by fixation with 4% paraformaldehyde. Immunohistochemistry of BrdU and IddU was previously described (Boehm et al., 2010). When the difference in the labeling time between IddU and BrdU is 1.

S6. Measurement of SHH signaling activity using a luminescent system (details)
The Gli-responsive element (Sasaki, H., et al., Development 124, 1313-1322, 1997 After normalization of emerald luciferase activity to red Phrixothrix luciferase activity, we visualized and quantified SHH signaling activity in the limb bud using MetaMorph software.

S7. In silico experiment (a modified vertex dynamics model)
To determine if spatially-biased growth or deformation anisotropy is the greater contributor to limb-specific morphology and limb elongation, we performed in silico experiments for two different situations.
(I) The situation in which the spatial pattern of tissue growth rate is the same as that in normal development but the anisotropy is low over the whole limb bud. (II) The situation in which the spatial pattern of anisotropy is similar to that in normal development but the tissue growth rate is uniform. In order to obtain more significant differences in morphology, the simulation time interval was set at 24 hours:  (Fig. 8). To simulate morphologies for both situations (I) and (II), we developed a modified vertex dynamics model ( Supplementary Fig. S4A). In our model, each piece of tissue is represented by a tetragon formed by linking four vertices. Each vertex moves so as to decrease energy of the whole system, called the Hamiltonian (denoted by H ): where i x is the spatial coordinate of vertex i ( N i , , 1  ) and t is the time. Depending on the Hamiltonian adopted, the dynamics of vertices change, leading to a different final morphology of the vertices network.
In our simulation, we used the Hamiltonian given by: where the first and second term on the right hand side are the constraints of volume growth rate and tensor components at vertex i , respectively. ) (i F and ) (i F jk are the matrix representations of the deformation gradient tensor F at i (for a given coordinate system) and its ) , ( k j -component, respectively. ) (i F is defined as the average of four linear transformations for the deformations of four adjacent squares ( Supplementary Fig. S4B). ) ( det target i F is the target volume growth rate. In simulation (I), it is set to be the same value as the wild type at each location (i.e., )

S8. Calculation of deformation anisotropy in the experiment of inhibiting cell proliferation with Trichostatin
We first calculated the deformation gradient tensor F from the relative positional changes of four DiI markers; concretely, the tensor was obtained by minimizing the mean square error, where ) , ( i i Y X and ) , ( i i y x are the position vectors for each marker i before and after deformation, respectively, when the origins of coordinates were set to the centroids of tetragons before and after deformation, respectively. From the resulting F , the deformation anisotropy was calculated as the ratio of eigenvalues of the stretch tensor U . )) ( ( )) ( ( 4