Segmentation of the zebrafish axial skeleton relies on notochord sheath cells and not on the segmentation clock

Segmentation of the axial skeleton in amniotes depends on the segmentation clock, which patterns the paraxial mesoderm and the sclerotome. While the segmentation clock clearly operates in teleosts, the role of the sclerotome in establishing the axial skeleton is unclear. We severely disrupt zebrafish paraxial segmentation, yet observe a largely normal segmentation process of the chordacentra. We demonstrate that axial entpd5+ notochord sheath cells are responsible for chordacentrum mineralization, and serve as a marker for axial segmentation. While autonomous within the notochord sheath, entpd5 expression and centrum formation show some plasticity and can respond to myotome pattern. These observations reveal for the first time the dynamics of notochord segmentation in a teleost, and are consistent with an autonomous patterning mechanism that is influenced, but not determined by adjacent paraxial mesoderm. This behavior is not consistent with a clock-type mechanism in the notochord.


Introduction
The segmented vertebral column is the hallmark of vertebrate species. In many vertebrates, each ossified metameric unit starts out as a chordacentrum, a small, mineralized, ring-like structure around the notochord. The chordacentrum is then expanded and later is comprised of a barrel-like centrum (or vertebral body) (Figure 2-figure supplement 1) and protruding neural and hemal arches that enclose and protect the spinal cord throughout the body axis as well as the inner organs. In amniotes, such as mammals and birds, all these elements as well as adjacent muscle and sub-cutis derive from transient embryonic segmented structures termed somites. The sclerotome occupies the bulk of the somite, whereas skin and muscle are derived from the smaller derma-myotome. The cells within the sclerotome responsible for producing and mineralizing the organic aspect of bone (osteoid) are termed osteoblasts. In all vertebrates, somites are formed rhythmically and sequentially along the embryonic axis in a process governed by a multicellular, genetic oscillator termed the segmentation clock. Sequential waves of oscillating gene expression move through the precursor cells of the pre-somitic mesoderm (PSM) from the posterior to the anterior, where their arrest prefigures the timing and location of each newly forming somite boundary (Masamizu et al., 2006;Aulehla et al., 2008;Soroldoni et al., 2014;Shimojo and Kageyama, 2016). The segmentation clock contains a network of oscillating genes that shows differences between species, but members of the Hes/her family of transcriptional repressors appear to form a core negative feedback loop in all species examined (Krol et al., 2011). Mutations that disrupt the mouse clock (e.g. Hes7 [Bessho et al., 2001]) cause defective somitogenesis and corresponding defects to the centra and arches of the spinal column. Similarly, surgical perturbation of somites in chick embryos produce defects in the centra and arches (Stern and Keynes, 1987;Aoyama and Asamoto, 1988). Engineered changes to the mouse clock that shorten the period of oscillations correspondingly produce more vertebral bodies (Harima et al., 2013). Thus, the current model in amniotes proposes that the segmentation clock is the major source of patterning information for the somites, which in turn provides the segmented organization of sclerotomal derivatives that form the individual centra and arches of the spinal column. Remarkably, although somitogenesis is similar in vertebrate classes, whether the cellular lineage(s) and mechanism of segmental patterning of the chordacentra are homologous remains unclear (Fleming et al., 2015).
It has been proposed that in the teleosts zebrafish (Fleming et al., 2004) and salmon (Grotmol et al., 2003;Wang et al., 2013) the notochord and not the sclerotome is the initial source of bone matrix for chordacentra formation, while in other teleost species the classical sclerotomederived osteoblasts have been suggested as the main drivers of chordacentrum formation (Inohaya et al., 2007;Renn et al., 2013). Chordacentra are first observed as they mineralize as rings around the notochord, forming sequentially along the axis in a process that begins several days after somitogenesis and muscle segmentation is completed and lasts several weeks. This difference in developmental timing between the formation of segmented muscle and segmented skeleton in teleosts raises a registration problem between these mechanical elements, which amniotes do not need to solve. The zebrafish fused somites/tbx6 (fss) mutant, which has unsegmented paraxial mesoderm, has ectopically positioned neural and hemal arches, but was reported to form normal vertebral centra (Fleming et al., 2004;van Eeden et al., 1996). This suggests that a segmented sclerotome is necessary for proper arch development, but not for segmentation of the chordacentra. However, the zebrafish hes6 mutant, which forms somites more slowly than wild type, makes correspondingly fewer centra (Schrö ter and Oates, 2010), suggesting that the segmentation clock can influence vertebral patterning. Indeed, gene expression studies suggest that the tbx6 mutant retains a dynamic segmentation clock in the posterior of the PSM, but specifically fails to output the clock's information in the anterior PSM (van Eeden et al., 1998;Holley et al., 2000;Nikaido et al., 2002). The tbx6 À/À and hes6 À/À phenotypes might be reconciled if, even in the absence of morphological somitogenesis, an active clock in the posterior PSM of the tbx6 mutant is sufficient to correctly pattern the chordacentra. Alternatively, the segmentation clock may not be required for segmental chordacentra formation, but may nevertheless be able to indirectly modulate the pattern.
In this paper, we investigate the developmental origin of segmental patterning of chordacentra in zebrafish, showing that the majority of chordacentra can form normally despite a disrupted segmentation clock. We image the segmentation dynamics of the notochord sheath cell layer directly and document a distinctive set of segmentation errors in mutants. To unify these observations, we describe a model for autonomous segmentation of the notochord and validate this by comparison of simulations to experimental data.

