Turing's turtles all the way down: A conserved role of EDAR in the carapacial ridge suggests a deep homology of prepatterns across ectodermal appendages

The scutes of the turtle shell are epidermal shields that begin their formation during the early stages of shell development. Like other skin appendages, turtle scutes are hypothesized to be patterned by reaction – diffusion systems. We have previously established ex vivo and in silico systems to study these mechanisms experimentally and have further shown that mathematical models can explain the dynamics of the induction of turtle scute primordia and the generation of final scute architecture. Using these foundations, we expand our current knowledge and test the roles of ectodysplasin and activin signaling in the development of turtle scutes. We find that these molecules play important roles in the prepatterning of scute primordia along the carapacial ridge and show that blocking Edar signaling may lead to a complete loss of marginal scute primordia. We show that it is possible to reproduce these observations using simple mathematical modeling, thereby suggesting a stabilizing role for ecto-dysplasin within the reaction – diffusion mechanisms. Finally, we argue that our findings further entrench turtle scutes within a class of developmental systems composed of hierarchically nested reaction – diffusion mechanisms, which is conserved across ectodermal organs


| INTRODUCTION 1.| Ectodermal appendages are a key vertebrate innovation
Morphological innovations have been linked to the success and adaptive radiation of animal clades during evolution.One of those structures that have been considered such novelties are ectodermal appendages (Wu et al., 2004), a class of placode-derived integumental organs common to all classes of vertebrates.While some types of ectodermal appendages, such as scales and teeth, are present across most extant and extinct clades, others, such as hair, feathers, horns, and mammary glands, tend to be rather clade-specific organs.A very distinct and conspicuous, yet hitherto understudied, type of ectodermal organs are turtle scutes.These are distinctively arranged structures that form the outer skin plates of carapace and plastron.So far, there are relatively few studies on turtle scute development, since their amenability to most experimental pipelines is, at best, limited (Moustakas-Verho et al., 2019).Yet, the availability of certain molecular and histological methods in conjunction with computational approaches has provided a suitable strategy to unravel partly the mechanisms underlying the development of these unique structures (Moustakas-Verho et al., 2014).Recent work has accumulated new data that allow to compare commonalities and differences between turtle scutes and other types of ectodermal appendages (Ascarrunz & S anchez-Villagra, 2022;Moustakas-Verho et al., 2014).
In general, ectodermal appendages emerge by specific and reciprocal crosstalk between committed epithelia and the underlying mesoderm (reviewed in Pispa & Thesleff, 2003, Biggs & Mikkola, 2014).This becomes evident in the formation of placodes, epithelial thickenings that arise by epithelial proliferation, while the mesenchyme begins condensing underneath.Finally, differentiation is induced in some central cells which cease dividing and begin secreting different kinds of extracellular matrix proteins that, finally, will give rise to the plethora of different types of appendage structures.These tissue-level changes are, again, orchestrated by mutual signaling interactions, involving for example, the canonical Wnt pathway, fibroblast growth factors (Fgfs), bone morphogenetic proteins (Bmps), hedgehog (Shh), and ectodysplasin (Eda) and its receptor (Edar) (reviewed in Widelitz & Chuong, 1999, Pispa & Thesleff, 2003, Biggs & Mikkola, 2014, Mogoll on et al., 2021).These molecular pathways are commonly accepted to generate, by their interaction, the biomechanical and cell density changes that give rise to Turing dynamics.Turing dynamics are local symmetry breaks (instabilities) that lead to the emergence of a spatially alternating pattern with a specific distance between neighboring maxima and minima.These then explain the emergence of periodic patterns (Glover et al., 2017;Inaba et al., 2019;Jernvall & Thesleff, 2000;Meinhardt & Gierer, 1974), both in the spacing of single ectodermal organs and in their substructures.In later developmental stages, tissue differentiation processes and secretion of specialized extracellular matrix proteins will account for, and amplify, morphological differences between ectodermal organs.During these processes, interactions of cell behaviors and biomechanical processes may play a major role in morphogenesis as well (Marin-Riera et al., 2018;Milinkovitch et al., 2013).

