Hydrodynamic stress and phenotypic plasticity of the zebrafish regenerating fin

How genetic information encodes a body shape is a century old scientific question (Johannsen, 1911). Over the past 100 years, a tremendous amount of knowledge has been gathered on the role of genes and epigenetics in guiding the morphogenetic programs (Ecker et al., 2018; Orgogozo et al., 2015). In addition, the impact of physics in the modulation of cell shape and tissue morphogenesis has progressively been recognized (Aegerter-Wilmsen et al., 2007; Buchmann et al., 2014; Heisenberg and Bellaïche, 2013; Vogel, 2013). A prime example of physical feedback on organ architecture is the adaptation of bone to altered load (Wolff, 1986). The locomotory organs, either on land or in water, are particularly exposed to external modulations, because they interface with forces and movements. However, continuous analysis of mechanotransduction during the development of complex vertebrate organs remains very challenging, thus limiting the progress of this research area. Consequently, it remains unknown to which extent morphogenesis can be influenced by the mechanical environment.

The caudal fin allows fish to generate propulsion power and sophisticated maneuverability (Flammang and Lauder, 2009). An advantageous feature of this appendage is its robust regenerative capacity throughout the entire fish life (Pfefferli and Jaźwińska, 2015; Sehring and Weidinger, 2020). The caudal fin is also experimentally convenient, as it is easily accessible for surgery, imaging and morphometric analysis during the regrowth process. Given its major contribution to swimming, the caudal fin is particularly susceptible to hydrodynamic forces at the time of regeneration. Overall, this appendage thus provides an ideal experimental model to study the interface between physics and biology.

In the zebrafish Danio rerio, the caudal fin has a fan-like shape with a superficially symmetrical dorsal and ventral lobe, and a central cleft (Fig. S1A,A′) (Bird and Mabee, 2003). The fin is a compound organ made of 16 to 19 principal rays with segmented bones, spanned by soft inter-ray tissue. During fin elongation, bones can branch dichotomously one to three times, so that the distance between adjacent rays is nearly constant throughout the fin surface (Fig. S1A,B′). This regularity of rays is thought to provide an optimal rigidity of the appendage during swimming. Together, the bone matrix, cartilaginous actinotrichia fibers, surface geometry and bifurcation pattern contribute to the biomechanical properties of this main locomotory appendage (König et al., 2018; Shoele and Zhu, 2009).

Upon amputation, each cut ray establishes its autonomous blastema that acts as an independent signaling center during regeneration (Jaźwińska et al., 2007; Pfefferli and Jaźwińska, 2015; Sehring and Weidinger, 2020). After a few weeks, the form of the fin is fully reproduced; however, new ray bifurcations are distally shifted relative to the original developmental position (Fig. S1D,D′), suggesting a morphogenic plasticity of this anatomical feature (Azevedo et al., 2011). Further evidence for a susceptibility of ray branching to modulations derives from a finding that the induction of ray bifurcation can be suppressed by physical separation of adjacent rays during regeneration (Marí-Beffa and Murciano, 2010; Marí-Beffa et al., 1999; Murciano et al., 2002, 2007). These observations indicate the role of ray-spanning tissue for branching morphogenesis. Consistently, the ventral- and the dorsal-most principal ray, which lack any soft tissue along one side, have never been observed to bifurcate in teleost fish (Schultze and Arratia, 2013) (Fig. S1A,D). However, at the stage of bone branching, no evidence has been reported that the inter-ray tissue would be a carrier or a source of specific growth factors for regulation of ray morphogenesis at the necessary space and temporal scales. How the inter-ray tissue regulates the bone bifurcation process remains unknown.

We previously showed that the fluid dynamics might provide cues for the ray–inter-ray distribution of sensory cells in the regenerating zebrafish fin, suggesting the impact of hydrodynamics on tissue patterning (König et al., 2019). In mammals, shear forces can be sensed through cellular mechanotransduction by osteocytes (Noda, 2011). Zebrafish osteoblasts have been suggested to have a similar role in mechanotransduction (Suniaga et al., 2018) and a swim-training experiment with zebrafish larvae revealed spatio-temporal changes in the animal's skeletogenesis (Fiaz et al., 2012). In other fishes, it has been proposed that external mechanical forces are being transduced by the inter-ray tissue to the underlying mechanosensing cellular components of the rays (Aiello et al., 2017), hence modulating the biomechanical structure of the growing locomotory appendage and guiding the ray pattern. Consequently, ray bifurcation is an attractive morphometric variable to study how the biomechanical properties of a locomotory appendage are refined by external mechanical forces.

