The fin-to-limb transition as the re-organization of a Turing pattern

A Turing mechanism implemented by BMP, SOX9 and WNT has been proposed to control mouse digit patterning. However, its generality and contribution to the morphological diversity of fins and limbs has not been explored. Here we provide evidence that the skeletal patterning of the catshark Scyliorhinus canicula pectoral fin is likely driven by a deeply conserved Bmp–Sox9–Wnt Turing network. In catshark fins, the distal nodular elements arise from a periodic spot pattern of Sox9 expression, in contrast to the stripe pattern in mouse digit patterning. However, our computer model shows that the Bmp–Sox9–Wnt network with altered spatial modulation can explain the Sox9 expression in catshark fins. Finally, experimental perturbation of Bmp or Wnt signalling in catshark embryos produces skeletal alterations which match in silico predictions. Together, our results suggest that the broad morphological diversity of the distal fin and limb elements arose from the spatial re-organization of a deeply conserved Turing mechanism.

D espite the remarkable diversity of animal shapes, the repertoire of genes used for morphogenesis is unexpectedly conserved between species-an observation termed deep homology 1 . How similar sets of genes can govern the generation of diverse morphologies during evolution, however, is still not well understood. The fin-to-limb transformation is a paradigmatic example of morphological evolution 1 . The distal skeletal patterns of vertebrate limbs and fins have changed multiple times during evolution (Fig. 1a), and their homologous relationships have been controversial [2][3][4][5][6][7] . According to comparative anatomy, digits are regarded as a novel structure of tetrapod limbs, and do not trace back to non-tetrapod sarcopterygian fins 2,8 (for example, Sauripterus 9 and Panderichthys 7 ). Yet despite the clear skeletal differences, recent molecular studies show unexpected similarities between the distal fins and limbs at the genetic level. For example, the digit-specific regulatory sequence of the murine Hoxa and d genes has recently been found in the genomes of the skate and the spotted gar, where they also drive similar expression in the distal fin/limb bud 6,10 . Thus a deep question remains: how can the skeletal arrangement change so markedly, when the well-known patterning genes do not?
Recent studies 11,12 have provided strong evidence that digit patterning in the mouse limb is driven by a Turing mechanism, which has long been suggested theoretically [13][14][15][16] . Specifically, it has been proposed that regulatory interactions between BMP, SOX9 and WNT form a Turing network that creates a periodic molecular pre-pattern specifying the positions of the digits (the BSW model) 12 . A Turing model can generate different types of patterns, such as spots and stripes with only slight changes in parameter values. Therefore, we explored whether the BSW model could explain the marked changes in the distal skeletal arrangement of fins and limbs.
In this study, we focus on the pectoral fin development of the catshark, Scyliorhinus canicula for two reasons: (a) its fin skeletal elements are formed by individual condensations 17 , which are similar to the condensation process of tetrapod limbs; and (b) its genome is less derived than that of teleost genomes 18 . We show that spots of Sox9 expression underlie the distal elements of S. canicula pectoral fin buds. In addition, by building a computer model, we demonstrate that such spot-like Sox9 expression can be explained by the BSW model with slight modification of its parameters. Together, our results suggest that the broad morphological diversity of the distal fin and limb elements arose from the spatial re-organization of a deeply conserved Turing mechanism.