| Waves of patterning in the development of different ectodermal organs
From this general recipe, various modifications and additions have been described which explain the striking differences between ectodermal organs.A first difference between the distinct classes of appendages consists in the specific nested arrangements of patterning processes.While a few types of ectodermal organs exist as unique structures (horns, nails), many types typically occur in geometrically spaced arrays.These specific spatial patterns can be traced back developmentally until their first placodal stages, during which they seem to be initialized through signals emanating from a morphologically and spatially distinct signaling center.In birds, interactions between Fgfs and Bmps along the dorsal midline lead to a one-dimensional line of periodic mesenchymal condensations (Ho et al., 2019;Inaba et al., 2019;Lin et al., 2006).This midline will then serve as a center from whence propagate signaling waves ("traveling" waves) that, in turn, initialize local feather bud patterning in a centrifugal sequence (Bailleul et al., 2019;Inaba et al., 2019).
There has been some debate about the specific mechanisms that create these waves, with some scientists suggesting mechanical stresses and cellular mechanosensing, or differences in cell density as triggering local symmetry breaks in a coordinated fashion (Ho et al., 2019;Shyer et al., 2017).In teeth, interplay between Fgfs and other epidermally expressed factors that is constrained by Bmp expression in the surrounding non-odontogenic mesenchyme (Peters & Balling, 1999), defines the area overlaying the mandibular arch where the dental lamina will give rise to teeth.Molars are initialized sequentially in an anteriorto-posterior gradient, suggesting a slowly traveling inductive signaling wave (Sadier et al., 2019) which starts from a punctuate, rather than longitudinal, center.These dynamics may also be central to establishing the inhibitory cascade leading to mathematically predictable size ratios among molars (Kavanagh et al., 2007).
Mammalian hair buds, on the other hand, do not appear to emerge along any well-synchronized wave front, but by inducing Turing dynamics almost simultaneously (or in a very fast spreading wave) in many places (Glover et al., 2017;Painter et al., 2021).This might also, unlike in flighted birds, but possibly in a manner reminiscent of the feather bud arrays in flightless birds, cause a less ordered pattern of primordia to emerge (Curantz et al., 2021).On the other hand, it has been suggested that this symmetry break tends to occur earlier in areas around the nascent mammary glands than in the remainder of the integument (Painter et al., 2021).
The initiation of integumental denticles within the elasmobranch skin has been found to show similarities to feather buds (Cooper et al., 2018).In sharks and rays, denticle patterning appears to start from two dorsolateral lines that are later characterized by larger, differently spaced and distinctly shaped denticles.It has been shown that, similar to the avian dorsal midline, antagonistic patterning by BMPs and FGFs is involved in forming the initial primordia along antero-posterior lines.Reptile scales, too, emerge from two dorsolateral lines of spots (Di-Poï & Milinkovitch, 2016).
Intriguingly, we see a similar pattern in developing carapacial turtle scutes, despite their overall divergent features: before the primordia of the more central scutes emerge during development and expansion of the carapacial tissue, the future carapacial ridges (CRs) already show two lateral lines of primordial spots (Cherepanov, 2014;Moustakas-Verho et al., 2014).These spots, which mark the positions of the future marginal scutes, represent a sequence of partially overlapping expression domains of Shh, Bmp, and Gremlin.Although it is currently unknown what mechanisms cause the formation of these patterns, it has been shown that experiments interfering with SHH and BMP signaling will prevent them from forming (Moustakas-Verho et al., 2014).Since no punctuate expression of this prepattern is visible in softshell turtles that have lost the ability to develop scutes in evolution, it has been proposed that this lateral prepattern is required for the formation of the correct carapacial scutation.A mathematical model based on the assumption that the CR serves as a prepattern from which centripetal activation waves will induce the formation of the more central scutes has been able to recreate both natural, environmentally and experimentally caused scute anomalies (Moustakas-Verho et al., 2014;Zimm et al., 2017).Thus, this comparison between different types of appendages strongly suggests a potential ancestral developmental module: a one-dimensional prepattern creating the array of secondary primordia by traveling inductive waves.These lead to local symmetry breaks and the coordinated onset of Turing dynamics, ultimately establishing a spatially regular pattern.
Appendages emerging from secondary primordia tend to be different from the ones emerging directly from the inductive spots, both in a morphological and developmental perspective.In catsharks, primary integumental denticles assume stellate shapes and are overall larger than the more arrow-shaped secondary denticles (Cooper et al., 2017;Cooper et al. 2018).Similarly, squamate reptiles are covered by differently sized scales (Di-Poï & Milinkovitch, 2016, Chang et al., 2009), often corresponding to different waves of patterning during development.In mammals, fur consists of larger and smaller hair types that surround the former (Duverger & Morasso, 2009).These smaller hairs emerge from secondary inductive waves that form around the primordia of the larger hairs, although presence of the latter is not strictly required for initialization of the former (Cheng et al., 2014;Sick et al., 2006).In turtles, secondary patterning leads to the emergence of primordia for the larger scutes, which is presumably caused by inductive waves spreading from the CRs (Moustakas-Verho et al., 2014).Marginal scutes, which issue from the signaling centers along the CR, are substantially smaller.
Following the formation of the array of secondary primordia by inductive waves, further tertiary, or even quaternary, patterning occurs.These patterns tend to be more organ-specific, as their deployment typically involves concomitant local tissue changes, usually growth and tissue differentiation, that take place according to the specific rules of morphogenesis.Notoriously, this is the case for feathers, where subsequent waves of inductive signaling and reaction-diffusion systems have been proposed to pattern different feather types, and then their substructures, such as feather branches (barbs) and smaller ramifications (barbules) branching off the rachis (Chen et al., 2015;Harris et al., 2005;Jung et al., 1998;Prum & Williamson, 2002;Wu et al., 2004).
Even in the carapacial dermis in turtles, another morphogenetic process succeeds the establishment of the array of scute placodes.Emanating from those placodes, centripetal waves define the characteristic scutation pattern by morphogenetic signals that stabilize at their mutual collision sites, forming the inter-scute boundaries (Moustakas-Verho et al., 2014).Although overlapping signaling pathways are believed to be involved in these subsequent inductive events, the described anatomical differences between the resulting structures, that is, size and shape, suggest that these patterning modules are not just repeated iterations of the very same mechanisms.
Thus, we find that the development of ectodermal appendages represents a nested sequence of more or less independent, deployable patterning modules, most of which involve reaction-diffusion systems and growth.From an evolutionary perspective, these modules may have been pivotal for the success of ectodermal organs, as they represent a highly flexible, and adaptive, developmental system capable of readily producing a great variation in shape and size.Their hierarchical organization endowed them to form in a stereotypical and relatively robust manner, while adopting specialized tasks by recruiting modified versions of secondary and tertiary patterning modules during developmental stages of overall tissue growth and differentiation (Table 1).