Throughout ontogeny and regeneration, zebrafish experience an important transition in flow regime (McHenry and Lauder, 2006; Müller et al., 2000). The flow regime of a flapping fin can be defined based on its Reynolds number (Re=uL/v), which is defined by the average external flow speed (u in m s⁻¹), the fin characteristic length (L in m) and the kinematic viscosity of the fluid (ν in m² s⁻¹). Larvae and juvenile zebrafish operate mostly in the viscous regime (Re<300), where viscous shear forces acting tangentially to the surface dominate. Ontogenic growth induces the transition of zebrafish through the intermediate regime (300<Re<1000) and into the inertial one Re>1000)_, where forces perpendicular to the fin dominate. They are at least one order of magnitude larger than viscous shear forces and they offer the main contribution to propulsive forward forces (thrust). During regeneration, the gradual expansion of the fin surface with supportive rays implies related hydrodynamic changes, where the ratio of inertial versus viscous forces increases progressively, as the length (and thus Re) of the appendage increases. Hydrodynamic conditions during regenerative growth could provide positional information to regrowing tissues and thus guide the differential growth rates and bony rays branching patterns.

Measurements of spatially and temporally resolved hydrodynamic forces can be obtained by combining surface reconstructions of the analyzed object with high-resolution velocity maps of the surrounding flow of water (Dagenais and Aegerter, 2021). Flow fields can be quantified by tracking at very high resolution tracing particles seeded into the liquid, an experimental approach called particle tracking velocimetry (PTV) (Dracos, 1996; Maas et al., 1993; Pereira et al., 2006). Such reconstructions of fluid velocity fields have allowed researchers to evaluate the vortex wake and propulsive forces of synthetic foils (Bohl and Koochesfahani, 2009; David et al., 2012, 2017; Godoy-Diana et al., 2008; Green et al., 2011; Kim and Gharib, 2010; Marais et al., 2012; Muir et al., 2017; Shinde and Arakeri, 2014; Triantafyllou et al., 2000) and of elaborate bioinspired fin replicas (Dagenais and Aegerter, 2020; Dewey et al., 2012; Esposito et al., 2012; Ren et al., 2016a,b; Tangorra et al., 2007, 2010). Hence, properly tuned biomimetic hydrofoils can provide useful information about the external mechanical forces acting on a caudal fin.

In this study, we used caudal fin regeneration experiments in free swimming zebrafish in combination with hydrodynamic force quantification on synthetic hydrofoils to test the hypothesis that fluid forces are responsible for morphological plasticity in regenerating fins. We varied the biomechanical properties of zebrafish caudal fins via surgical and teratogenic manipulations to create artificially narrowed fins. We then quantified the morphometric properties of these fins to observe phenotypes in the overall shape and in ray branching patterns after complete regeneration of the narrowed appendages. In parallel, we used flexible synthetic fins, which reproduce the flow regime of caudal fins, to quantify the fluid forces acting on their surfaces. We examined the different force components (parallel and perpendicular to the surface) on biomimetic hydrofoils of varying dimensions. By using synthetic fins of increasing length, we mimicked the elongation process observed during regeneration, whereas the use of different widths allowed us to reproduce the effect of narrowing the caudal fins.

We predicted that modifying the geometry and biomechanical properties of the stump would alter the hydrodynamic forces acting on the forming tissue (Dagenais and Aegerter, 2021; Feilich and Lauder, 2015; Low and Chong, 2010; Lucas et al., 2017; Matta et al., 2020; Valdivia y Alvarado and Youcef-Toumi, 2005), and thus modulate the structure of the regenerated appendage. Therefore, varying the surface area of the caudal fin through the ablation of specific fin rays should induce changes in the bifurcation pattern and final length of the rays after regeneration. Furthermore, we expected that quantifying the surface distributions of hydrodynamic forces on biomimetic hydrofoils will highlight parallels with the morphological properties of the fin. Correlating these biological and physical data, should thus allow us to identify force components that are compatible with this purported mechanical feedback.

The following strains of zebrafish (Danio rerio) were used in this study: wild type AB (Oregon strain), Tg(-2.2shha:gfp:ABC) (Ertzer et al., 2007), referred to as shh:GFP; (OlSp7:mCherry^zf131) (Spoorendonk et al., 2008) and osterix:mCherry. Animals were maintained at 26.5°C.

