Modeling and live imaging of mechanical instabilities in the zebrafish aorta during hematopoiesis

All blood cells originate from hematopoietic stem/progenitor cells (HSPCs). HSPCs are formed from endothelial cells (ECs) of the dorsal aorta (DA), via endothelial-to-hematopoietic transition (EHT). The zebrafish is a primary model organism to study the process in vivo. While the role of mechanical stress in controlling gene expression promoting cell differentiation is actively investigated, mechanisms driving shape changes of the DA and individual ECs remain poorly understood. We address this problem by developing a new DA micromechanical model and applying it to experimental data on zebrafish morphogenesis. The model considers the DA as an isotropic tubular membrane subjected to hydrostatic blood pressure and axial stress. The DA evolution is described as a movement in the dimensionless controlling parameters space: normalized hydrostatic pressure and axial stress. We argue that HSPC production is accompanied by two mechanical instabilities arising in the system due to the plane stress in the DA walls and show how a complex interplay between mechanical forces in the system drives the emerging morphological changes.


Experimental background
Zebrafish is a primary source of knowledge in HSPC development studies due to its advantages in fast reproduction cycle and high transparency in the optical spectrum 50 . Despite the tremendous popularity of this model organism, data with high temporal and spatial resolution on its dorsal aorta morphogenesis has become available only in the second decade of the twenty-first century 34 . Below we briefly review the known data on the development of the zebrafish dorsal aorta during the hematopoiesis (when the endothelial-to-hematopoietic transition takes place).
During the time period relevant for EHT (~ 24-72 h post-fertilization) 12,13,34 , zebrafish DA is a long tube formed by a monolayer of flat endothelial cells (so-called squamous epithelium) (Fig. 1). The embryonic dorsal aorta is filled with fluid and immersed in the surrounding tissue matrix. While these characteristics of the DA remain unchanged, it undergoes a series of drastic morphological changes 11 .
At 24 hpf, the embryonic heart starts beating and initiates blood circulation. From that moment on, DA rapidly increases its diameter from ~ 24 to ~ 32 µm at 42.5 hpf. Then at around 40 hpf, the aorta develops a periodic pattern of alternating thinner and thicker regions (see Fig. 1). Its appearance coincides with the initial extrusion events of future HSPCs. After peaking at 42.5 hpf, the DA diameter starts to decrease. At 65 hpf, the aorta recovers its original cylindrical shape and diameter 11 .
Interestingly, the HSPC extrusion rate peaks between 42.5 and 52 hpf precisely when the aorta diameter starts to decrease 11 . High-resolution fluorescent microscopy also allows following fates of individual cells and reveals details of the EHT at cellular level. An endothelial cell destined to become a hematopoietic one migrates toward the ventral part of the dorsal aorta (it is still unknown when precisely cell fate differentiation takes place). Once there, the cell drastically changes its morphology. Starting as a flat elongated cell with virtually equal apical and basal membranes areas (characteristic diameter ~ 10 µm) and submicron thickness, it undergoes a strong anteroposterior contraction and acquires a round plate-like shape, buckling outside of the dorsal aorta lumen toward the cardinal vein 11,13,34 . The process is finalized by the formation and closure of an actin ring around the emerging cell. The actin ring or actomyosin ring is the most archetypical structure of cytokinesis, which is an essential part of cell division. The ring contains actin filaments cross-linked by the motor protein Myosin II 51 . Neighboring endothelial cells also assist the process by accommodating their shape and ensuring aorta integrity to prevent blood hemorrhage.