| Conserved patterning roles for molecular signals
From a signaling perspective, most of the discussed mechanisms recruit a conserved set of pathways.Wnts/ β-Catenin are often found to activate different downstream signals in nascent placodal organs (Kawakami et al., 2001;Noramly et al., 1999); Fgf and Shh are typically associated with activatory and promitotic functions, while they are antagonized by pro-differentiation BMPs (Åberg et al., 1997;Dassule et al., 2000;Jernvall et al., 1998;Jung et al., 1998;Rishikaysh et al., 2014).Central for maintaining the balance between these activatory and inhibitory pathways is EDA and its receptor EDAR (Mikkola & Thesleff, 2003).They have been proposed to exert positive and negative feedback on BMPs, Wnts, Fgfs, and Shh (Pummila et al., 2007).Mutants of the Edar pathway have been found to exhibit strikingly consistent phenotypes in different ectodermal organs.Usually, they will have smaller, sparser, and simpler ectodermal appendages, such as those observed in hair (Botchkarev & Fessing, 2005).While a sparser pattern presumably represents a failure to locally achieve symmetry breaks, smaller organs might result from their insufficient stabilization, leading to smaller placodes.Simplification can, for instance, be seen in teeth of Edar mutant mice, where formation of secondary cusps fails (Kangas et al., 2004), despite conspicuous variation between mutant phenotypes (Charles et al., 2009).For certain ectodermal organs, blocked Edar signaling can even lead to their complete absence (Mikkola, 2008;Mikkola, 2009).Overexpression of the Edar pathway, on the other hand, may cause unusually densely spaced and morphologically altered organs.Increasing Edar signaling in rodent molars has led to the formation of supernumerary cusps (Harjunmaa et al., 2012;Tucker et al., 2004).Thus far, Edar effects have been described in almost all types of ectodermal appendages, including the previously mentioned reptile and fish scales, and feathers (Di-Poï & Milinkovitch, 2016;Harris et al., 2008;Houghton et al., 2005).Thus, it is intriguing to study the role of Edar in the development of turtle scutes.
Previous work (Loredo et al., 2001;Moustakas-Verho et al., 2014) has allowed us to establish the spatial relationships between expression domains of central signaling molecules within the CR.BMP, Shh, and the BMP-inhibitor Gremlin have been reported to form a sequence of alternating and only partially overlapping expression domains, suggesting mutually antagonistic interactions.Furthermore, studies in developing feather buds point to a similar sequence of domains in which these morphogens are expressed (Ohyama et al., 2001, Harris et al., 2005, Noramly et al., 1999, Wu et al., 2018, Patel et al., 1999, schematically pictured in Figure 1), possibly suggesting homology.It is noteworthy that this alternating spatial distribution of three factors requires some anterior-posterior cue, because the pattern is directional and asymmetrical (in other words, an A-B-C-B-A pattern could arise autonomously, without any additional information about global direction, but an A-B-C-A-B-C pattern cannot.Specifically, this input has to be consistent and reproducible, ruling out randomly distributed cues).This pattern is absent in softshell turtles, where signaling molecules are either expressed in a faint line along the CR or not at all (Moustakas-Verho et al., 2014).
T A B L E 1 Overview of proposed nested systems of regular spatial patterns across ectodermal appendages.

Placodes along lateral lines
Prepattern not observed 1. Placodes along dental lamina 1. Placodes along dorsal midline Given the ability of EDA to positively and negatively regulate a wide range of signaling factors involved in orchestrating ectodermal appendage development, we sought to investigate its role in CR patterning.To this end, we assessed the phenotypic effects of interfering with Eda signaling in vitro and in silico.