Fin amputations of adult fish (females and males; 5–12 months old) and chemical treatments were performed as described previously (König and Jaźwińska, 2019). See supplementary Materials and Methods for further details. We devised three complementary approaches to obtain narrower fins compared with normal shape after uniform amputation (Figs 1 and 2). The first one is a surgical procedure, termed stepwise amputation, where a uniform amputation is followed by a more proximal amputation (hereafter termed ablation) of the three dorsal and three ventral outermost rays (Fig. 1C and Fig. 2C). The second one is a non-invasive transient treatment of the blastema with retinoic acid (RA, Fig. S2), which has been shown to deform the regenerating fin by contracting the wound epidermis (Ferretti and Géraudie, 1995; Géraudie et al., 1994; White et al., 1994). For this experiment, we used double transgenic fish shh:GFP, demarcating the signaling center of the wound epithelium (Armstrong et al., 2017) and osterix:mCherry, highlighting activated osteoblasts (Knopf et al., 2011; Thorimbert et al., 2015). These transgenic fish permit live imaging the early dynamics of the regeneration process showing the transcription of the sonic hedgehog (shh), which is required for ray bifurcation, together with the presence of bone-forming activated osteoblasts. Following a pulse treatment with chemicals, where solvent DMSO was used as control, fins of both groups were allowed to regenerate in a natural way for 21 days (Fig. S2A–F). In the third approach, we surgically manipulated caudal fins to create artificially narrower regenerates by ablating a number of principal rays in several different patterns (Fig. 1D–G and Fig. 2D–G). In the clipped experimental group, the blastema of the three most ventral and the three most dorsal rays were repeatedly ablated throughout regeneration (Fig. 1D and Fig. 2D). In the trimmed one, five rays on each side were similarly prevented from contributing to the regenerated locomotory appendage (Fig. 1E and Fig. 2E). In the quadruplets (Fig. 1F and Fig. 2F) and doublets (Fig. 1G and Fig. 2G) regenerates, groups of four or two rays, respectively, were surgically isolated by the ablation of the neighboring rays. Following amputation, the fins were allowed to regenerate for 21 days (Fig. 2). In the clipped, the trimmed, the quadruplets and the doublets regenerates, some portions of the fins were prevented from regrowing during the whole regeneration process in order to maintain the narrowed geometries. In those cases, the ablation of specific rays was repeated every 3 to 4 days. Experimental design with the timeline post-amputation for all the surgical and chemical treatments are provided in Fig. 2Ai.

All animal experiments were approved by the cantonal veterinary office of Fribourg.

Images for the following 9 experimental groups (Figs 1 and 2) were analyzed: intact fins (I, n=31, Fig. 1A and Fig. 2A), uniform amputations (U, n=24, Fig. 1B and Fig. 2B), 0.1% DMSO-treated control uniform amputations (O, n=4, Fig. 1B, Fig. S2A), retinoic acid-treated uniform amputations (R, n=5, Fig. 1B, Fig. S2D), stepwise amputation (S, n=5, Fig. 1C and Fig. 2C), clipped (C, n=33, Fig. 1D and Fig. 2D), trimmed (T, n=6, Fig. 1E and Fig. 2E), quadruplets (Q, n=8, Fig. 1F and Fig. 2F) and doublets (D, n=7, Fig. 1G and Fig. 2G). To compensate for the variability in ray number present in zebrafish caudal fins, only the seven outermost dorsal and ventral rays were considered for ray-level analyses.

Four lengths were measured for each given fin x of our experimental set (Fig. 1A and Fig. 2A): (1) ⁠, the length of ray i in fin x; (2) ⁠, the length of ray i until its first bifurcation point; (3) ⁠, the distance between the distal ends of dorsal and ventral rays I; (4) ⁠, the peduncle width. In addition, we define B_x as the set of non-ablated bony rays of fin x. Using these measures, we calculated five morphometric properties of fin x: (1)⁠, the bifurcation position for ray i of fin x; (2) ⁠, the longest non-ablated ray of fin x; (3) ⁠, the shortest non-ablated ray; (4) , the cleft depth; and (5) , the width of fin x.

Using these properties, we define the three measures used in Fig. 4A–C, normalized to the fish peduncle's width to compensate for any variation in animal size and divided by the average value in intact fins, as follows: relative fin length (Fig. 4A); relative fin width (Fig. 4B); and relative cleft depth (Fig. 4C). Fig. 4D shows the bifurcation position as defined above. The two measures used in Fig. 5 are calculated per ray i of fin x as follows: normalized ray length (Fig. 5A) and fraction of unbifurcated rays (Fig. 5B).

Statistical significance was calculated using Student's two-tailed t-tests for real values and using a Fisher's exact test for binary values. Significance is represented as *P<0.05, **P<0.01, ***P<0.001. Image analysis were performed manually using Fiji (Schindelin et al., 2012) and the GNU Octave programming language (Eaton, 2012) with custom scripts was used to process the data. See supplementary Materials and Methods for further details.

We used a three-dimensional, three components particle tracking velocimetry approach (3D-3C PTV) to quantify flow in a custom-made recirculating system composed of two large tanks connected to both sides of a 8.5×8.5×12 cm³, flow chamber (Fig. 3A,B). See supplementary Materials and Methods for further details. A pump pushed water from the upstream tank into the flow chamber through a flow straightener to obtain a laminar stream of controlled fluid velocity. 3D fluid velocity vector fields were recorded as previously described (Lai et al., 2008; Pereira et al., 2006; Pothos et al., 2009).

Reconstruction of the particles' positions and velocities was performed using the V3V software (Lai et al., 2008; Pothos et al., 2009). A 2D Gaussian fit of the particles intensity distributions in each image was used to identify the 2D positions maps. Triplets of position fields were captured simultaneously by the three cameras and superimposed to compute the 3D positions by triangulation. Our temporal resolution resulted in an average displacement of about 10 pixels between subsequent time frames. 3D fluid velocity vector fields were then calculated between two consecutive position fields using a particle tracking algorithm based on a relaxation method (Lai et al., 2008; Pereira et al., 2006; Pothos et al., 2009). Raw velocity fields were filtered to remove outliers using three different filters: a velocity range filter (where vectors whose components exceed ±0.25 m s⁻¹ are removed), a Gaussian smoothing filter (σ=1.5 voxels) and a spatial median filter (with neighborhood size=5 mm). Filtered 3D velocity vectors were then interpolated onto a regular Cartesian grid with an isotropic spatial resolution of 0.75 mm (Fig. 3E).