Results and discussion
Continuous model of the dorsal aorta. Models based on the concept of energy minimization have found much success in describing self-organization processes, including morphogenesis 1,21,52,53 . Here, we develop a micromechanical model of the dorsal aorta. We define the elastic energy of the system and, importantly, we show how the regulation of mechanical properties and stress generation assists the molecular mechanism driving EHT. As shown in Fig. 1, morphological changes of the DA arising during the EHT process are reminiscent of deformations observed in biological and non-biological systems with similar geometry. For example, analogous deformations occur in ordinary thin cylindrical shells under axial load 17,19 , swelling gels 54 , and tubular lipid membranes 25 . Despite the similarities between tubular lipid membranes (TLMs), whose stability was extensively studied in our previous works 25,27,55 , the DA has a number of features distinguishing it from lipid tubes.
A major difference is a nonzero shear modulus of the DA. Unlike lipid membranes with virtually non-existent shear modulus, zebrafish DA, composed of endothelial cells, resists all types of deformations (shear, bending, and stretching/compression) because of the rigidity of cytoskeleton of individual cells and cell to cell adhesion 7 . In mathematical terms, this means that the deformation energy stored in the system depends not only on the final shape of the tube but also on all the displacement field components (unlike TLM deformation energy, which is a function of radial displacements only). Thus, we model the DA as a two-dimensional solid cylindrical shell (such description of monolayer epithelium as a two-dimensional sheet is quite common 8,9,52,56 ) and parametrize its surface as: where R is the tube radius in relaxed state and u = [u r , u ϕ , u z ] is a displacement field characterizing the DA surface in cylindrical coordinates. Since the DA has a nonzero shear modulus, the vector u has both radial and tangential components.
Another essential feature of the DA is that it finds itself in a stressed state for most of the developmental process (like many other biological tissues in physiological conditions 57 ). The stress should have both inhomogeneous and homogeneous components. The inhomogeneous contribution may be associated with the difference in the EC expansion/contraction rates and blood flow, whereas homogeneous one may originate from hydrostatic blood pressure and different growth rates of the aorta and surrounding tissues. We start with a simplified symmetric model accounting for only homogeneous stress and then discuss its possible modifications. Basing on the observed shape of the DA and its displacement relative to the surrounding tissues, we have already hypothesized that axial load induced by the tissue growth rate mismatch plays a crucial role in aorta morphogenesis 11 . Different growth rates of two adjacent tissues can lead to the situation when the faster growing tissue is under (1)  www.nature.com/scientificreports/ confinement, and the growth rate difference causes effective compression. Confined growth of tissues is actively studied and attributed to the formation of the shapes of different organs 20,22 . Here we assume that the dorsal aorta is subjected to the homogeneous stress σ zz along its axis together with hydrostatic blood pressure and write the energy of the DA as follows: where ε ik is a two-dimensional nonlinear deformation tensor (please, see Supplementary materials Eq. S4), and µ are two-dimensional analogs of Lame coefficients, κ is a bending rigidity, �H ′ = (∂ 2 ϕ u r + R 2 ∂ 2 z u r − ∂ ϕ u ϕ )/R 2 is a renormalized average curvature variation 58 . The integration is performed over the DA surface S. V stands for the DA volume and P stands for the pressure difference in the system (for more details see Supplementary materials Eq. S6).
Homogeneous mechanical stress induced by the external tissues and hydrostatic pressure (terms σ zz ∂ z u z and V P ) leads to a homogeneous deformation of the cylindrical shell, maintaining its symmetry. A total displacement field (caused by both thermal fluctuations and external forces) of the DA surface can be represented as: where the homogeneous part is given as u 0 = ε 0 ϕϕ R, 0, ε 0 zz z . By substituting ε 0 ϕϕ and ε 0 zz in Eq. (2) and minimizing the equation with respect to these quantities, one can obtain that the field u 0 corresponds to the following linear deformations of the membrane: We assume that both control parameters PR and σ zz are significantly smaller (in absolute value) than any of the Lame coefficients. Thus, we ensure that deformations (4) are small enough to be considered as perturbations of the initial state. Now we can find the free energy of the system u ′ , as a function of the additional virtual displacement u ′ , taking stressed cylinder with energy �(u 0 ) as a reference state: After performing all the substitutions in Eq. (5), we use the relations σ zz , �PR ≪ , µ and linearize the model. Thus, we retain only terms of the first order in σ zz and P and terms of the first and second orders in the components of the displacement field u ′ and its spatial derivatives: where ε ′ ik is the linearized deformation tensor 58 . Before proceeding to the stability analysis, let us justify the assumed smallness of displacement fields u and u ′ . Biological tissues under physiological conditions often exhibit nonlinear elastic properties, meaning that material constants should be redefined at every step of the deformational process. Such approach is often used to describe living tissues 57 and is known as mechanics of incremental deformation 59 . In addition to elastic deformation, the morphogenesis of the DA involves genetically preprogrammed processes leading to the variation of the size and mechanical properties of the cells. We use an adiabatic approximation and assume that the system evolves through the sequence of equilibrium states 3,21 . This approach allows us to use our linearized model for the stability analysis at various steps of the process, even though DA dimensions change substantially during the considered period. Critical deformations, driving the instabilities can be considered as small perturbations with respect to the stretched but still highly symmetrical cylindrical state.
Providing the DA has the length L we expand the displacement field associated with these perturbations as: where k m = 2πm L , n and m are wave numbers and A j n,m , (j = 1..3 ) are complex amplitudes of the field harmonics. After substitution of Eq. (7) in Eq. (6) and integration the linear term in Eq. (6) disappears to obtain: The matrix M defines the shape of the DA. In the high symmetry phase, the matrix M should satisfy the inequality Det(M) > 0 , which means that all modes of the displacement field u ′ vanish at equilibrium, and DA preserves the shape of a straight circular cylinder. Points of the parameter space �PR, σ zz , which violate this inequality for specific n and m wave numbers, correspond to the adjacent phases with critical modes of the displacement field lowering DA symmetry. When building the model, we have not made any assumptions concerning symmetry of these phases; therefore, analyzing M is sufficient for finding all instabilities changing the cylindrical shape of the aorta. In the following sections, we examine critical modes of the DA and determine the stability region of the system.