| Experimental modification of Activin and Edar signaling in the developing CR
Trachemys scripta elegans (red-eared slider) eggs were collected from commercial turtle farms in Louisiana, USA, and processed for ex vivo culture to test molecular function and in situ hybridization to examine gene expression as described in Moustakas-Verho et al. (2014, 2019).For the inhibition of ectodysplasin signaling, stage Greenbaum 15 (G15, Greenbaum 2002) and G16 embryo explants were cultured for 5 days with 2 μg/ml EctoD2, a blocking anti-EDA mouse monoclonal antibody previously shown to phenocopy complete Eda deficiency in mice, and to be successful in inhibiting feather bud formation in chicken skin (Ho et al., 2019;Kowalczyk-Quintas et al., 2014).Embryo explants were also cultured under identical conditions with an agonist anti-Edar antibody (mAbEdar1) able to rescue Eda-deficiency in mice (Kowalczyk et al., 2011), or with an isotype-matched control antibody (Aprily2).To examine the effect of Activin A on turtle scute development, stage G15 and G16 embryo explants were cultured for 5 days with 500 ng/ml activin A protein (courtesy of Marko Hyvönen; Harrington et al., 2006) or bovine serum albumin (1 μg/μl).Effects on scute patterning were analyzed by expression of Shh following Moustakas-Verho et al. ( 2014).

| RESULTS I
Following culture, control G15 and G16 explants formed regularly spaced marginal scute primordia along the developing CR (Figure 2a; n = 14).Explants cultured with mAbEdar1 antibody did not show significant differences from controls, suggesting that this antibody does not cross-react with turtles (n = 9; data not shown).Explants cultured with EctoD2 blocking antibody ranged in effects from complete loss of Shh expression and marginal scute primordia (Figure 2b) to primordia and Shh expression domains that are reduced in size and number, but with regular spacing (Figure 2c,d; n = 18).
G16 explants cultured with Activin A did not show significant differences from controls (n = 10; data not shown).G15 explants cultured with Activin A protein showed effects opposite to those of explants treated with EctoD2: Shh expression domains were expanded and irregular (Figure 2e-g; n = 8).The effect of Activin A on Shh expression in the marginal scute primordia resulted in patterns similar to those obtained after treatment with the FGF receptor inhibitor SU5402 (Moustakas-Verho et al., 2014).Furthermore, these effects also coincided with regard to developmental timing as both Activin A and SU5402 were active in G15 but not G16 explants.

| Mathematical modeling of CR pattern formation by a reaction-diffusion mechanism
One of the most widespread mechanisms involved in the formation of regular arrays of ectodermal appendages are reaction-diffusion systems, typically comprising a set of F I G U R E 1 Key morphogens form a conserved periodic pattern in turtle scute primordia, including the carapacial ridge, and in the dorsal row of feather bud primordia.Reviewing data published in Moustakas-Verho et al. ( 2014) and several studies about morphogen expression in feather buds (see Harris et al., 2005;Noramly et al., 1999;Ohyama et al., 2001;Patel et al., 1999;Wu et al., 2018), we reconstruct a conserved pattern of alternating and partially overlapping expression domains of BMPs (Bmp2 in the turtle), Shh and Gremlin.The double line in the simplified turtle scute primordial pattern represents the epidermis.Green rectangles in the schematic drawings of turtle and chicken embryos on the left correspond to area depicted on the right conserved morphogen pathways.Thus, we implemented a network of key morphogens (Wnts, FGFs, BMPs, Edar, Activin) into a simple mathematical toy model, using a system of partial differential equations similar to those used in well-established gene regulatory network (GRN) models (Moustakas-Verho et al., 2014, Meinhardt & Gierer, 1974; for a more detailed and didactic dissemination of this type of models in the specific case and in developmental biology in general, see also the supplementary file).In our model, the CR outline was simplified as a one-dimensional stripe of 300 connected cells, with closed boundary conditions.Within this stripe, all morphogens i were free to diffuse at their specific diffusion rates D(i).
Changes in morphogen concentration were, at each simulation step, calculated as follows: with c(i) referring to the concentration of gene product i at a given time, μ(i) being the specific degradation rate and u p (i,j) and u n (i,j) being two functions describing the specific, positive or negative, interaction strength between two morphogens i and j.The former function assumes a value of 0 whenever it would be negative, the latter assumes 0 whenever its value would be positive otherwise.Interaction means here that morphogens will affect each other's expression rates in a specific way.We initialized the simulations by setting the concentration of Wnt to 0.5 plus a small stochastic noise term between 0 and 0.1, while all other morphogens were set up with concentrations of zero.This noise term was included to account for ubiquitous molecular heterogeneity.We allowed only morphogen concentrations between 0 and 10 throughout the simulations.
Since pattern formation only occurs for specific combinations of parameters, we first performed a broad screen by exploring the pattern formation capacities of 50,000 randomly parameterized networks.We only demanded a negative feedback between BMP and FGF, as such an interaction has often been suggested (Neubüser et al., 1997) and is key to giving rise to Turing dynamics.We observed that the vast majority (about 98%) of networks were not able to generate spatially stable periodic patterns.Among the remainder, we picked networks with a sufficiently large period length whose dynamics were reminiscent of a Turing mechanism (e.g., the perturbed pattern was seen to spontaneously rearrange into "emptied" spaces).Subsequently, we pruned these networks by removing network connections that were functionally neutral to the emerging spatial pattern.We then went on to see whether any of those networks would show the specific properties of the real CR patterning network that we identified experimentally.
To test whether the model could account for the observations of our experiments (and experiments described in Moustakas-Verho et al. ( 2014)), we added and removed gene product concentrations in the model.This was done by either adding a fixed value or removing a percentage of the present concentration, that is, essentially implementing ectopically increased degradation.Since the effect of the chemicals used in the experiment may have been spatially heterogenous (e.g., due to differential tissue permeabilities), we also altered the gene product concentrations in a stochastic manner.Finally, we settled on one example network with the desired properties, although we found several other, mostly topologically similar, networks that exhibited comparable dynamics.