We designed a set of biomimetic hydrofoils that recapitulate the most relevant features of our experiments on caudal fin regeneration (Fig. 3C, Table 1). Three foils with increasing lengths (60%, 80% and 100%) were used to study how forces vary during morphogenesis. Foil A1, the shortest of these foils, replicates the fin geometry right after uniform amputation by following the proportions of a zebrafish caudal from its peduncle to the amputation plane. Foil A2 corresponds to an intermediate length, slightly proximal to the level of ray bifurcation. Foil A3, the longest one, represents a regenerated fin [at 21 days post ablation (dpa)], such that ray bifurcations would in principle occur during the transition from A2 to A3. To incorporate our observations with altered fin width, we designed foil B to be 70% wider than our reference foil A3. Finally, foils C and D were created to assess the influence of the cleft on the distribution of hydrodynamic forces, in particular when the distal edge undergoes a transition from flat, to slightly rounded and to bilobed during regeneration.

The six different hydrofoil designs (Fig. 3C) were produced as follows. A rigid hollow cast was 3D printed, liquid polydimethylsiloxane (PDMS) was poured inside the cast, a rod was inserted at the base of the hydrofoil and the material was cured for 36 h at 58°C. The resulting flexible membranes (Fig. 3D) have a thickness of 0.55 mm and a Young's modulus of approximately 0.8 MPa (Puri et al., 2018). The hydrofoils were used independently in the flow chamber. They were mounted horizontally inside the chamber and actuated through a dedicated opening on the back wall by a servomotor fixed outside the chamber (Fig. 3B). The flapping movement was a pitching sinusoidal motion with its rotation axis parallel to the z-axis and perpendicular to the camera plate. The exact position of the foil was determined by image analysis. In each time frame, 25 points were manually drawn on the midline of the fin and fitted with a quadratic function. The 3D surfaces of the hydrofoil were then reconstructed from this fitted midline. The reconstructed fin surface is at most 0.7 mm away from the real surface because particles located at the fluid-solid interface cannot be resolved. This distance is smaller than the spatial resolution, allowing the shear stress patterns to be properly captured at the surface.

We use two dimensionless parameters to determine the flow regime acting on our flapping fin-like hydrofoils and to encompass the propulsion dependence on oscillations: Re and the Strouhal number (St=fA/u), where f is the flapping frequency and A the amplitude of the oscillations. More theoretical background on these parameters and the characterization of flow regimes can be found in Prandtl (1952) and Whitaker (1968). Given that v=10^–6 m² s^–1 for water at room temperature, we tuned the values of the other parameters (f=3 Hz, u=0.055 m s⁻¹) and the characteristics of the hydrofoils (Fig. 3C, Table 1, 0.015≤L≤0.025 m, A inferred by image analysis) to fit biologically relevant values of Re and St when compared with those measured for zebrafish caudal fins (Mwaffo et al., 2017; Palstra et al., 2010; Parichy et al., 2009; Taylor et al., 2003; Triantafyllou et al., 2000).

Time-resolved stress maps on the surfaces of flexible fins were obtained as recently described (Dagenais and Aegerter, 2021). Briefly, the viscous stress tensor at time t is a function of the spatial derivatives of the velocity field: ⁠. This tensor was evaluated on each vertex of the 3D cartesian grid using a finite central difference scheme. τ_ij was then projected on the surface of the solid fin to obtain a stress vector, which is the total force per unit area generated by the fluid on the foil. Uniformly spaced unit vectors normal to the surface of the foil, and pointing outwards, were calculated based on the vertices of the 3D reconstructed fin, which consists of 25 nodes in the dorso-ventral direction, 25 nodes in the anterior-posterior direction and 3 nodes in the left-right direction. A coordinate system local to each node was defined as (dv, ap, lr), where dv is the local dorso-ventral axis, which is coplanar with the foil's side surface and parallel to the base-tip midline, ap is the anterior-posterior axis, which is perpendicular to dv and also coplanar with the side surface, and lr is the left–right axis, which is normal to the side surface. Each component of the stress vector acting on a surface with a unit normal vector (oriented outwards) is expressed according to ⁠. The cartesian stress vectors (S_x, S_y, S_z) were finally decomposed into biologically relevant components, by projecting them onto the (dv, ap, lr) directions, which gives a total of six different stress components. These are the three components on the side of the foil in the local coordinate system, i.e. ⁠; ⁠; ⁠, as well as on the distal tip surface of the foil, i.e. ⁠; ⁠; . As an example, the distribution of is shown in Fig. 3F. The components on the dorsal and ventral surfaces were not considered in this analysis because of their negligible involvement during morphogenesis.