Critical modes with rotational symmetry.
After the matrix M has been obtained, numerical stability analysis is relatively straightforward, providing material constants characterizing the system are known. Unfortunately, measuring their values in vivo is extremely challenging in zebrafish DA. Currently, we have no experimental data on material constants at any point in time relevant to EHT. Moreover, they are also likely to change throughout the whole developmental process. Therefore, an analytical study of the DA stability is currently the most rational option. We perform it below by introducing several assumptions that do not affect the generality of the result obtained.
We start with the case when the system preserves its rotational symmetry and assume n = 0, which corresponds to so-called corrugation modes. Thus the matrix elements M 12  Since DA radius is much smaller than its length (R << L), we can consider the wave vector k m as a variable with a continuous spectrum of values. Solving the Eq. (10) with respect to σ zz and minimizing the result with respect to k m (by its absolute value), we obtain the following expression for the critical wave vector: As we shall demonstrate later, the critical radial corrugation mode with this wave vector is responsible for the destabilization of the DA. Using the relations σ zz , �PR ≪ , µ , one can obtain a corresponding critical value of σ zz : +2µ is a two-dimensional Young modulus and ν 2D = +2µ is a two-dimensional Poisson's ratio. Since the thickness of the DA wall is small compared to the DA diameter, we can assume that κ ≪ E 2D R 2 , meaning the expression in the brackets is approximate to 1 . As Eq. 12 shows, negative longitudinal stress σ zz , arising due to the difference in the growth rates of the aorta and surrounding tissues, may cause a mechanical instability, that leads to the formation of the corrugation pattern. On the other hand, the stretching stress term − PR 2 , corresponding to blood pressure dilating the tube, tends to stabilize the system. It is also worth noting that, according to Eq. 12, two factors may contribute to the eventual loss of the DA stability. First, in Eq. 12, the increase of the aorta radius decreases the value of the second term associated with the rigidity of the system. Moreover, the dilatation may also weaken the hydrostatic pressure, thus, effectively decreasing the forces which stabilize the system. Simultaneously, during the development of the embryo, the growth rate mismatch leads to an increase in the negative longitudinal stress, which ultimately exceeds the stabilizing forces and causes the instability in the DA. Later in the text, we discuss in greater detail the development of the instability in vivo.
Equation (12) is reminiscent of the classical expression for the stability of the axially loaded cylindrical shell: where h stand for the wall thickness, E is Young modulus, and D is the cylindrical rigidity of the shell 17 . www.nature.com/scientificreports/ Instabilities violating rotational symmetry. In this section, we examine the stability of the system with respect to modes violating rotational symmetry. Euler buckling is a well-known instability arising in axially loaded tubes and rods leading to the bending of their main axis 18 ; it has also been found in adult arteries 52,60,61 . The author of the work 61 studied the buckling of an artery subjected to hydrostatic blood pressure and axial tension. According to their results, the artery may buckle even under stretching tension σ zz > 0 , providing a sufficient pressure delating the tube �P > 0 is applied. The wave vector of this instability is reported to be k m = π/L . Buckling instability with the same wave vector develops in axially loaded rigid rod with pinned ends 18 . As proven by experimental observations, this type of instability is not provoked in the DA during hematopoiesis because of the surrounding tissues 11,34 , like, for example, the rigid notochord, which is in close contact with the DA apical part. Nevertheless, since our model can be applied to other tubular structures, it is interesting to compare the stabilities of the buckling mode and the corrugation mode studied in the previous section. To do that, we solve the equation Det(M) = 0 , taking the wave number n = 1. In this case, the determinant of the matrix M, Eq. (9) is an eighth-order polynomial with respect to the wave vector k m , which has the following form: where k m = 2πm L and coefficients a i are functions of the control parameters P and σ zz . The determinant (14) does not contain any free terms and thus vanishes at k m = 0 . This trivial solution represents the zero-energy Goldstone mode that corresponds to the tube translation as a whole and has no influence on the stability of the system. To calculate the critical value of stress σ zz leading to the instability that breaks the rotational symmetry of the system, we introduce a normalized wave vector k = k m R . Following the work 52 , we use the long tube approximation (L >> R) and assume that the normalized wave vector k has a small yet finite value. In this case, the critical σ zz value can be estimated from the equation a 1 k 2 + a 2 k 4 = 0 . Taking into account the smallness of the wave vector ( k << 1 ) and using the assumption �PR ≪ µ, , we obtain the following expression for the critical load: Now let us consider another type of instability invoked by the pressure difference P . This time we assume that after destabilization, the tube preserves its translational symmetry and thus substitute k m = 0 in the matrix M . If |n| > 1 , the following inequality governs the stability of the system: Modes with |n| ≤ 1 should be considered separately since translation and rotation of the tube as a whole break the inequality (16) but do not break the tube stability. Solving Eq. (16) with respect to P and taking into account the relations κ ≪ R 2 , µR 2 , we obtain the expression for the critical value of the pressure difference: From the inequality (17), one can easily obtain that the instability is associated with a two-fold degenerate mode |n| = 2, k m = 0 and occurs when the pressure difference reaches the critical value: According to this classic result 17 , the rigid cylindrical tube cross-section subjected to a sufficient hydrostatic compression takes an oval shape. Analogously to the previously discussed Euler instability, however, this transversal buckling is irrelevant for the DA, since the pressure difference caused by the blood flow is obviously positive. Observed ovality of the DA cross-section is due to an anisotropic compression applied by the muscles located on the sides of the aorta and not due to the transverse buckling with spontaneous symmetry breaking caused by the blood pressure 11 .
Since EHT process is associated with deformation and, importantly, buckling of future HSCPs, it is interesting to shortly discuss the stability of individual cell within the theoretical framework we introduced. Obviously, the stability region is different from that of the whole tube. To estimate it, we model the cell as an approximately flat square plate with sides of length a = R. Then the critical stress value is given by the following classic equation: where k = l c π/R and the coefficient l c depends on the type of boundary conditions. Based on the experimentally observed shape of the extruding cells and behavior of their borders during EHT 11,13,34 , we believe that pinned boundary conditions are the most suitable ones; for the cell with pinned borders l c = 1 . Equation (19) can be obtained from the equation M 11 = 0 after performing the limiting transition from cylindrical to planar geometry, i.e., R → ∞ . Equation 19 shows that buckling instability can be initiated by the compressive stress or changes in the rigidity κ of the extruding cell. The stress may be associated with passive processes such as inhomogeneous growth or active ones such as actomyosin contraction. This result will be used afterward.
Stability domain of the cylindric phase. Let us compare the tube stability with respect to the corrugation mode, the transversal buckling mode, and the long-wave buckling mode of the system; to do that we introduce small dimensionless parameters: the normalized pressure p = �PR/E 2D , the normalized axial stress www.nature.com/scientificreports/ σ = σ zz /E 2D , and the normalized bending rigidity γ = κ/E 2D R 2 . For the sake of clarity, we consider the tube made of an incompressible material, taking ν 2D ≈ 1 (unlike in 3D, where incompressible material corresponds to ν = 1/2 ). Then by keeping only the lowest non-vanishing terms in γ and κ we rewrite Eqs. (12), (15), and (18) as: with σ corr , σ buck being normalized axial stress values leading to the corrugation and buckling instabilities in the DA, accordingly, and p trans being normalized pressure leading to the transverse buckling of the DA cross-section. Equation (21) shows that the critical value of the axial stress σ zz leading to the Euler buckling instability is proportional to the R/L ratio squared, whereas the stability of the corrugation mode is independent of the relative length of the tube (see Eq. (20)). Figure 2 shows the stability domain of the system in the parameter space p, σ for the tube with radius to length ratio R/L = 1/50 , which agrees with the observations in zebrafish embryo. In the considered approximation, critical values of the control parameters corresponding to the Euler instability lie along the straight horizontal line (depicted in red). The corrugation mode loses stability along the straight line with a slope ≈ − 1 2 , which crosses the y-axis at the point σ = −γ . And finally, vertical straight line p = −3γ 2 corresponds to the transversal buckling instability associated with the mode |n|= 2, m = 0.
Contribution of surrounding tissues to the stability of the system. In the previous sections, we have analyzed how the material constants characterizing the tube influence the stability domain of the high-symmetry cylindrical phase. The influence of the surrounding tissues has only been considered indirectly by restricting the symmetry of the low-symmetry phase. In this section, we show that the stability of the DA strongly depends on the properties of the surrounding medium. The simplest way to consider an interaction between the tube and the outer tissues is to introduce the term s C 2 u r 2 dS to the energy (6), where C is a phenomenological radial pinning coefficient 60,62 . Since in general for tubular shells immersed into a rigid matrix, the tangential pinning is much weaker than the radial one 60,62,63 , for the sake of simplicity, in the energy (6), we keep only the term www.nature.com/scientificreports/ proportional to the square of radial displacement. Let us note that according to the experimental data between 30 and 60 hpf, individual cells migrate along the DA surface, which also justifies why we discard the tangential pinning in the energy (6). This modification of the energy leads to the emergence of the single additional term C in the element M 11 of the matrix (9). Consequently, the inequality (17) describing the stability of the system with respect to transversal buckling modes is modified to: If the pinning effect is relatively weak CR 4 ≤ κ , the transversal buckling instability is associated with the mode |n|= 2, m = 0, like in the free tube. However, when the interaction with the medium is significant CR 4 /κ ≫ 1 , the tube loses stability with respect to modes with higher wave numbers |n|> 2. The expression for the critical σ zz value (12) corresponding to the tube corrugation is also modified and takes the following form: where E ′ 2D = E 2D + CR 2 is an effective Young modulus. Thus, the previous estimation of the critical stress value leading to the corrugation remains relevant, the only difference being an effective increase of the Young modulus of the tube, caused by the introduction of the pinning contribution into the model. On the contrary, the buckling modes with ( n = ±1, k = ±2π/L ) become non-critical because of an additional constant term that appears in the polynomial (14).
An analogous result was obtained by the authors of the work 52 where a similar system was studied. The authors showed that in an artery surrounded by a rigid tissue matrix, instability is associated with corrugation modes (in their paper, this type of instability is referred to as "varicose instability") and not with classical Euler buckling modes, like in the work 61 . Despite some similar results regarding the influence of the environment on the possible critical modes of the tube, our work has several significant differences from the paper 52 . First of all, in Ref. 52 , artery is not subjected to a real axial longitudinal mechanical stress. Instead, a homogeneous isotropic tension is introduced; this tension is analogous to the liquid-liquid interface energy in the theory of lipid membranes. The authors of the work 52 justify the use of such, unconventional for solids, type of energy by the complex structure of arteries. They model arteries as elastic tubes covered by a monolayer of epithelial cells from inside. According to this reference, epithelial cell division and growth in confinement caused by surrounding tissue results in an analog of an isotropic interface energy. While conceptually this idea is analogous to the compression by growth rate mismatch, we believe that mathematical description, developed in Ref. 52 , is inapplicable for the DA description. As the experimental data shows [11][12][13]34 , the ECs deform inhomogeneously, which is crucial for the EHT. Thus, a concept of conventional anisotropic stress typical of solids should be used.
Bridging theory and experiment. Now, let us analyze the development of the zebrafish dorsal aorta between 25 and 65 hpf-period relevant for the definitive wave of hematopoiesis-within the model framework, and validate our model by comparing its predictions to new and published experimental data. While focusing on the changes in the mechanical properties of the system accompanying hematopoiesis, we also briefly discuss microenvironmental cues leading to these changes. We consider DA morphogenesis as a movement of the system in the parameter space, p, σ , of the normalized pressure difference and the normalized longitudinal stress respectively. As it was previously mentioned, we assume that the system passes through a succession of equilibrium states.
In order to reinforce experimental data on the values of material constants characterizing the system, we perform a zeroth-order estimate. This parameter provides a better context for how the model reacts to changes in these values in physiologically relevant ranges and to give an idea of how material constants can be deduced for future experimental studies. The parameter γ = κ/E 2D R 2 can be estimated within the Foppl von Karman theory for thin membranes 18 . According to this theory, the bending rigidity κ is expressed as E 2D h 2 12(1−ν 2D 2 ) , where h is the membrane thickness 62 . We choose the DA radius to be R ≈ 14µm , which corresponds to the period when it loses stability with respect to corrugation mode (see, Fig. 3). Assuming wall thickness is h ≈ 500nm , and 2D Poisson's ratio is ν 2D ≈ 0.9 , we obtain γ ≈ 1/50 . Choosing such an estimate of 2D Poisson's ratio, we follow the approximation that epithelial cells are nearly incompressible 64 . And finally, the pinning coefficient that describes interaction with surrounding tissues only renormalizes the Young modulus when dealing with corrugation mode, making the system more robust to bending (see Eq. 20). In the phase diagram (Figs. 2 & 3), even a two-or three-fold increase of the estimated normalized bending rigidity only slightly changes the slope of the blue line, which separates corrugated and cylindrical phases, and lowers the line down. The phase diagram does not change qualitatively; the system can switch between phases while staying in the same quadrant of the parameter space, meaning that positive hydrostatic blood pressure and compressive longitudinal stress can drive the process.
The trajectory describing the system development can be subdivided into the following steps (Fig. 3): 1. At around 24 hpf right before the start of a heartbeat, the pressure difference in the system should remain negligible, and the dorsal aorta has a shape of a cylinder with a diameter of 24 ± 0.8 µm 11 . With the initiation of the blood flow and the rise of the hydrostatic blood pressure, the DA diameter starts to increase with the www.nature.com/scientificreports/ expansion of the EC apical/basal membranes. The increase of the aorta diameter and pressure difference contributes to the increase of the normalized pressure difference p . While the difference in the growth rates of aorta and surrounding tissues may cause axial compression, it is still insufficient to destabilize the system and cause symmetry breaking deformations. 2. At around 40 hpf, the aorta develops the periodic pattern with sinusoidal modulation of DA diameter. This means that the axial stress, caused by the interaction with the surrounding tissues, eventually overpowers the stabilizing effect of the hydrostatic pressure and provokes corrugation instability with the vector k m (depicted by the black arrow in Fig. 3). 3. The DA average diameter peaks at approximately 42 hpf, reaching 32 ± 0.9 µm. From that moment on, the EC expansion can no longer compensate for DA surface area decrease due to the extrusion of the future HSPCs. This process gradually relieves the longitudinal stress in the DA walls, which corresponds to the decrease of σ . The aorta diameter reduction corresponds to the decrease of the parameter p , even at constant blood pressure. 4. Eventually, the decrease of σ brings the system back to the high symmetry phase. At approximately 65 hpf, the aorta restores its initial shape and diameter (those at 24 hpf) and HSPCs production stops.
Comparing experimental data to modeling results, one can notice that observed deformations are inhomogeneous with respect to the angular coordinate (see Fig. 3a). Multiple studies report that corrugation has the maximum amplitude in the dorsal part of the aorta and is absent in the ventral one, which is only natural considering that DA is confined by the rigid notochord at the top and by the muscle tissue at the sides. However, when performing stability analysis, we assume that elastic constants characterizing the aorta and its surroundings do not depend on the azimuthal angle. This assumption is crucial for finding an analytical solution to the problem. We indirectly account for the inhomogeneities by selecting instabilities relevant for the zebrafish DA and extrapolate the results on the experimental data. While inhomogeneous rigidity of surrounding tissues indeed plays an important role in the process, allowing cells to extrude from the DA in the vicinity of the cardinal vein so they can enter the blood flow, qualitatively results of our analysis still hold, and the predicted succession of shape phase transitions agrees with experimental observations.
To better understand how mechanical forces orchestrate the HSPCs production, let us compare DA shape evolution in wild-type (WT) and mutant embryos. In Refs. 12,45 , it has been shown that dorsal aorta of silent-heart (SH) mutants with no blood flow and mutants with 50%-reduced blood flow are severely flattened and have decreased cross-section area with respect to WT embryos. Reduction of the average DA radius has also been observed in mutant embryos with disrupted mechanosensing 12 . Based on these results and ability of the ECs www.nature.com/scientificreports/ to sense shear stress 48,65,66 , we speculate that in WT embryos, blood flow-induced shear stress activates genetic mechanisms leading to the expansion of the individual cells, thus overall aorta radius increase. Importantly, aorta of the SH mutants is thinner than that of mutants with 50%-reduced blood flow and with disrupted mechanosensing. This allows us to assume that, while the shear stress mainly plays a signaling role, the hydrostatic blood pressure directly stretches the tube and helps it preserve a more circular cross-section. Translational symmetry of the DA also points toward the secondary role of the blood flow-induced shear forces in deforming the aorta directly. If the DA shape had strongly depended on the shear stress value, then the DA radius, initiation of the corrugation instability, and cell extrusion should have been highly inhomogeneous along the DA length, which is not observed in WT embryos. Corrugation instability signaling that the longitudinal stress has reached a critical value appears in all the previously considered cases, regardless of the blood flow, which is also in accordance with our model; the lack of hydrostatic blood pressure means lowered stability of the DA to the longitudinal stress originating from the growth rate mismatch (see Eq. 12). In the developed theoretical framework, the development of mutant embryos is described by trajectories with a steeper slope in the p, σ space since the decrease of both the blood pressure and the DA radius leads to the decrease of the normalized pressure p . The steeper the trajectory, the faster it reaches the corrugation phase (see Fig. S2), which is in accordance with experimental data 12,45 . However, it is important to note that mutants rarely survive the whole period of 24-65 hpf, meaning that in reality, corresponding trajectories would be incomplete. Now let us focus on the cellular level and discuss the roles of the blood flow and mechanosensing in the EHT. It has recently been shown 12,45 that both the blood flow inhibition and abnormal mechanosensing lead to premature initiation of the extrusion process. Moreover, the process itself is disrupted in the mutants. Cells are extruded not only outward but also inward the aorta. Upon the extrusion, they have an abnormal ovoid shape instead of a cup-like one. Many cells burst when leaving the aorta. Because of premature extrusion, ECs are not subjected long enough to the molecular and mechanical environmental cues, including shear stress, necessary to specify and maintain HSPC identity 12,44,67 , which hinders the HSPC production.
Disruption of actomyosin machinery normal functioning affects cell extrusion similarly to the absence of blood flow 11,12,45 ; the cells buckle abnormally (inside and outside of aorta), do not undergo anterior-posterior constriction like in WT embryos, many of them are extruded only partially and burst. Thus, even though some sort of cell buckling occurs in mutant embryos, actomyosin ring formation around the cells, which are ready to leave the aorta, and its subsequent contraction are essential for the successful extrusion.
Based on the abovementioned results, one can readily make two assumptions. First, gene expression controlling assembly and contraction of the actomyosin ring is probably activated by the shear stress. That is why blocking the blood flow itself, suppressing the ability of ECs to sense it, and disrupting normal actomyosin machinery, all have very similar effects on the EHT process.
Second, buckling of the individual cells toward the sub-aortic space preceding the cell extrusion is another mechanical instability driving the HSPC production. The instability is activated by the plane stress in the aorta walls and not by the blood pressure. The authors of Ref. 13 expressed an analogous idea by stating that both the deployment of pushing forces by the adjoining endothelial cells and the actin ring contraction drive the EHT process. Figure 3a shows that wave vector k m of the DA corrugation instability (Eq. (11)) is larger than characteristic cell size and, thus, the wave vector k of the buckling instability (Eq. (19)) in an individual cell. One can conclude that buckling mode in most ECs comprising the aorta is more stable than the critical corrugation mode in the DA. Nevertheless, some cells, particularly ones undergoing EHT, do buckle even in embryos with abnormal actomyosin machinery. Thus, it is natural to assume that the differentiation starts prior to the cell extrusion and changes mechanical properties of future HSPCs, making them less stable to buckling than other cells in DA walls. However, when and how ECs 'decide' to become HSPCs is yet to be discovered.
While there are some discrepancies in experimental data, results of several studies suggest 11,12,34 a correlation between formation of corrugation pattern and initiation/peak in EHT events. Considering this fact, the formation of the corrugation pattern may be viewed as an indication of the stress reaching a critical value that is sufficient for the extrusions of future HSPCs. However, a connection between these two events requires further experimental studies to be confirmed. Such an interpretation agrees with the blood flow reduction leading to premature cell extrusion. As we have shown, blood pressure stabilizes the system serving as a counteracting force to the compression induced by the growth rate mismatch between the aorta and the surrounding tissues. Thus, in the absence of blood flow, longitudinal stress in the aorta walls reaches the critical value earlier, leading to both the formation of the aorta corrugation pattern and the buckling of the individual cells.