Results
The first periodic expression of Sox9 is a distal row of spots.
To understand the dynamics of skeletal patterning in S. canicula pectoral fin buds, we first examined a time course of Sox9 expression using optical projection tomography (OPT) 19 . Our data revealed that the formation of distal nodular radials can indeed be captured in detail by Sox9 expression, which is initiated as a single row of spots along the anterior-posterior axis (Fig. 1b). In more detail: at early stages Sox9 starts to be expressed in the basal elements ( Fig. 1b-i), and in the posterior-distal region (bracket in Fig. 1b-ii). Subsequently, a curved row of spots develop (arrowheads in Fig. 1b-iii), which is initially more continuous in the posterior region, but gradually also breaks up into spots (arrowheads in Fig. 1b-iv and v). These Sox9 spots can be identified as the second row of distal nodular elements of the final skeleton (see Supplementary Fig. 1a-c for the later stages and the detailed annotation). A previous study reported a roughly similar expression pattern for Sox8 (ref. 17), but distinct spots were hard to discern as a three-dimensional imaging technique (such as OPT) had not been used. Our data therefore reveals that the first stage of radial patterning is a dynamic specification of a spot-like pattern, in contrast to the stripy Sox9 pattern of the mouse limb bud. We thus focused our study on control of this spot pattern (rather than subsequent expression of Sox9 proximally or distally) for two reasons: firstly, we are interested in the initial symmetry-breaking process responsible for the overall radial arrangement, and secondly because previous studies suggest that the mechanism of patterning the distal periodic elements shows molecular differences from those controlling more proximal elements 12,20 .
Out-of-phase patterns of Bmp and Wnt expression with Sox9. If the patterning of the S. canicula pectoral fin was controlled by a Turing system similar to that controlling mouse digit patterning 12 , Bmp and Wnt might be expressed or active in a pattern out-of-phase with Sox9 (Fig. 2a). We thus examined expression of Bmp-and Wnt-related genes in the S. canicula pectoral fin buds. Firstly, we found that Bmp2 was expressed only in the distal fin edges, whereas in mice it displays the strongest out-of-phase pattern with Sox9 (ref. 12) ( Supplementary Fig. 2a). Instead, Bmp4 was expressed in the fin mesenchyme and indeed has a pattern complementary to Sox9 ( Fig. 2b and Supplementary  Fig. 2d; white and black arrowheads indicating a row of expression gaps where Sox9 has a row of spots). Next, we examined genes related to Wnt signalling ( Supplementary  Fig. 2b,c,e). Wnt5b had a pattern out-of-phase of Sox9 ( Fig. 2c Sox9/Bmp4 or Sox9/Wnt5b ( Supplementary Fig. 2f). Thus, although several differences were found, the overall relationship between Bmp, Sox9 and Wnt is conserved from fish to mammals.
A dynamical model of S. canicula fin development. To confirm if a BSW Turing network could reproduce the early spot pattern of the S. canicula fin, we built a realistic computational model using a similar approach to our previous mouse limb model 22 .
In particular, we obtained a time course of pectoral fin morphologies, created a series of two-dimensional (2D) triangular meshes, and calculated hypothetical tissue trajectories which represent possible growth maps ( Fig. 3; Supplementary  Fig. 3). A crucial step was to determine how to align the chronological series of fin shapes ( Supplementary Fig. 4a). To constrain this configuration, we required real fate map data, and thus despite the very slow growth of S. canicula fins we performed carbon-particle-based fate mapping 23 (which required a minimum of 30 days to observe sufficient displacement of labelled tissue). By comparing real tissue displacements with in silico predictions, we could derive a realistic computational growth map, in which the posterior part of pectoral fin bud expanded more than the anterior part ( Supplementary Fig. 4b-e). Interestingly, the asymmetry in growth observed along anteriorposterior axis is consistent with the fate maps observed in chick limb buds 24,25 .
In silico modelling of the spot-type Sox9 expressions. Because of the overall conservation of the distribution of Sox9, Bmp and Wnt, we next explored whether the S. canicula distal fin elements could also be specified by the BSW model 12 , which is expressed by the following partial differential equations: where S, B and W are abstract variables representing the amounts of Sox9, Bmp and Wnt expression, respectively, k 2 to k 9 are kinetic parameters (Fig. 4a), D B and D W are diffusion constants of B and W, respectively, a B and a W are constant production terms of B and W, respectively, and b is a global coefficient that controls the speed of pattern appearance. This system is composed of one non-diffusive molecule, Sox9 (S) and two diffusive molecules, Bmp (B) and Wnt (W). Our numerical simulations revealed that the model formed spots instead of stripes when Wnt production was significantly higher than Bmp production ( Supplementary  Fig. 5a). When we simulated with this condition in the fin growth model, a uniform distribution of spots emerged that had no resemblance to the real Sox9 expression patterns ( Supplementary  Fig. 5a). Previous work in the mouse has shown that distal Hox genes and fibroblast growth factor (FGF) signalling provide spatial modulation of the Turing network to sculpt the Sox9 pattern into the normal digit arrangement 12 . We thus hypothesized that these molecules could play in a role in shaping the Sox9 expression into a curved row of spots at a certain distance from the distal fin edge. In the mouse model, Hoxd13 restricts the domain where the Turing instability occurs. Because in S. canicula, Hoxa13 instead of Hoxd13 is significantly expressed in the distal fin buds 17,26 , we first examined Hoxa13 expression with OPT, but found that the Hoxa13 expression domain did not overlap significantly with the distal expression of Sox9 (Fig. 4b). Hence, we ruled out a similar role of Hoxa13 in S. canicula fin bud.
We therefore asked whether Fgf alone could control the Sox9 spot pattern. We simulated an Fgf gradient by assuming that the ligand is produced at the distal fin edge and diffused towards the proximal part (Fig. 4c). The shape of the resulting gradient was roughly similar to the expression domain of an Fgf target gene, Dusp6 (ref. 27) (Fig. 4c). Next, we assumed that the Fgf gradient modulated the BSW model by regulating the same parameters as in the mouse model-repressing k 4 and boosting k 7 -which made the system pass through the Turing space from proximal to distal ( Fig. 4c; and see Methods for the equations). When the BSW model was simulated under the influence of the Fgf gradient, a curved row of Sox9 spots was formed at a certain distance from the ectoderm (Fig. 4d). More specifically, the dynamic pattern shared two features with the observed time course: (a) it started at the anterior and posterior ends (which are also the more proximal positions) and gradually extended distally, and (b) it initially showed some connected regions of expression, which then broke up into a series of spots (compare Fig. 4d with Fig. 1b). To evaluate the role of growth in the model, we also simulated the BSW network on a static fin model, and found that without growth the Sox9 pattern broke into spots more slowly than with growth ( Supplementary Fig. 5b), suggesting that growth may contribute to reliable spot separation. In addition, both Bmp and Wnt showed strong expression in the distal region and a series of expression gaps which correspond to the spots of Sox9 (Fig. 4e), consistent with the experimental data. The relatively shallow predicted interdigital Wnt distribution was also consistent with the real expression pattern of Wnt5b. Therefore, the model qualitatively reproduced the expression patterns of Sox9, Bmp4 and Wnt5b in S. canicula fin buds. Our computer model reflects the normal Sox9 patterning, but could it correctly predict the main features of experimental perturbations? A clear prediction of the model is that if Fgf signalling is reduced, the position of Sox9 spot expressions will move closer to the distal fin edge (Fig. 4f). To test this prediction, we treated S. canicula embryos with the Fgf receptor inhibitor SU5402 (ref. 28), and confirmed the efficiency of inhibition by qPCR of the target gene Dusp6, Supplementary Fig. 5c). The resulting Sox9 pattern showed some variability (losing its periodic form and losing expression in the anterior fin-discussed further in the legend accompanying Supplementary Fig. 5d) but it was frequently shifted distally, consistent with the prediction (Fig. 4g; Supplementary Fig. 5d for details). Thus, the row of Sox9 expression spots appears to be positioned by Fgf signalling.
Experimental tests for in silico model predictions. We wished to test two other molecular perturbations: inhibitions of Bmp and Wnt, to see if our model predictions would match with in vivo experiments. Firstly, we performed numerical simulations with decreasing values of k 2 , which represents inhibition of Bmp signalling. The simulation showed two features: the distal-most Sox9 spots failed to form, and those which did form were smaller (compare Fig. 5b with Fig. 5a). To carry out experimental perturbations, we treated embryos with inhibitors about 4 days before the distal Sox9 expression appears, and checked the efficiency of inhibition by qPCR analysis and in situ hybridization on target genes ( Supplementary Fig. 6a,b,c). Consistent with this prediction, S. canicula embryos treated with a Bmp inhibitor LDN-193189 (ref. 29) showed a loss of some or all of the Sox9 spots (compare Fig. 5d and e). To assess the skeletal patterns of Bmp inhibitor-treated embryos, we carried out longer treatments on embryos and cultured them more than 1 month. Interestingly, this long-term treatment sometimes resulted in an expansion of the apical ectodermal ridge-like structure and the width of the pectoral fin buds ( Supplementary Fig. 6d), suggesting that the previously reported negative effect of Bmp on the chick AER 30 is also conserved in catshark fin buds. Cartilage staining clearly showed that the posterior nodular elements were lost and sizes of the remaining spots were smaller (black arrowheads in Fig. 5h) than in control fins (Fig. 5g), as seen in the simulation.
We next examined Wnt inhibition. In the model, this perturbation was performed by decreasing the Wnt production term. The simulated Sox9 spots became partially fused into continuous regions, and those spots which did form were larger than in the control simulation (Fig. 5c). Consistent with this in silico prediction, treatment with a porcupine inhibitor C59, which inhibits Wnt secretion (thus Wnt production) 31 also resulted in a partial or complete fusion of Sox9 into a continuous domain parallel to the distal fin edge of S. canicula embryos (Fig. 5f). The later cartilaginous patterns of treated fins also showed continuous or larger condensations (bracket and arrowheads in Fig. 5i)-again reflecting the simulation result. We thus found a high consistency between the model predictions and the phenotypes of in vivo experimental perturbations, and also a remarkable similarity in response to these inhibitors between catshark and mouse. Taken together, this suggests that the distal elements of S. canicula pectoral fins and mouse digits share a deeply conserved Turing system.