These stress maps were evaluated at 30 time points equally distributed over 3 flapping periods, and the results were time-averaged over these 3 periods. Spatial averaging was also performed following the dorso-ventral symmetry of the hydrofoil to reduce experimental uncertainties. This averaging procedure was performed for both the absolute values and the signed values of the stresses, yielding a total of 12 different stress measures to evaluate. Signed values provide a measure of the resulting internal tension (i.e. of the excess of contraction or stretching at the surface). Owing to the symmetric shape, spatial averaging yields to 0 for the stresses and ⁠. To determine the tension in the dorso-ventral direction, we calculated ⁠, where l is the position along the dorso-ventral axis and Δl the spatial step size between nodes on the surface of the foil. With this, positive values of and (with units μN m⁻¹) indicate contraction, whereas negative values indicate stretching. has a period average equal to zero, due to the flapping motion which is symmetrical in the left–right direction. This component is therefore left out of the analysis. Hence, 11 stress components were precisely quantified: ⁠; ⁠; ⁠; ⁠; ⁠; ⁠; ⁠; ⁠; ⁠; ⁠; ⁠.

A detailed analysis of experimental uncertainties and limitations in resolution for our experimental set-up was previously published (Dagenais and Aegerter, 2021). See supplementary Materials and Methods for further details.

Using the proportions of uninjured fins (Fig. S1A), we classified the 25 nodes of our stress maps along the dorso-ventral axis as follows: 2 nodes for the dorsal rim, 8 for the dorsal lobe, 5 for the cleft, 8 for the ventral lobe, 2 for the ventral rim. Using this definition, we can formalize averages of the stress components in different regions of the foil as follows: ⁠, where h is the hydrofoil, which S(i) is measured on, region refers to the lobe, cleft or rim and i_end and i_start are the beginning and end points of the above defined regions.

Stress components regulating growth rates should be compatible with higher growth rates in shorter fins. This implies that the force distribution should either decrease with increasing fin length (A1<A2<A3) if the signal is an activator of growth, or increase if the signal is an inhibitor of growth. In other words, a monotonic dependence needs to be present in order to lead to growth deceleration, which we formalize with our ‘growth deceleration’ mathematical criterion: [ave(A1_{rim+lobes+cleft})>ave(A2_{rim+lobes+cleft})>ave(A3_{rim+lobes+cleft})] for a growth activator or [ave(A1_{rim+lobes+cleft})<ave(A2_{rim+lobes+cleft})<ave(A3_{rim+lobes+cleft})] for a growth inhibitor.

Growth-inducing mechanical signals should also be compatible with a different average amplitude between the lobes and the cleft of the fins to induce a bilobed shape. This difference (signal in the lobes minus signal in the cleft) should be positive if the signal promotes growth (activator) and negative for an inhibitor of growth. Moreover, the difference in mechanical signal between the lobes and the cleft should be higher in a wide fin compared to a narrow fin, to be consistent with the measured reduction of cleft depth in narrower fins (Fig. 4C). We formalize both conditions with our ‘cleft reduction’ mathematical criterion: |[ave(B_lobes)–ave(B_cleft)]|>|[ave(A3_lobes)–ave(A3_cleft])|.

Stress components regulating the bifurcation of rays should be compatible with the average induction of bifurcations, represented by hydrofoils A2 and A3. This implies that the force signal would need to increase from shape A2 to A3 if the signal activates bifurcation or decrease from A2 to A3 if it inhibits bifurcation. Because narrower fins are associated with a distalization of bifurcations (Fig. 4D), this force signal should also be lower in narrower fins (A3 versus B) for an activator or vice-versa for an inhibitor. We formalize both conditions with our ‘bifurcation position’ mathematical criterion: [ave(A2_lobes+cleft)<ave(A3_lobes+cleft)<ave(B_lobes+cleft)] for a signal promoting bifurcation or [ave(A2_lobes+cleft)>ave(A3_lobes+cleft)>ave(B_lobes+cleft)] for a signal inhibiting bifurcation.

Bifurcation-inducing mechanical signals should also be consistent with the absence of bifurcations in the outermost rays of regenerated fins. We formalize this condition with our ‘unbifurcated outer rays’ mathematical criterion: ave(A3_lobes+cleft)−ave(A3_rim) where a positive value indicates an activator and a negative value an inhibitor.

For visualization purposes, averages were normalized to the corresponding values on hydrofoil A2.

To study how varying the surface area of the caudal fin modulates the structure of the regenerated appendage, we quantified the morphometric properties of three sets of artificially narrower fins after regeneration (Fig. 2A,Ai).

We first compared uninjured fins, uniform amputation, stepwise amputation and RA pulse-treated fins. Fins of all groups were allowed to regenerate in a natural way for 21 days (Fig. 2B,C, Fig. S2). At 21 dpa, morphometric analysis of fin regenerates indicated that the fins recreated more than 80% of the original intact length, being close to the termination of regeneration (Fig. 4A). In the stepwise group, we observed a significant 20% fin narrowing (Fig. 4B), a 20% shallower cleft (Fig. 4C) and a comparable distalization (Fig. 4D) with respect to uniform amputations. In each fin of the RA group, some of the rays have fused following the contraction of inter-ray tissue (Fig. S2F), leading to a chemical ablation of rays. Interestingly, we measured a 50% fin narrowing (Fig. 4B), a 20% shallower cleft (Fig. 4C) and a significant further distalization (Fig. 4D) in this group compared with the uniform group. Analysis of single ray lengths showed that treatment with RA caused the 4 dorsal and the 5 ventral outermost rays to be significantly shorter (Fig. S2G). For the central rays, there was no significant shortening, indicating that the shallower cleft is mainly caused by changes in the outer rays.