| RESULTS II
We observed that inhibition of BMP would lead to overall blurred expression patterns.In experiments where BMP inhibitors were added to the culture medium (Moustakas-Verho et al., 2014), blurred or homogeneous stripes of Shh expression were observed, which is in line with the analogous experiment in our network model.
Gradually decreasing EDAR expression created random patterns in silico, and stronger inhibitions decreased the number of stripes until the pattern was spatially homogeneous, but would undergo temporal oscillations.A similar effect was observed after removing FGF (not shown) and adding Activin ectopically.Although arguably weaker, corresponding effects were seen in the above-described experiments.Interestingly, we found that a very high rate of EDAR removal would lead to another, albeit substantially sparser, stripe pattern.Although we decided to emulate the experimental settings by using random terms in the differential equations, very similar results could be achieved when we changed morphogen concentrations in a spatially homogeneous way.In all cases, very strong interference with pattern development prevented the formation of stripes altogether.We display our chosen GRN and the patterns resulting from its normal and perturbed dynamics in Figure 3.

| DISCUSSION
This study is the first to report the phenotypic effects of experimentally interfering with EDAR and Activin signaling in the developing turtle carapace.In accordance with observations in other ectodermal appendages, we find that blocking Edar signaling has an inhibitory effect on the formation of scute primordia within the carpacial ridge.Furthermore, we found a partly stochastic disruption of the regular pattern after experimental enhancement of Activin, a signal that has been previously shown to stimulate EDAR expression (Laurikkala et al., 2002) which interferes with scute patterning on the CR.In addition, it has been shown to be critical for the morphological complexity of teeth (Harjunmaa et al., 2012).However, this effect was not binary, with a multitude of different phenotypes (Charles et al., 2009).
Whereas ectopic increase of Activin at an early developmental stage led to more severe phenotype changes, typically a stark reduction of CR primordia numbers in a randomized manner, increased Activin at a later stage did not seem to affect the patterning of scute primordia.This suggests that the function of Activin signaling is restricted to specific developmental stages.In EDAdeficient mice, tail hair never form, but administration of recombinant EDA at birth, followed 4 hr later by an inhibitor of EDA, was sufficient to induce formation of numerous and permanent tail hair, suggesting that placodes can be stabilized and development of ectodermal appendages can be irreversibly induced in a matter of hours (Swee et al., 2009).Thus, if exogenously added Activin acts by upregulating EDAR expression in this system, its action may become indiscernible if applied after the time of natural EDAR activation by EDA.We have previously shown FGF signaling to be important for scute patterning along the CR (Moustakas-Verho et al., 2014); we suggest that EDAR may play a similar, albeit possibly more significant, role during scute development.This would be consistent with suggested roles of EDAR in the development of other ectodermal organs (Botchkarev & Fessing, 2005;Di-Poï & Milinkovitch, 2016;Harris et al., 2008;Houghton et al., 2005;Kangas et al., 2004).Full removal of a periodic pattern can be considered the most severe phenotypic consequence of pattern reduction and simplification.Thus, the role of EDAR in the CR seems consistent with its function in other ectodermal organs, such as teeth.As the punctuate pattern along the CR can be considered homologue to primary patterns (i.e., relying on mechanisms shared with the first wave of patterning in other ectodermal organs), we can compare its removal by EDAR inhibition with the complete lack of the first wave of guard hairs in the Tabby (Eda-deficient) mutant mice (Mou et al., 2006).Unfortunately, the limited range of experimental procedures that can be performed in turtles has not allowed us to observe the effect of blocking EDAR on downstream patterning, such as the formation of the more central scutes.We speculate that, in analogy to the simplification and incomplete differentiation of secondary patterns, such as hairs and tooth cusps (Kangas et al., 2004;Pispa et al., 2004), in the respective edar mutants, one might expect the formation of fewer and incompletely keratinized central scutes.
We found that a simple toy model can reproduce these nontrivial consequences of inhibition of EDAR, alongside other experiments.We suggest that the randomization and simplification of patterns we observe upon inhibiting EDAR both in the model and in experiments can be explained by an insufficient stabilization of the regular pattern created by Turing dynamics.In absence of sufficient edar pathway activity, it appears that molecular noise drives pattern organization, as we see chaotic patterning even when we model the reduction of EDAR concentration without a random term.This suggests that, in this case, the low noise levels added in the initial gene expression of Wnt were sufficient to destabilize pattern formation.Such an irregularity of phenotypes can be compared to the diversity of tooth shapes in mice with different mutations in the Edar pathway (Charles et al., 2009).Ultimately, it may stem from the nonlinear relationship between morphogen concentration and phenotype, making the latter highly sensitive to some ranges of morphogen concentrations, but very robust to others (Green et al., 2017).
In extreme cases, we observed oscillations of spatially homogeneous morphogen concentrations in the model.This finding can be interpreted as EDAR signaling providing the necessary activation needed for the transition In each case, a random value a between 0 and b was multiplied with delta in each cell and time step and either removed from, or added to, the current concentration of the targeted gene product, while ensuring no value was below zero.For BMP: b was À30.0 and À40.0, respectively, from above.For Edar: b was À7.5, À7.77, À8.0, from above.For Act: b was 2.12, 2.2, 2.225, 3.0, from above.We indicated the different degrees of interfering with signaling in the graphic by different numbers of "+" and "À" above the respective graphs.Higher values of b for Act and values of b between À7.8 and À100 for EDA(R) typically led to the loss of patterning and the emergence of mostly oscillatory, spatially homogeneous, behavior between two different regimes of signaling dynamics, one stable in time but unstable in space, the other unstable in both.We do not claim that the GRN we found by parameter exploration represents the actual interactions of morphogens in the CR; yet it may capture the relationship between the edar pathway and the key signaling kernel consisting of BMP and its activators that is pivotal for the emergence of Turing dynamics.
It is intriguing to compare the absence of CR pattern upon blocking the edar pathway to its loss in softshell turtles.Although we do not know yet all developmental differences that were critical at the evolutionary origins of those taxa-those may now be superseded by further modifications in developmental systems by drift-we suggest that preventing the stabilization of activators is a parsimonious way to ablate the complex pattern of scute primordia.However, as there are always multiple ways of simplifying complex patterns, alternative hypotheses cannot be ruled out.
Regarding the developmental mechanisms underlying the formation of the punctuate CR prepattern itself, we suggest three hypotheses, as summarized in Figure 4: 1. Regular alternating patterns are a common feature of classic Turing mechanisms (Meinhardt & Gierer, 1974).Since Turing mechanisms have already been established as central to the subsequent scute formation, it would be relatively parsimonious to suggest that those signaling circuits that give rise to Turing dynamics at later stages may have become involved in patterning at earlier stages.Conversely, the pattern spacing in the CR is very different from the more medial scute primordia pattern, which points to a largely independent mechanism.Interestingly, fossil specimens of prehistoric turtle carapaces display a large amount of diversity in the central scutes, which is contrasted by an overall conserved pattern of marginal scute rows (Ascarrunz & S anchez-Villagra, 2022).Thus, the respective developmental mechanisms underlying central and marginal scutes are likely to involve at least some important differences.2. Alternatively, the CR prepattern might simply arise from tissues derived from somites that underly the forming CR (Cherepanov et al., 2019).As the regular somitic pattern is well established at the onset of scute development, it might be rather simple to transpose the same pattern onto the dermis.This hypothesis is currently supported by a direct local correspondence between somitic segments and each CR scute primordium (Cherepanov et al., 2019;Moustakas-Verho et al., 2014).However, due to the relatively early inception of CR patterning, only early developmental stages of somitic derivates can be considered candidates.3. Third, we may not discard the option that an independent wave may have acted along the nascent CR, presumably exhibiting clock-and-wave front dynamics (Cooke & Zeeman, 1976).Although this hypothesis might appear prima facie more baroque and less parsimonious, clock-and-wave front mechanisms actually emerge rather easily and readily in random network topologies (Cotterell & Sharpe, 2010;Hagolani et al., 2019;Salazar-Ciudad et al., 2000), suggesting that they may be rather common.Alternatively, an anterior-posterior wave could launch a classic Turing mechanism in a temporally ordered manner.In addition, such a hypothesis would provide an elegant explanation of the consistent antero-posterior polarity of signaling molecule domains within the CR.In other, better understood, ectodermal appendage developments, for example, feather bud initiation, a temporal antero-posterior delay has been observed (Jung et al., 1998), which may be in line with this hypothesis, without excluding the others.Such a temporal sequence in which CR spots emerge has not been described yet in turtles.
Another line of theory has attempted to link posterior biases in scute anomalies to local developmental instabilities, which may, ultimately, be linked to developmental sequences ("dovetail syndrome," Ewert, 1979).We have previously shown that a model can explain these and other spatial biases without presupposing a morphogenetic wave of antero-posterior development, but merely by acknowledging that the central expansion of the carapace in lateral direction imposes a higher susceptibility to local developmental noise (Zimm et al., 2017).In other words, concomitant signaling and morphogenesis make central scute positions less stable than lateral ones.While there are arguments supporting each hypothesis, more studies are needed to add evidence.Examining spatial repartition of Edar mRNA or protein expression during the onset of CR patterning might be informative (Mou et al., 2006).
We have discussed that a hierarchical structure of nested morphogenetic modules, each capable of unleashing Turing-style pattern formation dynamics, is conserved across different ectodermal appendages.One can speculate whether this particular developmental structure is in part to be held responsible for the success of these organs.Stacking developmental modules with commonalities in their generative mechanisms and GRNs is an easy way to increase modifiability and adaptability (Dassow & Munro, 1999).On the other hand, subsequent mechanisms tend to be increasingly organ-specific, building upon and increasing the specific deviations to the common theme that were laid down in earlier stages of development.This is chiefly because the patterning coincides and mutually interacts with morphodynamic tissue growth in 3D and tissue differentiation.For instance, some ectodermal organs tend to grow out from the initial epithelial plane, starting from the placode stage, whereas others tend to grow down into the underlying mesenchyme (Pispa & Thesleff, 2003).
Overall, we observe that the relative similarity of deployed developmental mechanisms between different ectodermal appendages decreases with increasing developmental stages (cf.nested morphogenetic processes; see Table 1, also cf.Wu et al., 2004).In this sense, we can discern the upper part of an organ class-specific "developmental hourglass" (von Baer 1828; Cordero et al., 2020;Irie & Kuratani, 2014;Raff, 1996), in which the relative amount of morphological disparity follows a temporal sequence of developmental events, each of which contributes to amplification of morphological differences.The early stages of placode formation would then correspond to the "phylotypic stage" or "bottleneck" of the hourglass, as the involved mechanisms can be considered rather stereotypical across vertebrate clades.Defining a "developmental hourglass" specifically for a well-defined, albeit hyper-diverse, organ system based on a sequence of nested mechanisms has the advantage that it might run less into some of the problems whole-embryo hourglass hypotheses tend to be plagued with, such as establishing stage homologies across clades, heterochrony of key events and heterogeneous bases of comparison.On an evolutionary scale, this might mean that organs whose development passes through very similar initial stages, will have accumulated developmental differences at later stages that explain their phenotypic disparities.Thus, ectodermal organ development may have been critical in vertebrate evolution due to its ready potential of increasing phenotypic variation without relying on developmental complexity, but on recycling and nesting of a successful morphodynamic theme.
In the developmental context, inputs often consist of (1) a number of cells with a specific spatial arrangement, (2) a network of genes that interact, via their products, and (3) some initial gene product concentrations in the cells.More complex models may also consider feedback between gene products and cell or tissue behaviours (such as cell division, apoptosis, cell adhesion, ECM secretion, cell shapes, etc.), and subcellular structures might be included.The complexity of a model should typically be sufficient to address a specific set of hypotheses, but not more complex than necessary.Model outputs should correspond qualitatively to their inputs.
In this study, we built a model to understand the outcome of two experiments in which ectopic changes in two morphogens (Eda and Act) lead to a substantial and unintuitive phenotypic result.Thus, we included a network of morphogens, including Act and Eda, that would be in line with some of our current knowledge, and some informed assumptions, about their interactions.We included Bmp, Fgf/Shh and canonical Wnt pathways, since there is a rich literature about their interactions in developing ectodermal appendages, and because they are used as reporter of patterns within the CR.In actual cells, interactions of these signalling pathways would usually consist in the binding of gene products to the promoter of their target genes, affecting the expression of the latter.Proteins may either stay within a cell, be presented on its surface, or be secreted and diffuse in the extracellular space, mediating interactions between cells.These different localizations are essential for morphogenesis, since juxta-or paracrine cell communication via secreted morphogens and cell surface receptors can give rise to spatial patterning within tissues.Since our study focuses on a spatial morphogen pattern, without considering morphological tissue changes, we used a fixed set of cells in our model.We also simplified the spatial context by not considering different tissue compartments (namely, epithelium and mesenchyme) and implementing a 1-dimensional row of cells.Note that what we call here "cells" in the model does not necessarily correspond to single cells in reality, but may be "translated" as equally sized minute tissue portions.