Discussion
In this study, we have provided experimental and theoretical evidence that a Bmp-Sox9-Wnt Turing network represents a new example of deep homology-underlying skeletal patterning all the way from sharks to mammals. The molecular details are not identical (for example the Bmp4 ligand is the stronger candidate in the catshark, while it is BMP2 in the mouse), however the most striking feature of our results is that in both species the same basic interactions are seen between Bmp, Wnt, Sox9 and Fgf. Furthermore, the experimental perturbations of Bmp and Wnt signalling closely mirror both the model predictions and the results from mouse experiments 12 . In S. canicula fin buds, Sox9 forms a spot-like pattern, which is different to the stripe-like pattern in mouse digits. However, our computer simulation reveals that the BSW network of the mouse digit patterning 12 is also able to explain this spot pattern with just quantitative adjustments of its parameters. Therefore, relatively minor changes to the underlying deeply conserved network may be enough to trigger dramatic changes in skeletal arrangement.
Interestingly, teleosts seem not to use the BSW network (in zebrafish sox9a and b are expressed uniformly across the fin bud with no periodic pattern 32 , and bmp2a expression overlaps with sox9s (ref. 33)), and do not pattern their radials in the same manner (they produces a uniform endochondral disc, which is subsequently perforated to make the final skeletal pattern 34 ). Although convergent or parallel evolution is theoretically possible-with a Turing network developing separately in cartilaginous fish and tetrapods-the most parsimonious explanation is that the BSW network was lost (or significantly altered) in the teleost lineage.
Although most components of the mouse BSW model appear to have conserved roles in the catshark (Bmp, Wnt, Sox9 and Fgf), the exception to this is the distal Hox genes, which play no role in our catshark model. This decoupling of the BSW network with the distal Hox genes is suggested by our observations and previous studies 17, 26 that Hox expression domains do not overlap with the Sox9 expressing region significantly (Fig. 4b), and is consistent with recently published results highlighting the differences in Hox gene regulation between tetrapods and fish. In tetrapod limbs, Hoxa and Hoxd genes are regulated by two distinct genomic domains: a 3 0 domain regulating expression in the zeugopod and a 5 0 counterpart controlling autopod expression 3,31 . Fish fins, by contrast, do not have such strict relationship between Hox gene regulations and their anatomical regions. Although a bimodal regulation has been found in fish fins 3,6 , the expression pattern of Hoxa13b in zebrafish, for example, is almost uniform 35 . Similarly, misexpression of the distal Hox genes also causes very different results: misexpression of Hoxd13 or a13 in chick limb buds results in a truncation of zeugopod elements 36,37 , while a similar experiment in zebrafish pectoral fin buds causes an increase of cartilage condensation in the distal region 38 . Thus, much experimental data supports the idea that while the Hox gene regulation and anatomical modules are strictly coupled in tetrapod limb development, this modular regulation is less strict in fish fin development. Elaboration of Hox gene regulation, suggested by many other studies 3,10,38,39 , may be relevant to coupling the interaction between Hox genes and the BSW network in the digit patterning.
We have focused here only on the distal nodular bone formation in S. canicula-because it is the first periodic pattern to form in the fin bud-and thus the mechanism of the proximal stripe formation remains to be addressed. Because the distal nodular elements are each connected with a proximal stripe element in adult catshark fins, one possibility is that the proximal elements are formed just by elongation of the distal Sox9 expression spots. Our Wnt inhibition experiments question this idea, as even when the distal Sox9 expression becomes continuous, the stripe elements are still formed (though they are thicker and fewer in number), suggesting that formation of the stripe elements is not totally dependent on the patterning of the distal nodular elements (Fig. 5i). This semi-independent nature of the distal and proximal elements implies that additional unknown molecular controls might contribute to patterning of the proximal regions. Nevertheless, the very periodic nature of this pattern suggests that even if different molecules are involved, Sox9 is likely patterned by a Turing-type mechanism, which may be related to the BSW model.
Finally, we propose that the changes in the fin and limb skeletal arrangement may have involved a change in the role of Fgf gradient for organizing the Sox9 expression patterns. In S. canicula Fgf appears to act as a positional cue 40 -positioning the row of Sox9 spots at a certain distance from, and therefore parallel to, the distal fin edge-while in the mouse it appears to align the digital stripes perpendicular to the distal limb edge, and to control the wavelength 41 (Fig. 6). In the model, this dynamical difference can be explained by at least two parameters ( Supplementary Fig. 7). One is the ratio between Wnt and Bmp production terms, which can change spots to stripes-also demonstrated by the Wnt inhibitor treatment in Fig. 5. The other is the inhibition of Sox9 by Wnt (k 3 ), which affects the position of Turing space. Decreasing k 3 results in a shift of Turing space to the distal edge, allowing it to form a pattern in the distal domain. In our computer model, these Wnt-related parameters are enough to change the role of Fgf gradient from positioning spots to aligning stripes. Biologically, FGF and WNT are known to have a synergetic repressive activity to Sox9 expression 42 . Therefore, we could speculate that this synergetic activity of Wnt and Fgf might be relatively stronger in the distal mesenchyme of catshark fin bud than in the mouse digit forming region. Although the real fin-to-limb transformation must have involved more complex processes, including fin/limb shape changes, anterior-posterior patterning changes 43 , the loss of actinotrichia proteins 44 , Hox gene regulation 3,10,38,39 and others, our simple BSW model is nevertheless able to capture some key qualitative features of this morphogenetic change.
In conclusion, our study reveals that the morphological diversity of the distal fin and limb elements can be explained as the re-organization of a Turing patterning process. It highlights how relatively small regulatory changes can lead to major re-arrangements of the skeleton, and also emphasizes the difficulty of assigning homologous relationships between the distal elements of fins and limbs.