Second, given that chemically ablating rays showed an enhanced phenotype, we quantified fins with mechanically ablated rays. At 21 dpa, the clipped and trimmed fins restored their expected length, compared with uniform and stepwise amputated fins (Fig. 4A and Fig. 2D,E). In addition, these ablations were able to further narrow the relative width of the regenerates to a similar extent as the contraction observed following RA treatment (Fig. 4B). In the clipped group, we observed a shallow cleft depth comparable to both RA treatment and stepwise amputation (Fig. 4C), but a distalization phenotype intermediate to these two groups (Fig. 4D). In the trimmed group, we detected a further reduction in cleft depth and an increased distalization with respect to clipped fins (Fig. 4C,D). Analysis of single ray lengths showed that the reduced cleft depth in the clipped group was caused by a significant reduction of lengths of the outer rays (Fig. 5A), similar to the observations in the RA treated fins. In the trimmed group, the outer rays also seemed to contribute most to the observed reduction, although they were not significantly shorter (Fig. 5A).

Third, we reduced to the maximum the width of the regenerates by separating the rays of the stump into quadruplets and doublets (Fig. 1F,G and Fig. 2F,G). Comparison of relative ray length at 21 dpa to uniform amputation confirmed that all bony rays in the quadruplets/doublets experiments regenerated to their expected length (Fig. 5A), with the exception of the two central doublets. We observed an apparently random uneven length between various pairs of these central doublets, likely a result of spontaneous breakage of one of the regenerates due to the fragility of these thin structures.

Quantifying the fraction of bifurcation events for individual ray per experimental group showed that all central rays of mechanically altered fins bifurcate during regeneration (Fig. 5B and Fig. 6H), while the outermost principal rays very rarely do so (Fig. 5B). Inhibition of bifurcation was confirmed at the molecular level in transgenic fish shh:GFP, which normally display a split shh expression pattern proceeding ray splitting (Fig. S1E–H). We found that the shh expression domain of the new unbifurcated outer rays surrounding the ablated regions did not split (Fig. 6A–F).

To investigate how hydrodynamic forces are modulated by variations in the geometry of the fin, we measured the surface distributions of the 11 possible components of hydrodynamic forces on 6 biomimetic synthetic foils with different morphological properties. To identify forces compatible with a role in mechanoregulation, we classified each stress component as a candidate activator, candidate inhibitor or as uncorrelated with respect to four key morphometric properties observed in our regenerated fins.

Growth deceleration, the first criterion which formalizes enhanced regeneration speed in proximal cuts, finds and as putative activators, ⁠, ⁠, ⁠, ⁠, ⁠, and as putative inhibitors as well as and as uncorrelated stress components (Fig. 7A). Cleft reduction, which is based on the observed shallower cleft depth in narrower fins (Fig. 4C), highlights ⁠, and as putative activators, no putative inhibitors and the 8 other stress components as uncorrelated (Fig. 7B). The bifurcation criterion combines the average positioning of bifurcations at 83% of the total length at 21 dpa (Fig. 4D) with the observed distalization of bifurcations in narrowed fins (Fig. 4D) to classify and T^t_dv as putative activators and all remaining stress components as uncorrelated (Fig. 7C). Finally, formalizing the unbifurcated outer rays (Fig. 5B) identifies ⁠, ⁠, ⁠, ⁠, ⁠, and as putative activators, and as putative inhibitors as well as and as uncorrelated (Fig. 7D).

In order to be a candidate growth regulator, a mechanical signal must consistently act as an activator (or an inhibitor). This was only the case for which would act as a growth activator according to both criteria (Fig. 7A,B). Similarly, for a candidate regulator of bifurcation, only was a suitably consistent signal for activating bifurcations in regenerating fins (Fig. 7C,D).

Having identified the dorso-ventral stress at the tip of the fin and the resulting tension as candidates for the regulation of morphological changes during growth, we studied their properties in more detail.

The period-averaged absolute value of the dorso-ventral shear stress at the tip indicates how the putative growth activating signal is distributed over the fin tip (Fig. 8A). The profile displays a bilobed signature in A1, A2, A3 as well as B (Fig. 8A), with higher signal in the locations of the prospective lobes and lower magnitude in the center as well as on the rims. The magnitude of decreases as the fin model elongates (A1 to A2 to A3), and increases as the fin widens (A3 to B) (Fig. 8A). Moreover, the difference in stress magnitude between the lobes and the cleft is more important for the wider foil compared with the narrower version (Fig. 8A).