Conclusion
We have studied the evolution and instabilities of the zebrafish dorsal aorta shape during the definitive wave of hematopoiesis. A new micromechanical model of the DA has been developed. Based on our analytical results, new 4D confocal microscopy data, and previously published experimental results, we argue that conversion of ECs, comprising the aorta, into HSPCs is accompanied by two mechanical instabilities arising in the system due to the compressive longitudinal stress in the DA wall. The first instability leads to the appearance of the periodic pattern of alternating thinner and thicker regions; the second one leads to the buckling of the individual ECs destined to become HSPCs. We show that, while it is known that the blood flow-induced shear stress serves as a signal controlling gene machinery, it is hydrostatic blood pressure that directly stabilizes and stretches the dorsal aorta, affecting its shape during hematopoiesis. Thus, mechanical forces not only serve as signals promoting and synchronizing HSPCs production but also to assist the process by activating shape instabilities directly.  68 were maintained, crossed, raised and staged as described previously 69,70 . All animal experiments described in the present study were conducted at the University of Montpellier according to European Union guidelines for handling of laboratory animals (http:// ec. europa. eu/ envir onment/ chemi cals/ lab_ anima ls/ home_ en. htm) and were approved by the Direction Sanitaire et Vétérinaire de l'Hérault and Comité d'Ethique pour l'Expérimentation Animale under reference CEEA-LR-13007.
Microscopy. Time-lapse imaging was performed using Nikon Spinning Disk at 20X magnification as described in previously 71 . Embryos were anesthetised with tricaine (0.016%) and mounted on a glass cover dish with 0.7% low melting agarose and covered with standard E3 medium supplemented with tricaine and 1-phenyl-2-thiourea (PTU) (0.003%) to prevent pigment formation. Temperature was maintained at 28 °C by placing the dish in a temperature-control chamber during time-lapse acquisitions. Images were analysed using ImageJ.