Disruption of the segmentation clock in double and triple mutants
The zebrafish segmentation clock's core pacemaker circuit consists of her1, her7 and hes6, which display partial and overlapping redundancy (Schrö ter et al., 2012). The analysis of mutant combinations was previously not possible because her1 and her7 are located~10 Kb apart on chromosome 5. We generated a her1;her7 double mutant by injecting a TALEN construct directed against her1 in the her7 mutant (Choorapoikayil et al., 2012) and a novel hes6 mutation also using a TALEN approach (Figure 1-figure supplement 1). The her1 TALEN allele has two consecutive premature stop codons (TAA TAA). The predicted Her1 protein from the TALEN-induced mutation lacks the basic DNA-binding domain and the HLH dimerization domain, suggesting that the protein product has no functionality. It has been previously shown that the hu2526 allele is a her7 null mutant resulting from a stop codon in the HLH domain (Schrö ter et al., 2012). Importantly, the phenotype of her1;her7 double mutants is consistent with the phenotype of her1, her7 double morphants and also of the b567 deletion allele (Henry et al., 2002;Oates and Ho, 2002).
To assess patterning in the PSM, we used her7 ( Figure 1A-D'), her1 and deltaC expression (Figure 1-figure supplement 2), which show wave-like expression domains, as markers for the oscillation of the clock, and mespb, which shows expression in cells along the anterior border of two newly forming segments, for the clock's segmental output in the anterior PSM at the 10-somite stage ( Figure 1E-H; Figure 1-figure supplement 3A). In tbx6 À/À , her7 still oscillated posteriorly ( Figure 1B,B´), but mespb was not expressed ( Figure 1F), as expected (Oates et al., 2005;Durbin et al., 2000). In her1;her7 double mutants, her7 was expressed throughout the PSM ( Figure 1C,C'), but lacked oscillatory waves, and mespb expression in the anterior PSM occurred in a diffuse domain, lacking segmented stripes ( Figure 1G; Figure 1-figure supplement 3B), consistent with previous results using anti-sense knock-down reagents (Oates and Ho, 2002). Expression of other markers of segmental patterning in the anterior PSM (papc, ripply1, ripply2) also lacked segmental stripes (Figure 1-figure supplement 3C-H). This reveals that in her1 À/À ;her7 À/À , the segmentation clock is severely disrupted throughout the PSM, and that a correspondingly disordered, non-segmental output is made in the anterior PSM. Triple mutants for her1, her7 and hes6 have the same expression patterns as her1;her7 double mutants for all markers analyzed (Figure 1-figure supplement 2). Hence, and in order to simplify breeding, we focused further analyses on her1 À/À ; her7 À/À . Finally, we examined her1 À/À ;her7 À/À ;tbx6 À/À and found that her7 expression in the posterior PSM lacked oscillatory waves and neither her7 ( Figure 1D,D') nor mespb ( Figure 1H) was expressed in the anterior PSM. From these results we conclude that her1 À/À ;her7 À/À ;tbx6 À/À resembles the simple addition of the tbx6 and her1;her7 double mutant phenotypes; namely a severely disrupted segmentation clock that lacks any output in the anterior. In addition, bright field images at the 18-somite stage ( Figure 1I-L') were taken. Somite boundaries in wild type embryos are periodic and sharp ( Figure 1I,I'). Somite boundaries could not be discerned in tbx6 or her1;her7;tbx6 triple mutants ( Figure 1J,J' and L,L'). In her1 À/À ;her7 À/À ( Figure 1K,K'), partial boundaries were visible, but were lacking regular shape and periodic arrangement. This indicates that periodic morphology is disrupted in the mutant paraxial mesoderm. We next evaluated the presence of periodic patterns in muscle pioneers by analyzing en2a expression at the 20-somite stage ( Figure 1M-P'). The normal periodic pattern is lost in tbx6 À/À ( Figure 1N,N'), her1 À/À ;her7 À/À ( Figure 1O,O') and her1 À/À ; her7 À/À ;tbx6 À/À ( Figure 1P,P'), indicating that although these cell types are present in these mutants, there is no overt segmental pattern emerging directly from the PSM.
Normal centra form in the absence of periodic paraxial patterning After having established that her1;her7 double mutants, tbx6 single mutants and her1;her7;tbx6 triple mutants display severe disruption of the segmentation clock and its output in the paraxial mesoderm, we examined to what extent these early paraxial defects were reflected in vertebral bodies of the adult. If the previously reported ability of the fss/tbx6 mutant to form normal centra was due to the remaining segmentation clock activity in the posterior PSM, then we expected to see a strong disruption of vertebral bodies in both her1 À/À ;her7 À/À and her1 À/À ;her7 À/À ;tbx6 À/À adults. Mineralized bone was visualized using Alizarin Red (AR) staining of adult skeletons. We observed  Figure 1. Disruption of the segmentation clock in tbx6, her1;her7 and her1;her7;tbx6 mutants. (A-D') In situ hybridization for segmentation clock marker her7. (B and B') her7 oscillates in the posterior PSM of tbx6 À/À , but does not oscillate in her1 À/À ;her7 À/À (C and C´) or her1 À/À ;her7 À/À ;tbx6 À/À (D and D'). (E-H) In situ hybridization for segmental output marker mespb. mespb is not expressed in tbx6 À/À (F) or her1 À/À ;her7 À/À ;tbx6 À/À (H), but is weakly expressed in her1 À/À ;her7 À/À , albeit not in segmental stripes (G). (I-L') Somite boundaries in the paraxial mesoderm. In tbx6 À/À (J), her1 À/À ;her7 À/À (K) and her1 À/À ;her7 À/À ;tbx6 À/À (L) mutants, boundaries lose periodic order. (M-P') Spatial distribution of muscle pioneers marked by in situ hybridization with en2a. In tbx6 À/À (N), her1 À/À ;her7 À/À (O) and her1 À/À ;her7 À/À ;tbx6 À/À (P) muscle pioneers lose segmental pattern. A-H' are dorsal views of 13.5 hpf (10 somites) embryos, I-P' are lateral views of 18-19.5 hpf (18-20 somites) embryos. a -anterior, p -posterior. Scale bar in A is 100 mm and applies to A-G. Scale bar in I is 150 mm, applies to I-L and in I' is 100 mm, applies to I'-L'. Scale bars in M and M' are 150 mm and 100 mm respectively, and apply to M-P and M'-P' respectively. duplications, fusions and abnormalities of the neural and hemal arches throughout the axis in every tbx6 À/À , her1 À/À ;her7 À/À and her1 À/À ;her7 À/À ;tbx6 À/À adult, consistent with a loss of pattern in the sclerotome. However, the majority of centra in her1 À/À ;her7 À/À and her1 À/À ;her7 À/À ;tbx6 À/À adults were remarkably well-formed, defined as being cleanly separated from neighboring centra and similar to wildtype in their basic hourglass shape ( Figure 2F-H). Furthermore, we confirmed that all centra were well-formed in all her1, her7 and hes6 hetero-and homozygotes (Hanisch et al., 2013), in all heterozygous double and triple crosses, and in hes6;her7 double mutants. In addition, centra were also well-formed in mutants where the segmentation clock's oscillating cells slowly desynchronize due to a loss in Delta-Notch signaling, such as in beamter/deltaC and after eight/deltaD (Durbin et al., 2000) ( However, in some mutant combinations we did observe localized defects, such as a bend in the axis or occasional malformations and fusion of neighboring centra, scattered along every her1;her7 and her1;hes6 double mutant, and her1;her7;hes6 and her1;her7;tbx6 triple mutant skeletons, as well as in 80% of tbx6 À/À skeletons ( Figure 2F-H and Figure 2-figure supplement 4R,U). This shows that previous reports of normal segmentation of the centra in tbx6 À/À were incomplete (Fleming et al., 2004;van Eeden et al., 1996), and indicates that in the absence of periodic order in the early paraxial mesoderm, formation of the centra is error-prone. Given the proximity of the developing chordacentra to the notochord, and the suggestion that the notochord serves as a linear template for the vertebral column (Gray et al., 2014), our findings argue instead that formation of periodic chordacentra may arise from a separate segmentation mechanism intrinsic to the notochord.
Segmental entpd5 expression in notochord sheath cells is the key step for chordacentrum mineralization Secreted Entpd5 has previously been shown to be required for bone formation in zebrafish, and to be co-expressed with osterix (osx) in craniofacial osteoblasts (Huitema et al., 2012). When examining entpd5 promoter activity outside the craniofacial area, we observed a striking segmented pattern in the sheath cells along the notochord ( Figure 3A), well before the onset of chordacentrum mineralization is first observed in the anterior notochord at 6 dpf (Morin-Kensicki et al., 2002). Photoconversion at 3 dpf ( Figure 3B) of Kaede expressed from an entpd5:kaede transgene showed that these rings arise from an earlier ubiquitous expression domain by de novo synthesis ( Figure 3C), making entpd5 the earliest axial segmented marker known for zebrafish. We found entpd5 to be segmentally expressed first in the dorsal region of the sheath cell layer at the anterior, and new, periodically forming rings of entpd5 expression formed sequentially from anterior to posterior along the axis ( Figure 3C and D). entpd5+ sheath cells precisely predicted the position of the mineralized chordacentra and were located centrally to the mineralized matrix as revealed by co-staining with AR ( Figure 3E and F). In classical osteoblasts (cells able to produce bone matrix), as are found in the cleithrum or parasphenoid bones of the head, entpd5 is co-expressed with the osteoblast regulator osterix/Sp7 (Figure 3-figure supplement 1A) whereas within the developing vertebral column, osterix expression is first observed after 17 dpf when neural and hemal arches begin to form (Spoorendonk et al., 2008). entpd5 mutants form normally segmented osteoid, but do not     In wild type (n = 16), entpd5+ segments are added in an orderly manner from anterior to posterior (black lines). (B) tbx6 À/À (n = 4), (C) her1 À/À ;her7 À/À (n = 4), and (D) her1 À/À ;her7 À/À ;tbx6 À/À mutants (n = 4) also form entpd5+ rings in an anterior to posterior manner, but in an error-prone fashion. Defects such as gaps in the segmental pattern (increased intervertebral spaces), or small vertebrae followed by fusions (black boxes) can be seen scattered along the axis. Red asterisks (*) in her1 À/À ; her7 À/À kymogram ( The patterning mechanism in the notochord sheath is influenced, but not determined, by paraxial segmentation We reasoned that a better understanding of the occasional defects observed in the tbx6 À/À ,her1 À/À ; her7 À/À , and her1 À/À ;her7 À/À ;tbx6 À/À chordacentra might provide insight into the mechanism of chordacentra segmentation. To describe the dynamics of these defects we first imaged entpd5 expression in the sheath cells along the axis at intervals of 2 days in each mutant and recorded the distribution in a kymogram ( Figure  We next examined the dynamics of entpd5 expression rings in different mutant backgrounds. In the case of tbx6 (n = 4), her1;her7 (n = 4), her1;her7;tbx6 (n = 4) and hes6 (n = 6) mutants ( Figure 4B-D, Figure 5A  , the chordacentra formed in an anterior to posterior direction as in . However, we observed two types of defects: in the first, a ring was initially not formed, leaving a transient gap ( Figure 5A and longer spaces in the kymograms), which was then modified by the subsequent intercalation of a new entpd5 expression ring two or more days later (segment 1' in Figure 5A, 19 dpf). The intercalated ring was thinner than neighboring elements, likely reflecting an earlier stage in ring development, and remained smaller and distinct from neighboring segments. In the second type of defect, a ring was added (on schedule or out of schedule) in the middle of a normal intervertebral distance, creating a distance shorter than expected between two rings. In this case, the notochord sheath cells of the small segment fused with one or two of its neighboring segments ( Figure 5B, black boxes in kymograms). In the cases where two sequential defects are inserted by any of the two processes described above, we observed a transient bending of the axis, which was later modified by fusing segments or leaving two or three contiguous small vertebrae ( Figure 4C, her1 À/À ;her7 À/À kymogram and Figure 4-figure supplement 3). The phenotypic defects we see are different to the scoliosis phenotypes reported in mutant leviathan/col8a1a animals with defects in collagen deposition (Gray et al., 2014) or in mutant stocksteif/cyp26b1 animals with increased retinoic acid levels (Spoorendonk et al., 2008). In these cases, the first pattern of mineralization appears correctly segmented, and defects arise because of axial bending of the notochord and bone overgrowth, respectively. In contrast, we have directly observed the initial emergence of defective entpd5+ rings, suggesting that the periodic patterning of the notochord has been affected.
Segment defects in the paraxial mesoderm occur sequentially along the axis, consistent with a disrupted clock-type mechanism. In contrast, the intercalation defects observed in the mutant axes  form out of schedule with the sequence of notochord segmentation. We asked if the intercalation of an entpd5 ring was associated with an error in the initial local spacing of the rings along the notochord. We measured the distance between the entpd5+ rings bordering the location where the intercalated segment was added at the time point immediately before its appearance in her1;her7 mutants (n = 4 animals; 18 intercalations, red squares; Figure 5C) and compared this to equivalent distances in control embryos (n = 16 animals, 288 segments, black crosses). The distance preceding an intercalation was either similar to or larger than expected from the controls, but never smaller (n = 18, 1.4 ± 0.3 fold increase, mean ± SD). Thus, the intercalatory rings formed at a range of distances from their earlier neighbors, suggesting that the defects arise as a response to an error in positioning the earlier rings and are not a simple delay in entpd5 expression. These distinctive intercalations are difficult to reconcile with a clock-type mechanism for segmenting the notochord.
The chordacentra defects described above may reflect some influence of the disrupted paraxial mesoderm pattern on the notochord. Previously, it was shown that hes6 mutants show weak myotome boundary defects in the posterior tail at low penetrance at 1.5 dpf (Schrö ter and Oates, 2010), and the hes6 mutant generated in this paper also show these defects (9/67 (13.4%)). Nevertheless, the existence and distribution of skeletal defects had not been examined. In hes6 À/À skeletons we observed from one to four centra defects, restricted to the caudal vertebrae, in AR stains of adult bone (4/15, 27%) ( Figure 6A and B) and in entpd5 expression at 28 dpf (8/41, 20%) ( Figure 6C-E). Additional analysis showed that 2/5 (40%) hes6 mutants with a chordacentra defect revealed by entpd5 expression also presented a myotome defect in the same region. The remaining 60% (3/5) had no myotome defect, but still showed a chordacentra defect. This partial overlap of myotome and chordacentra defects in hes6 mutants, as well as the coordinated reduction in total myotome and vertebral number observed previously (Schrö ter and Oates, 2010) and in our work suggests that the pattern of the paraxial mesoderm can influence the segmentation of the notochord.
A reaction-diffusion model of axial patterning in the zebrafish These findings can be synthesized in a model of axial segmentation in which the notochord possesses an intrinsic segmentation mechanism, likely within the sheath cells, that does not depend on the paraxial segmentation clock to produce periodic entpd5 rings and subsequent mineralization. This mechanism is proposed to act directly upstream of entpd5 expression, but is nevertheless sensitive to information from the paraxial mesoderm, likely from the myotome structure, which can bias the position of a ring to enable coupling of the early-developing myotome with the later-forming skeleton in wildtype.
To assess the plausibility of this hypothesis and to investigate what kind of mechanism has these properties, we developed a theoretical description that formalizes these ideas and incorporates key experimental findings. We describe the intrinsic patterning mechanism operating in the notochord sheath cells as a reaction diffusion system with two components, an activator and an inhibitor (Murray, 1993) (see the Theory section in Materials and methods and Figure 7-figure supplement 1). This theory is capable of producing an autonomous pattern, sequentially adding segments from anterior to posterior ( Figure 7A, Figure 7-figure supplement 2 and Video 1). The cues provided by the paraxial mesoderm pattern are introduced as a distribution of sinks for the inhibitor that are of the same order as the mechanism's intrinsic wavelength (Figure 7). The choice of sinks instead of sources, together with vanishing initial conditions across the notochord except for a perturbation localized at the anterior, are required to preserve the sequential character of notochord segmentation. Although other more complex descriptions that improve robustness to noise are possible (Materials and methods), this simple theory successfully accounts for the sequential formation of regular segments in the presence of sinks as observed in wild type ( Figure 7B, Figure 7-figure supplement 3 and Video 2). The intrinsic patterning mechanism allows for a range of pattern wavelengths, Figure 5 continued (red dots) and compared to the equivalent axial position in wild type (WT, n = 16) (black crosses). The distance preceding an intercalation in her1;her7 mutants was either similar or larger than wild type. DM, distance measured; E entpd5+ segment. All scale bars are 100 mm. DOI: https://doi.org/10.7554/eLife.33843.019 providing the plasticity for the sinks to pin segments to specific locations, altering the length of each segment. We conjecture that the mutants do not affect the intrinsic notochord patterning mechanism, but change the features of the sink distribution. The potential of sinks for biasing the pattern can interfere with ring formation, revealing a process that can intercalate a ring into a mispatterned gap in the sequence. Both reducing the sink strength and increasing the noise in the positioning of sinks can induce defects in the intrinsic patterning mechanism (Figure 7-figure supplement 4). Thus, the spatially disordered and variable sized myotomes observed in the mutants (Figure 2A-E) are described in the theory as large fluctuations in the positions and amplitudes of inhibitor sinks ( Figure 7C and D). The her1 À/À ;her7 À/À situation is described by strong sinks with large position errors ( Figures 2C and 7B, Figure 7-figure supplement 5 and Video 3) while tbx6 À/À and her1 À/À ;her7 À/À ;tbx6 À/À are described by sinks with reduced amplitude and shorter wavelength, accounting for the smaller and scattered myotome fragments ( Figures 2B,D and and 7D, Figure 7-figure supplement 6 and Video 4). In the framework of the theory, tbx6 À/À and her1 À/À ;her7 À/À ;tbx6 À/À are described by the same set of parameters. The mechanism has some flexibility and is compatible with a range of wavelengths. For example, the larger myotomes  Figure 7. A reaction diffusion theory accounts for key experimental findings. A sink profile (blue) describes cues from myotomes that bias the position of segments. The Entpd5 pattern is given by the concentration of an activator (green) that is regulated by an inhibitor (red). (A) The system is capable of autonomous pattern formation in the absence of sinks. (B) Wild type condition is described by regularly placed strong sinks, according to the output of a functioning segmentation clock. (C) In her1 À/À ;her7 À/À strong sinks are misplaced due to a malfunctioning segmentation clock, causing segments to be also misplaced and giving rise to defects. (D) tbx6 and her1;her7;tbx6 mutants are characterized by weaker segmentation clock output and fragmented and scattered myotome boundaries, here described by weaker sinks with a shorter wavelength. (E) The hes6 mutant is here characterized by a sink profile wavelength that is 6% larger than wild type. Parameters: a = 10 À3 , b = 10 À2 , t= 0.   In the simulations of the wild type, there is always a correspondence between the position of the sink and the position of the peak of activator ( Figure 7B and E); this situation is also found between the positions of the myotome boundary and the entpd5 expression ring in experimental wild type animals (Figure 7-figure supplement 8A-F). In the simulations of tbx6, her1;her7 and her1;her7; tbx6 mutants this strict correspondence between sink and activator is lost ( Figure 7C and D); activator peaks occur both together with sinks and in between them. In the case of hes6, the strict correspondence between sink and activator may be lost in the last tail segments. To test this prediction of the model, we examined the distribution of myotome boundaries and entpd5 expression rings in tbx6 and her1;her7 double mutants. We observed that the disorganized myotome boundary fragments in the mutants had lost strict correspondence with the entpd5 rings (Figure 7-figure supplement 8G-R). This lack of spatial correspondence between paraxial structure and axial expression is in agreement with the model, and it supports the hypothesis that the sink is associated with some feature of the myotome boundary.
Lastly, in order to compare the quantitative output of the model against the biological data, we counted the number of vertebral segments as well as the number of segmentation defects (small vertebrae or fusions) in mutant fish using entpd5 expression at 28 dpf, when the vertebral bodies have developed into a mature form (Figure 7-figure supplement 9A,B). All mutants had a greater variability in segment number than wild type (n = 48), where segment number ranged from 29 to 31. In her1 À/À ;her7 À/À (n = 38) the segment number increased up to 37. The increase in her1;her7 was slightly ameliorated by the additional removal of tbx6 in the triple mutants (n = 47). As expected, hes6 mutants (n = 41) had a lower number of segments than the wild type, ranging from 28 to 31 segments. When the same fish were scored for segmentation defects (Figure 7-figure supplement 9B), her1 À/À ;her7 À/À presented a higher frequency and a broader distribution of defects than tbx6 mutants (n = 46). In her1;her7;tbx6 mutants, the number of defects per embryo decreased almost to the distribution seen in tbx6 mutants. The trends in these measurements are captured by the theory, as shown by the number of peaks in the pattern and (

Discussion
In order to dissect the interrelationship of paraxial mesoderm segmentation and chordacentrum formation in teleosts, we here use the zebrafish to manipulate the segmentation clock in the paraxial mesoderm, while using a novel marker for axial segmentation to examine chordacentrum and vertebral body patterning at early and late stages.
The zebrafish notochord consists of a core of vacuolated cells responsible for turgor pressure, and an outer, squamous epithelial layer, the sheath cells, that is adjacent to the notochord's basal lamina or sheath, which provides mechanical stability (Melby et al., 1996;Yamamoto et al., 2010) ( Figure 3E). Previously, we described the enzyme ectonucleoside triphosphate/diphosphohydrolase 5 (Entpd5) as essential for ossification in zebrafish (Huitema et al., 2012). Entpd5 mutants fail to mineralize osteoid, leading to a complete absence of bone. entpd5 is co-expressed with osterix in all osteoblasts of the early head skeleton, the cleithrum and the operculum (Huitema et al., 2012), but we here show that entpd5 is also expressed in notochord cells (Figure 3), and that, at 4 dpf, the expression in the axial mesoderm becomes restricted to segmentally organized rings in the notochord sheath cells. osx is not expressed at these early stages in the notochord sheath cells, nor do osx mutant present with an axial mineralization phenotype in zebrafish or Medaka (Figure 3-figure supplement 1; Yu et al., 2017). Combined, these results identify the osx-, entpd5+ sheath cells as those osteoblasts responsible for initial chordacentra formation, and indicate that the segmental patterning of the notochord is under way by 4 dpf in these cells. Mechanistically, the segmental expression of Entpd5 in the notochord sheath provides the enzymatic activity to mineralize the fibrous sheath of the notochord. It is possible that the respective contribution of notochord sheath cells in chordacentrum formation versus sclerotome-derived osteoblast function in centrum formation is different in other Video 3. Simulation of the her1;her7 mutant condition. The sink profile (blue) showing the noisy spatial distribution for the her1;her7 mutant condition in the top panel and corresponding activator (green) and inhibitor (red) patterns in the bottom. Parameters as in Figure 7 of the main text. DOI: https://doi.org/10.7554/eLife.33843.023 Video 4. Simulation of the tbx6 mutant condition. The sink profile (blue) representing the noisy spatial distribution and the reduced amplitude of sinks in the tbx6 mutant in the top panel and corresponding activator (green) and inhibitor (red) patterns in the bottom. This simulation also represents the her1;her7; tbx6 mutant. Parameters as in Figure 7 of the main text. DOI: https://doi.org/10.7554/eLife.33843.024 Video 5. Simulation of the hes6 mutant condition. The sink profile (blue) representing the longer spatial wavelength of sinks in the hes6 mutant in the top panel and corresponding activator (green) and inhibitor (red) patterns in the bottom. Parameters as in Figure 7 of the main text. DOI: https://doi.org/10.7554/eLife.33843.025 teleost species (Fleming et al., 2015;Yasutake et al., 2004;Kaneko et al., 2016), and this can be tested by generating entpd5 mutants in some of those species. In zebrafish, our findings settle a long-standing debate on which cells are functionally important for chordacentrum formation.
From an evolutionary perspective, it has been suggested that the chordacentrum is a novelty found only in actinopterygians (Arratia et al., 2001). It seems safe to state that chordacentrum formation is conserved at least among teleosts, while the situation for holosteans (Amia and the gars) is less clear. Whether chondrosteans such as Acipenser, which lack centra but exhibit mineralized neural and hemal arches are different deserves further studies in order to bridge the gap between paleontological, anatomical and molecular studies (Arratia et al., 2001;Laerm, 1976;Laerm, 1979). Recent work has shown that clock-based mechanisms are involved in body segmentation in vertebrates as well as invertebrates (Sarrazin et al., 2012) and there is wide agreement that amniotes segmentally pattern their muscles and axial skeletons at the same time during development using the 'segmentation clock'. However, from fossil and extant taxa it seems likely that neural and hemal arches are evolutionarily older than vertebral bodies (Cote et al., 2002), and while these structures are homologous across vertebrate classes (Williams, 1959), vertebral bodies are not. The demonstration of periodic entpd5 expression in the notochord sheath cells now provides a tool to study the earliest visible patterning events in the axial skeleton, and also allows an evaluation of the segmentation clock's input. To perturb the clock, we generated new alleles and new allelic combinations in tbx6 and the central clock genes, her1, hes6, and her7. Analysis using a number of different readouts demonstrate that the phenotype of her1;her7 double mutants is not exacerbated in the further absence of hes6. Strongest disruption of the clock is seen in the case of her1;her7;tbx6 triple mutants. The mutants we have used to interfere with the segmentation clock and its output leave the PSM without any overt signs of oscillation or segmental pattern. It is of course formally possible that some remaining clock-like activity exists that we cannot detect. However, our argument in this paper does not require that we have removed all segmentation clock activity, but rather that our perturbations of the segmentation clock and its output show highly ordered axial structures despite strong disorder in the early paraxial mesoderm.
An influence of the paraxial mesoderm on axial segmentation may be reflected in the differences between the mutants' distributions of defects. The strongest axial phenotype is found in her1 À/À ; her7 À/À , which possesses prominent myotome boundary fragments in the paraxial mesoderm, whereas tbx6 À/À and her1 À/À ;her7 À/À ;tbx6 À/À have a slightly weaker axial phenotype, which is preceded by a greater reduction of myotome boundaries. This is consistent with a dominant interfering effect of paraxial structures. The similarities in paraxial and axial phenotypes of tbx6 À/À and her1 À/À ; her7 À/À ;tbx6 À/À is explained by proposing that the effects of a disrupted segmentation clock (that can nevertheless output a disrupted pattern seen in her1 À/À ;her7 À/À ) are reduced in her1 À/À ; her7 À/À ;tbx6 À/À by removing the output function of tbx6 in the anterior PSM.
We have attributed the occasional defects in the mutant centra to the effects of fluctuations in the strength and position of sinks for the inhibitor component of the notochord's intrinsic segmentation mechanism. In the absence of any sinks, the notochord would in our model segment without any defects with a periodicity determined entirely by the internal dynamics. Thus, the defects can be viewed as a 'dominant interference' by the mis-patterned paraxial mesoderm. Unfortunately, due to the random and unpredictable position of the chordacentra defects seen in all mutants analyzed, cell transplantation experiments cannot serve as a suitable experimental strategy to test this. We have also contemplated the possibility that her1, her7 and tbx6 are also expressed in the notochord and may be directly responsible for the centra defects seen in the mutants. However, this appears highly unlikely, because no tbx6 (Wanglar et al., 2014), mespb or her1 expression (In situ hybridization or transgene expression at 5 dpf and 13 dpf, data not shown) has been detected in the notochord at somite stages or during chordacentra formation. A recent report indicates that overexpression of mespbb in the notochord sheath can regulate segment size, possibly via interfering with boundary formation (Wopat et al., 2018).
Another interpretation is possible in which the notochord's intrinsic mechanism is highly noisy, and the paraxial mesoderm has no influence. While any molecular patterning mechanism would be vulnerable to fluctuations in the embryo, since the notochord's segments are generated over a long time scale, with a period about 20 times slower than somitogenesis, it seems plausible that the large molecular numbers that could be synthesized in this interval could ensure a low error rate. Furthermore, the transfer of information from a segmentation clock-derived muscle pre-pattern to an autonomous but plastic mechanism in the notochord resolves the apparent discrepancy between tbx6 and hes6 mutant phenotypes, and explains how the animal coordinates segmented muscles and axial skeleton in biomechanical register despite the fact that they develop days to weeks apart in time. Such a coordination mechanism may be a widespread feature of the teleosts and potentially other vertebrates.
In summary, we have proposed a second mechanism for periodic segmentation of the axial skeleton of a vertebrate species, which exists in the sheath cells of the zebrafish notochord. The existence of intercalary defects argues that although there is a sequential generation of periodic entpd5 expression rings, this periodicity does not arise from a clock-type segmentation mechanism as found in the paraxial mesoderm. We have described this dynamic mechanism as a reaction-diffusion process that sweeps down the axis, generating segments autonomously and refining their position using information from the previously segmented myotomes. The genetic basis of the proposed autonomous notochord mechanism and the signal from myotome to notochord are not understood, and our model does not assume any specifics with regard to molecular identities, but it opens the door to the search for these molecules. entpd5 is the earliest known molecular marker for the sites of ossification, with plasticity in the presence of perturbation, suggesting that understanding the control of entpd5 expression may hold the key. Sequence-based reagent Talen hes6 this paper See Figure 1,

Genotyping
DNA was isolated through fin clippings and from embryos. Genotyping was performed as described in Table 1 and Table 2.

In situ hybridization and alizarin Red bone staining
Riboprobes were generated from either plasmids or PCR templates and in situ hybridisation was performed as previously described (Oates and Ho, 2002). Whole mount stained embryos were documented on an Olympus SZX10 stereoscope with a QImaging Micropublisher camera. Flatmounted embryos were photographed on an Olympus MVX10 stereoscope with an Olympus DP22 camera. Alizarin Red bone staining was performed as described previously (Spoorendonk et al., 2008) for fish between 8 weeks and one year of age.

Virtual time lapses
Embryos were kept in E3 at 28˚C until 7 dpf. At 7 dpf, ten embryos with fully developed swim bladder from mutant (tbx6 À/À , her1 À/À ;her7 À/À , her1 À/À ;her7 À/À ;tbx6 À/À , or hes6 À/À ) or transgenic entpd5: kaede lines were anaesthetized, photographed and then housed individually in the animal facility. Individuals were fed tetrahymena in combination with Gemma 75 for the first two weeks, followed by artemia and Gemma 150 for the following two weeks. Every second day, each individual was, anesthetized (described above) and photographed using a Olympus SXZ16 stereomicroscope (1.5X PlanApo objective) connected to a DFC450C Leica camera. The embryo was placed in a drop of E3 with anesthetic on the lid of a petri dish. Each picture was taken from the cleithrum to the posterior tip of the individual. Immediately afterwards, the embryo was returned to warm E3 without anesthetic. Embryos were returned to their specific tank in the animal facility only when they were completely awake and moving. This procedure was repeated until all entpd5+ segments had developed in the axis. The sedation and photography did not take more than two minutes per embryo, and did not compromise survival.

Imaging
To photograph somite boundaries, 18-19.5 hpf embryos were dechorionated and laterally aligned in conical depressions that fit the yolk, in an agarose pad (Sigma, 2% in E3) cast in a petri dish (Falcon, 50 mm x 9 mm) and topped up with E3. Photomicrographs were taken on an Olympus MVX10 microscope equipped with an Olympus DP22 camera. For imaging the notochord, embryos were mounted in 0.5% low melting point agarose in a culture dish with a cover slip replacing the bottom. Fluorescent imaging was performed with a Leica SPE 'live' Confocal Microscope, Leica SP8 confocal microscope and a PerkinElmer Ultraview VoX spinning disk microscope using a 10x or 20x objective with digital zoom. Usually, z-stacks with intervals of approximately 2 mm were captured and were then flattened by maximum projection in ImageJ. For photoconversion of whole entpd5:kaede embryos, the green Kaede fluorophore was photoconverted using a Leica fluorescence microscope by 15-30 min exposure through a UV light bandpass filter (360/40 nm, 100 W mercury lamp). For bone stains an Olympus S2 Â 16 microscope coupled to a Leica DFC420c camera was used. Photomicrographs were stitched with the pairwise stitching plugin in Fiji (RRID:SCR_002285) (Preibisch et al., 2009).

Theory
The theory describes notochord sheath cells pattern formation in terms of a one dimensional reaction diffusion system with two components, an activator U and an inhibitor V, see Figure 7-figure supplement 1A. The concentration U = [U] of the activator is reflected in Entpd5 concentration. The concentrations of both the activator and inhibitor species U x; t ð Þ and V x; t ð Þ depend on position x and time t. We propose a variant of the FitzHugh-Nagumo (FHN) model (Murray, 1993) where D U and D V are diffusion coefficients for U and V respectively, and k i are rate constants. The choice of the FHN model is based on its simplicity. There are positive linear terms for the activator and negative linear terms for the inhibitor in both Equations (1) and (2), and a single nonlinearity, the cubic term for the activator that limits growth and allows for the stabilization of steady states. The model is meant to represent a plausible mechanism rather than specifying the interactions of particular molecules. We reduce the number of parameters by transforming this theory to a dimensionless form. We first set the source term k 0 ¼ 0 since as discussed below we need to reproduce a sequential pattern formation. We introduce a lengthscale L that we will take as the system size, and a timescale T and concentration scale U 0 to be set below. In terms of these scales we define new variables x 0 , t 0 , u and v such that and replace in the reaction diffusion equations above dropping the primes for notational convenience We multiply both equations by T=U 0 to render them dimensionless Multiplying the inhibitor equation by k 1 =k 5 and rearranging terms where we have highlighted dimensionless groups in parentheses. We now select a timescale by setting and a concentration scale by setting and define dimensionless parameter groups With these definitions We additionally set k 4 ¼ 1 for simplicity and we call k 6 ¼ d In the rest of this work we consider this dimensionless form of the theory. Here a and b are dimensionless scaled diffusion coefficients of the activator and inhibitor species respectively, t is a relative timescale, and d is a dimensionless degradation constant of the inhibitor.
We first consider the homogeneous system q xx u ¼ q xx v ¼ 0. The resulting equations for the local reactions are where dots denote time derivatives. Introducing functions the nullclines of the system, defined by setting f u; v ð Þ ¼ 0 and g u; v ð Þ ¼ 0, are the curves in the u; The first one is an inverted cubic that goes through the origin and the second one is a linear function with slope d À1 controlled by the single bifurcation parameter d, see Figure 7-figure supplement 1B. Intersections of these two curves are the solutions to and define the fixed points of the system where u There is always a solution u 0 ; v 0 ð Þ ¼ 0; 0 ð Þ and for d>1 there are two additional solutions u AE satisfying The linear stability of fixed points is determined by the matrix where and derivatives are evaluated at the fixed point u * ; v * ð Þ. The condition for stability is that and For the fixed point with determinant and trace Given that t; d>0 this implies that the origin 0; 0 ð Þ is a stable fixed point if d < 1 and t < d: We consider in the following a situation in which the fixed point 0; 0 ð Þ is stable, setting the dimensionless timescale t ¼ 0:1 and degradation d ¼ 0:5. Under appropriate conditions, the dimensionless reaction diffusion theory Equations (18) and (19) can give rise to pattern formation. In particular we require that the activator diffuses slower than the inhibitor, a<b, and here set a ¼ 10 À3 and b ¼ 10 À2 . Because of the signs of the derivatives near the fixed point, the type of pattern predicted is as displayed in Figure 7-figure supplement 1C, where there is coexistence of activator and inhibitor (Murray, 1993).
In the presence of random perturbations distributed along the notochord the homogeneous state loses stability due to differential diffusion. A pattern may form out of this initial random background fluctuation, with segments forming almost simultaneously all along the notochord. Although segment formation is robust in this scenario and can accommodate a broad range of wavelengths, this is at odds with the experimental observation that ENTPD5 segments form sequentially from anterior to posterior.
We conjectured that one way to obtain a sequential segment formation is to start with an initial perturbation localized at the anterior, and vanishing concentrations all across the rest of the notochord. Since the anterior of the vertebrate axis is always more developmentally advanced than the posterior, such an anterior perturbation is a plausible hypothesis. A vanishing concentration along the notochord is important to ensure that patterning is not triggered until the wave of activator and inhibitor arrives at a given point. Therefore, we start simulations with initial conditions that have zero concentration for both u and v across the whole domain x 2 0; L ð Þ, except for a small perturbation near the origin x ¼ 0. The form of the initial condition is a smooth step u x; 0 ð Þ ¼ u 0 2 1 À tanh 5:0 x À 0:5 ð Þ ð Þ ð Þ Starting from such a small perturbation in the anterior, we observe that the system is able to form a pattern sequentially, from anterior to posterior, see Figure 7A, Figure 7-figure supplement 2 and Video 1. Thus, the theory proposed is capable of autonomous patterning of the notochord in the absence of input from the segmentation clock.
We next introduce the effect of the segmentation clock input into this otherwise autonomous patterning system as a spatial dependent degradation profile of the inhibitor This sink profile s ¼ s x ð Þ for the inhibitor has peaks at given positions along the x axis, describing the cues that the notochord patterning mechanism receives from myotomes. At positions where s x ð Þ is large, the inhibitor is locally degraded at a larger rate. Note that there is no source term for the activator since we set k 0 ¼ 0 above. This feature together with the choice of sinks instead of sources to describe the segmentation clock cues are motivated from the observation that Entpd5 segments form sequentially. The presence of sources would render pattern formation non sequential. The sink profile s x ð Þ is characterized by a sink strength S 0 , a wavelength l and sink wavelength variability s. The first sink is positioned at l=2 and the positions X i of consecutive sinks are determined by the wavelength l with an error drawn from a uniform distribution of width s. The sink profile is built from a combination of tanh ::: ð Þ functions to produce smooth peaks of steepness a and width dS In this work we fix the values a ¼ 100 We consider a system size L that we set to L ¼ 17:1 so that the wildtype condition makes 30 segments with the sink profile natural wavelength l ¼ 0:57. We normalize axes length scales to this value in all plots. For simplicity we assume that the activator and inhibitor are restricted to notochord sheath cells and we specify Neumann boundary conditions, that is derivatives at both ends are zero qu qx j x¼0 ¼ qv qx j x¼L ¼ 0 : As described above, here we ensure the sequential character of the patterning through an initial perturbation at the anterior and vanishing concentrations across the notochord. One may query the robustness of such scenario, since noise across the notochord hampers sequential patterning. An alternative hypothesis would be to postulate an additional wavefront that propagates through the notochord progressively turning on the reaction diffusion mechanism of Equations (41) and (42) as it goes. To illustrate this we consider an alternative dimensionless form of Equations (1) and (2). Turning back to Equations (9) and (10) we select a timescale and concentration scale setting and k 3 U 2 0 k 1 1: Introducing dimensionless parameter groups and setting for simplicity k 4 ¼ k 5 ¼ 1 and k 6 ¼ k we arrive at the dimensionless form In this alternative dimensionless form it is straightforward to decouple the reactions from diffusion by tuning the value of g. Thus, we can introduce a wavefront g x; t ð Þ that moves from anterior to posterior turning on the reactions in its wake. Such a wavefront could have a biological origin in a molecular maturation gradient invading the notochord from the anterior. Due to very slow dynamics before wavefront arrival, this would render the patterning mechanism more robust to noise across the notochord. Yet a different possibility is a scenario of patterning in a growing domain (Crampin et al., 1999), although here the tissue where the pattern forms exists previous to the establishment of the pattern. These alternatives are certainly interesting and the underlying mechanism patterning the notochord will be the subject of future work. However, in the absence of experimental data supporting other hypotheses, here we settle on the perhaps more parsimonious choice of vanishing initial conditions across the notochord.
In this work we solve the partial differential equations described above using a custom python code, see Source code 1. We discretize space with a discretization length Dx ¼ 0:01, and time discretization is chosen as Dt ¼ 0:9Dx 2 =2. We integrate the partial differential equations until the pattern reaches a steady state.

Data availability
Virtual time lapse data and theory code can be found at: http://icor-data.uni-muenster.de/. Source data files and source code have been submitted to eLife as additional material.