THE SPECIFIC IMPLEMENTATION
Here, we use a line of 300 cells with defined positions.In this setup, every cell, except the first and last one, has two direct neighbors.Since this asymmetry could give rise to artifacts, we defined the 1 st cell as neighbour of the 300 th cell, essentially converting the row into a ring (closed boundary conditions).Although the actual CR is not a closed ring, we considered this conformation as the mathematically most symmetric and simple solution.Each cell serves as a spatially defined container of gene products whose concentrations can change during simulations.Gene products, here coined i and j, were allowed to diffuse between neighboring cells, at a pace proportional to their intrinsic diffusion rate.Mathematically, we adjust their concentration between adjacent cells over time by moving a portion of the respective gene product concentration in the neighbouring cell with higher concentration to the neighbouring cell with the lower concentration.This portion i' is the product of the concentration difference and the specific diffusion rate D of the gene product under consideration: Here, we describe how the concentration c of gene product i in cell s changes from time time step t0 to time step t1: ( c(is)t0 → c(is)t1 ).Thus, we indicate the temporal and spatial markers as subscripts.
Gene products undergo also degradation, which is dependent on several biochemical factors.In the model, this is implemented by a term i'' in which a certain percentage of a given gene product concentration is removed at a specific rate μ: i'' = -c(is)t0 *μi Finally, gene products depend on the presence of other gene products (here: j) that either promote or inhibit their expression (promoter binding), modification (cleavage, secondary protein modification), activation (receptor binding by ligand), or persistence (protecting against degradation).In order to simplify this diversity of interactions, we only consider that a change in gene product concentration should be proportional to the concentration of the morphogen that it is affected by, and a particular interaction strength, which can be positive or negative or zero, and is defined for each pair of genes.The interaction strengths u(i,j) used in our specific model are given in the caption of Fig. 3.This behaviour can, mathematically, be expressed as a function f. i''' = f ( c(js)t0,u(j,i) ) We use a specific interaction function f which is derived from similar models (see references in the main text).This function which sums up the positive up and negative gene-gene interactions un has been chosen to emulate the usually non-linear dynamics of protein interactions.Note that we do not distinguish between different intermediate forms of gene products, or nested members of a messenger cascade that serve to amplify signals.
The three discussed components i', i'', and i''' will add up to define the gene product concentration at the subsequent time step: c(is)t1 = c(is)t0 + i' + i'' + i''' These three components that define change in gene product concentration (diffusion, degradation, reaction) can also be comprised and described together by a partial differential equation describing change at any given time for infinitesimal time steps dt.Since computation is not continuous, we define the minimal time step as a simulation parameter delta.While very small step sizes increase precision, they will lead to longer computation times.Too large steps sizes, on the other hand, may produce artifacts.
To account for the ubiquitiously present stochasticity of biological processes, we also use terms adding small random values, specifically in defining initial gene product concentrations.This is important for symmetry breaks to take place and may also ensure sufficient robustness of output by collapsing any potential instable pattern in the process.Experimental modifications are implemented by exogeneously adding, or removing, portions of gene product concentrations.
Simulation outputs consist of gene product concentrations per gene and cell.We also looked at their temporal change to filter out temporally instable, chaotic, or oscillatory behaviours.
The simplicity of the general structure of this model allows for a wide range of additions and modifications.In this study, we add or remove concentrations of specific gene products, formalized as an additive term (R) to the gene concentration changes per time step (to be added to the right side of the partial differential equation above): c(is)t1 = c(is)t1 + R Specifically, we consider both constant and stochastic changes.Stochastic changes are implemented by ranging the term R at random between 0 and a value b, 0 < R < b .This allows for b to become a model parameter, emulating larger or smaller, yet stochastic, changes.Note that b can also assume negative values, emulating the rate of removal of gene products.Overall, this way of modelling the respective experiments resembles increased protein degradation or increased production that is not subject to translational control.The provided code written in fortran90 (with .f90extension) allows to reproduce the in silico results used in this study and generate additional data ad libitum.It can be compiled with a suitable compilation program, such as (but not restricted to) gfortran: gfortran PROGRAM.f90-o EXECUTABLE Items written in capital letters are placeholders for the specific names of the program and its executable forms.To run the program, type the execution line described in the code header into the command line (bash):