The resulting internal tension (Fig. 8B) is the putative bifurcation regulator and represents dorso-ventral contractions of the tissue when positive and conversely stretching when negative. We observed a length-dependent transition from stretching (in A1) through neutral (in A2) to contraction (in A3). Moreover, we measured extremely small amplitudes of tension at the location of the outermost non-bifurcating rays as well as an overall lower tension in narrower fins (Fig. 8B).

Transitions from a uniform amputation (A2) to a rounded (foil C) and finally bilobed (foil D) shape show that both and transiently increase in shape C (Fig. 8).

The caudal fin yields powerful as well as complex locomotion and is thought to be one of the most important evolutionary innovations in the aquatic environment (Flammang and Lauder, 2009). In the present work, we explored the role of mechanical forces during its morphogenesis as recapitulated by regeneration. Through a variety of surgical and chemical treatments, we have shown the phenotypic plasticity of the biomechanical structure of the fin. RA treatment induced the fusion of numerous rays (Fig. S2F), thus producing fins with a reduced number of rays. The obtained morphologies were comparable to those resulting from surgically manipulated fins with a similar number of remaining rays (Fig. 4). In particular, we observed a progressive distalization of the ray bifurcations and reduction of the depth of the cleft in ablated regenerates. These changes could be adaptations aimed at improving the thrust or efficiency of the locomotory appendage by modulating the flexibility of the fin (Zhu and Bi, 2017). Such fine-tuning highlights the feedback mechanism between abiotic forces and locomotory organs. Consistently, the role of mechanical factors has also been linked with plasticity of the fish skeleton (Grünbaum et al., 2012). Several studies demonstrated that water velocity and exercise, which exert mechanical stress, influence body shape of fish, increase the number of bone-forming osteoblasts and result in higher levels of bone mineralization (Fiaz et al., 2012; Gunter and Meyer, 2014; Palstra et al., 2010; Suniaga et al., 2018). Although osteoblasts have been proposed as good candidates for the detection of mechanical loads (Witten and Hall, 2015), deciphering the mechanisms of cellular mechanotransduction, or how exactly cells transform mechanical stimuli into electric or chemical signals, constitutes a complex and unresolved problem (Suniaga et al., 2018).

Up to now, researchers have devoted most of their attention to the molecular and genetic processes related to positional information in the fin (Sehring and Weidinger, 2020; Wehner and Weidinger, 2015). Various grafting experiments have shown that position-dependent mechanisms are important for the control of the bifurcation process and the determination of individual ray lengths (Murciano et al., 2018, 2002, 2007; Shibata et al., 2018). Rays grafted into a new position within the fin adopt length characteristics of the host region following regeneration, although elements of positional information also exist within the rays, independently of their host environment (Shibata et al., 2018). Another striking sign of positional information is the growth rate dependence on the proximo-distal location of the amputation plane, where more proximal structures regrow faster (Lee et al., 2005; Monteiro et al., 2014; Uemoto et al., 2020). Our study further extends these findings and indicates that biological transformations occurring during fin morphogenesis may be induced by changes in the surface distribution of hydrodynamic stresses acting on the growing blastema.

Narrower fins created by the ablation of the outermost rays undergo a shortening of the remaining lateral rays and thus present a less-pronounced cleft (Fig. 4C). In these artificially thinner appendages, frequent bifurcation failure is observed in rays that normally bifurcate, most prominently in the outer rays deprived of their lateral neighboring rays, but also in the central rays (Fig. 5B). The fraction of unbifurcated rays supports previous reports that the outermost principal rays do not bifurcate (Marí-Beffa et al., 1999; Schultze and Arratia, 2013). The absence of shh expression domain splitting of the new outer rays adjacent to the ablated regions further highlights the substantial phenotypic plasticity of the blastema upon altering its mechanical conditions (Fig. 6A–F). However, our quadruplets/doublets experiments also showed that while these new outer rays rarely bifurcate, they occasionally do so (Fig. 5B and Fig. 6H). These results suggest that mechanotransduction of external forces through the inter-ray tissue contributes to splitting of shh expression and that the presence of inter-ray tissue on both sides of a ray is not strictly required for the induction of ray branching events. Moreover, the narrowing of regenerating fins induces an increased distal shift of the bifurcation points compared with uniform amputations (Fig. 4D), further enhanced in quadruplets and even more so in doublets (Fig. 5B). Overall, these results suggest that caudal fin morphology is modulated by the manipulation of its biomechanical properties during regeneration.

To test the hypothesis of a mechanical control of morphogenesis via an interaction of the regrowing tissue with external hydrodynamic forces, we measured fluid forces on the surface of hydrofoils reproducing the flow regime of caudal fins. This precise quantification demonstrated that the various fin shapes are subjected to different hydrodynamic force distributions. By screening for stress components compatible with the results of our regeneration experiments, we have identified the dorso-ventral shear stress at the tip of the fin as the only one, among all hydrodynamic forces acting on the hydrofoil, displaying changes in response to length and width that are compatible with a purported mechanical feedback (Fig. 7). In particular, its absolute value (⁠⁠) is consistent with a role as growth activator and its corresponding tension (⁠⁠) with an activation of ray bifurcations. This stress component is particularly relevant in the context of fin morphogenesis, as it corresponds to the viscous forces that can contract or stretch the tip of the growing tissue.