Methods
Animals. Experiments were performed in accordance with guidelines for animal experiments of Tokyo Tech and CRG, and experiments involving mice were approved by animal ethics committees of CRG (JMC-07-1001P3-JS). Catshark (Scyliorhinus canicula) eggs were provided by A. Tweedale (Bangor university) and Station Biologique de Roscoff, France, and they were incubated at 16°C in seawater and staged according to the standard staging system 45  paraformaldehyde in 1 M phosphate-buffered saline, dehydrated in a graded methanol series, and stored in 100% methanol at À 20°C.
Chemical treatments. S. canicula embryos at early stage 30 were removed from the egg shells and cultured in 6-well plate with 2B4 ml artificial seawater containing penicillin/streptomycin (Gibco). SU5402 (Sigma), LDN-193189 (Stemgent) and C59 (Merck Millipore) were dissolved in DMSO as stock solutions. The embryos were treated with 100 mM SU5402, 50 mM LDN-193189, 20 mM C59, or indicated concentrations and 1% DMSO during 4 days and fixed for in situ hybridization analysis. We also analysed gene expressions 2 days after the treatments, but found no significant differences ( Supplementary Fig. 6e). For qPCR analysis, the left and right pectoral fin buds of each inhibitor-treated embryos were dissected and pooled. RNA extraction was carried out with RNeasy Micro Kit (Qiagen), and cDNA synthesis was done with SuperScript III (Invitrogen). LightCycler 480 (Roche) and SYBR Green I (Roche) were used for the measurement of each gene expression amount. The qPCR primers were listed in Supplementary Table 1 x is any gene, DCT x i is a differential threshold cycle of gene x in sample i, and E x is PCR efficiency of gene x. For cartilage staining, embryos were cultured with the chemicals for 20 days and additional 10-20 days with the normal artificial seawater. Before alcian blue staining, embryos were permealized by xylene and staining was carried out with a standard protocol.
Fate map analysis. Eggs of S. canicula embryos around stages 26-28 were partially opened to label them with Indian ink (Pelican). To stop embryos from moving, eggs were cooled with iced seawater. Indian ink was injected into the pectoral fin buds with glass capillary, and the pictures were taken if applicable. After labelling the embryos, eggs were closed with plastic wraps and glue, and incubated in the artificial seawater containing penicillin/streptomycin (Gibco) for around 30 days. The labelled pectoral fin buds were scanned with OPT as described above. The images were manipulated with GIMP (the GNU Image Manipulation Program; http://www.gimp.org/).
In silico modelling. The technique for mouse limb modelling, which was implemented with Java 22 , was applied to the 2D fin growth model. Serial pictures of S. canicula pectoral fin buds were taken from stage 25 to stage 32 embryos ( Supplementary Fig. 3a), and their outlines were converted into spline curves. Then the shapes between each key stage were calculated by interpolation with one day temporal resolution (Supplementary Fig. 3b). Each fin shape had an independent triangle mesh implemented by gmsh 54 (Supplementary Fig. 3c). To carry out numerical simulations in the growing fin model, each mesh transmits information of species' concentrations to the next mesh. When a mesh was deformed to match the next shape, concentration of a triangle in the deformed mesh is split into the overlapping triangles of the next mesh (as previously described 22 ). This implementation also allows virtual fate map analysis (Supplementary Fig. 3d). We are happy to supply the code on request.
Mathematical analysis and numerical simulation. Mathematical analysis of Turing space and numerical simulations were carried out in conditions previously described 12,55 , in particular using the linear stability analysis described in White and Gilligan 55 . We considered the following general reaction-diffusion equations for Bmp (B), Sox9 (S) and Wnt (W): where D S , D B and D W are diffusion constants of S, B and W respectively.
To linearize them about the steady state (S*, B*, W*), we set where S 0 , B 0 and W 0 are constants, k is the wavenumber and s can either be a real number or a complex number. For |w| small, the equations (5)- (7) becomes where A is the Jacobian matrix at the steady state, and f, g and h are the partial derivatives of the indicated variables. For nontrivial solutions, the s is determined by the roots of the characteristic polynomial sI À A þ Dk 2 ¼ 0: Turing instability requires Re s k 2 ¼ 0 ð Þ ; Re s k 2 ð Þ40 for some k 2 : We implemented a simple linear model: where k 1 to k 9 are kinetic parameters representing regulatory interactions between genes. As described previously 12 , under D S ¼ 0, k 1 ¼ k 6 ¼ k 8 ¼ 0, k 2 40, k 3 o0, k 5 ¼ k 9 o0, we obtained the following inequality by solving (10) and (11) with Mathematica (Wolfram): k 4 o0; k 7 o0; k 3 o0; k 2 40; k 2 4 À k 3 k 7 k 4 ; D B 4 À D W k 2 k 4 k 3 k 7 : This inequality was used for determining the Turing space in Fig. 4c (parameter values were b ¼ 1, For numerical simulations, partial differential equations (PDEs) were solved by PDE solver written in Java with Huen method. Time step was 0.002. The finite volume method was used to calculate the amount of diffusion between neighbouring triangles. Zero-flux boundary condition was used in all simulations. Initial conditions were set as homogeneous steady states of each species. 1% of Gaussian multiplicative noise was added at each time step. Simulations were carried out from stage 29 to stage 31. Equations (1)-(3) were used in Supplementary Fig. 5a. Parameter values were b ¼ 1, k 2 ¼ 1, k 3 ¼ 1, k 4 ¼ 1, @W @t ¼ b a W À k 7 k F F S À k 9 W þ D W r 2 W À Á ð18Þ where F is the Fgf gradient defined below and k F is a constant. Parameter values where the decay rate, m F ¼ 0.1 and the diffusion constant, D F ¼ 600. a F is a local production from the fin edge where apical ectodermal ridge is formed in S. canicula fin buds. For the simulations on squares in Supplementary Fig. 7, the Fgf gradient was substituted to e À 3x (0rxr1).