USAGE OF THE MODEL CODE
./EXECUTABLE INPUTFILE GENE_TO_MODIFY MODIFICATION.
The specific INPUTFILE contained within the supplementary material provides all parameter values to run normal development.If desired, a gene number can be given along the amount of its exogeneous change.
e.g.: ./EXECUTABLECR_GRN0.log 4 -0.5 would run the compiled program with the parameters provided by the input file CR_GRN0.log and exogeneous inhibition of gene 4 (g4=Edar.cf.Fig.3) at a rate of 0.5 (arbitrary units).It will generate an output file containing cell positions and gene product concentrations at several developmental time points.For further, detailed questions, please do not hesitate to contact the authors of this study.

F
I G U R E 2 Experimental manipulation of turtle scute prepatterning.(a) Control Trachemys scripta G16 embryos cultured with Aprily2 show regular spacing in the patterning of scute primordia along the carapacial ridge as seen by the expression of Shh (arrow).(b-d) In G16 cultures with EctoD2 blocking antibody, scute primordia along the carapacial ridge are absent or reduced in size and number (arrows).(e) Control T. scripta G15 embryos cultured with bovine serum albumin (BSA) show regular spacing in the patterning of scute primordia along the carapacial ridge as seen by the expression of Shh (arrow).(f-h) T. scripta G15 embryos cultured with Activin A protein show fusions and absences of domains expressing Shh in scute primordia along the carapacial ridge (arrows).Anterior is toward the top

F
I G U R E 4 We suggest three hypotheses of how the prepattern in the carapacial ridge might emerge.(a) In classic Turing patterns, short-range activators and long-range inhibitors interact to break spatial symmetries and create a stable spot pattern with constant distances.(b) Alternatively, signals emanating from the underlying and already patterned somitic mesoderm might simply induce scute primordia along the carapacial ridge.(c) Frozen wave fronts: This mechanism presupposes the propagation of a symmetry-breaking wave front.A periodic clock (or a Turing mechanism) then creates a regular stripe pattern.