The spatial distributions of the dorso-ventral shear stress at the tip of the fin (Fig. 8) allowed us to put in perspective the results of our regeneration experiments with mechanically altered fins. On one hand, the enhanced shear stress in the lobes of short fins (A1, Fig. 8A) could explain the formation of a bilobed trailing edge from a flat amputation plane (Fig. S1C), induced by variable growth rates of individual rays (Uemoto et al., 2020). Note that this consistency with an overall bilobed shape was selected in our cleft reduction criterion only for A3 and B, yet it is fulfilled by the signal in A1 as well as A2. Growth rate would then gradually decrease following the reduction in stress intensity for longer appendages (A2 to A3, Fig. 8A), as observed in real fins as they approach their final size (Lee et al., 2005). Simultaneously, growth would reduce as the fin reaches its bilobed shape (C to D, Fig. 8A) because the stress signal decreases predominantly in the lobes upon the transition from a rounded to a bilobed shaped fin. The increase of the lobe signal in wide fins (A3 versus B, Fig. 8A) would account for the observed decrease in cleft depth in narrower fins (Fig. 4C). On the other hand, the transition from a stretching tensile force at the tip to a contraction force as the fin model approaches its final length (A1, A2, A3, Fig. 8B) suggests a process through which mechanical cues from the fluid could induce the bifurcation of bony rays. Lower contraction in the narrower fin model (A3 versus B, Fig. 8B) offers a plausible explanation for the observed distalization of bifurcation (Fig. 4D) as the length required to reach the hypothesized contraction threshold would be greater. Moreover, this lower mechanical signal would also predict the observed increase in the fraction of unbifurcated rays at the center of the fin (Fig. 5B). Our results thus suggest that dorso-ventral shear stress at the fin tip may modulate tissue growth rates, and that ray bifurcation may be triggered by tissue contraction above a certain intensity threshold.

Our results on fins of different shapes (A2, C, D, Fig. 8) give interesting predictions on the forces acting during the regeneration of a typical zebrafish caudal fin. The simultaneous emergence of the bilobed shape together with the lengthening of the regenerate presumably leads to the emergence of stress distributions more complex than those measured on our stereotypical hydrofoils. Measurements on intermediate combinations of shape and length, both on actual fins as well as hydrofoils, would potentially bring more insights in the relation between forces, bifurcation and fin shape.

Modeling caudal fins using synthetic hydrofoils is a powerful approach as it provides a description of the fluid force distributions on the fin surfaces, which could not be measured on real fins due to technical limitations related to the size of the zebrafish caudal appendage, the high complexity of the membrane deformations during swimming and the variability in swimming patterns. Indeed, even if the flow field can be recorded around real swimming fish using PTV experiments, robustly reconstructing force fields from such turbulent flows remains an enormous experimental challenge. Our approach using biomimetic hydrofoils allowed us to precisely quantify the changes in fluid forces induced by variations in the geometry of a fin under reproducible conditions. One inherent limitation of this approach is that the active membrane deformations (other than those induced by a pitching motion at the peduncle), which could also be affected by the removal of specific bony rays, were not mimicked by the synthetic foils.

Our study illustrates how the investigation of hydrodynamic stresses based on three-dimensional PTV experiments can provide important information about the mechanical interplay between a biological system and the surrounding fluid. We anticipate that hydrodynamic stresses will be integrated in a comprehensive growth model of the zebrafish caudal fin, as the most external and thus most upstream regulator. The works of Murciano et al. (2007), Shibata et al. (2018) and Uemoto et al. (2020) support this hypothesis as the grafted rays changed their bifurcation pattern or their length based on their new location, which could be explained by local differences in hydrodynamic stress. Understanding the role of the mechanical environment in the (re)establishment of positional information is of fundamental importance in the fields of regenerative medicine and tissue engineering.

By correlating our biological and physical data, we conclude that viscous shear forces acting on the distal edge of the fin and the resulting internal tension may contribute to the regulation of the fin's morphology. Our results indicate that the length and bifurcations of the bony rays are modulated by the geometry of the paddle, possibly via a change in hydrodynamic forces acting on the surfaces. The experimental results presented in this work thus corroborate the hypothesis of a coupling between mechanical forces and morphogenesis, and support the idea that mechanical cues are transmitted from the surrounding fluid to the fin tissue to be integrated in the genetic cascade responsible for growth processes such as ray bifurcation and the orchestration of the overall fin shape. The simplicity of the regeneration experiments combined with stress measurements on hydrofoils, and the set of results presented here pave the way for future studies on the plasticity of fins in relation to their hydrodynamic environment.

We thank V. Zimmermann for zebrafish facility management and husbandry; Sahil Puri, Siddhartha Verma, Flavio Noca and Petros Koumoutsakos for interdisciplinary discussions.