Stripe and spot selection in cusp patterning of mammalian molar formation

Tooth development is governed largely by epithelial–mesenchymal interactions and is mediated by numerous signaling pathways. This type of morphogenetic processes has been explained by reaction–diffusion systems, especially in the framework of a Turing model. Here we focus on morphological and developmental differences between upper and lower molars in mice by modeling 2D pattern formation in a Turing system. Stripe vs. spot patterns are the primary types of variation in a Turing model. We show that the complexity of the cusp cross-sections can distinguish between stripe vs. spot patterns, and mice have stripe-like upper and spot-like lower molar morphologies. Additionally, our computational modeling that incorporates empirical data on tooth germ growth traces the order of cusp formation and relative position of the cusps in upper and lower molars in mice. We further propose a hypothetical framework of developmental mechanism that could help us understand the evolution of the highly variable nature of mammalian molars associated with the acquisition of the hypocone and the increase of lophedness.


Results
Morphometric analysis in fully formed teeth. We reconstructed a 3D enamel-dentine junction (EDJ) model of each UM and LM mouse specimen, which we transformed into a circular 2D image using height from the cervical plane by a morphometric mapping method (Fig. 1). The difference between stripe and spot patterns was quantified by the average ratio of area to the squared perimeter for each object in a binary image. This stripelike pattern should be smaller than that of the spot-like pattern, as the perimeter of each object in a stripe pattern is longer than that of the circular spotted pattern object. To detect a 2D pattern difference, we compared UM to LM while changing the threshold for binarization to within 40-60% of total cusp height. This pattern indicator for UM was significantly smaller than LM, suggesting that UM tends to have a stripe-like pattern, but LM has a spot-like cusp pattern.
Mathematical analysis of mode doubling in a region growing system. To understand the pattern formation mechanism, we utilize the Turing model, which is widely used to model tooth development 12 . Before generating complex 2D models, we utilized a simple 1D model with a growing domain to understand the pattern of mode doubling 29 . The pattern of mode doubling can be classified into three categories: insertion (formation of a de novo peak between two peaks), splitting (formation of new peaks by division of preexisting peak), and mode tripling (splitting and insertion take place simultaneously). At first, we tried uniform domain growth with zero flux boundary conditions. We observed that new peaks always are generated from both edges of the region, which differs from the way pattern formation occurs in vivo (Fig. 2a). This is because the whole system is symmetric, and a zero flux boundary condition can induce deviations from ideal patterns, which makes the region become susceptible to the collapse of the pattern induced by growth 29 . Insertion and splitting both occur at the edges (Fig. 2b).
We then introduced the asymmetry of growth patterns into the model. When only one side of the region grows, new peaks always are generated near the growing boundary (Fig. 2c). We use reaction terms that have reverse symmetry in this case. Consequently, it remains difficult to see whether a new pattern is generated due to the splitting of a preexisting peak or by the insertion of new peaks between two peaks. The 2D simulation of the same reaction term without domain growth results in a stripe pattern (Fig. 2d).
To represent the difference in UM and LM morphology, we introduce reverse asymmetry into the system by adding a quadratic term qu 2 . It has been shown that the pattern formation process can be classified into splitting, insertion, and mode tripling depending on the reverse symmetry of the system 29 . We observed that the insertion pattern becomes dominant when a positive q term is introduced (Fig. 2e,f).
In all cases, the number of peaks is proportional to the domain size at a specific time. Therefore, we expect that the longer axis should have more peaks in a 2D system.

Computational modeling of molar morphogenesis.
To obtain data about how tooth germ size increases, we measured the length of the long and short tooth germ axes. We carefully isolated the UM and LM from the gnathic bones between embryonic day 14.5 (E14.5) and E18.5. During this period, the secondary enamel knots appear sequentially and regulate cusp patterning 30 . The long and short axes of the tooth germ were measured at various time points (number of samples, n = 4). To model the growth rates of tooth germs, growth was fit to a sigmoid curve [Eq. (1); Fig. 3]. Our simulation was constructed under the following framework; • The calculation of the elliptical domain is expanded by the growth function.
• Prepare two parameter settings that realize stripe and spot patterns.
• Morphogenesis of the UM and LM was simulated using parameter sets of stripe and spot patterns, respectively.
As a starting point, the model includes the concentration of an activator that is located in the center of the domain. This initial condition corresponds to the primary enamel knot that appears in the cap stage of development (Fig. 4). The model for UM indicates that the first activator concentration appears on the mesiobuccal side, which then progresses in the mesiolingual and distal direction, where it finally results in three wave-like structures (S1 Video). In the LM simulation, the first activator concentration appeared on the mesial side and then shifted toward the lingual direction (S2 Video). Another concentration then emerged on the distal side and split from the lingual to the buccal side, at which point the final concentration was added distally. These simulations of molar morphogenesis are consistent with earlier reports that described the order of molar emergence in UM 31 and LM 22,32 , respectively. The Turing model involved in a region growing system based on empirical data showed that UM and LM could be described by a parameter set that yields stripe and spot patterns, respectively. Our model successfully reproduces both UM and LM cusp patterns, as well as the actual process of tooth development.

Discussion
In this study, we tested the hypothesis that differences in UM and LM morphologies are explained by stripe and spot patterns using mice as a model. First, we tested the hypothesis with a "static" analysis of tooth shape. Results show that UM exhibits smaller index values (area/perimeter), which indicates a more complicated shape pattern compared to LM. Second, we tested the hypothesis with a "developmentally dynamic" analysis using the Turing model as a framework. The model revealed that the parameter settings for stripe vs. spot patterns matched the UM vs. LM morphologies. Furthermore, our tooth model reproduced the relative location of cusps as well as Turing patterns with domain growth have not been well studied until recently because the growth of material is not a common phenomenon in physical or chemical systems. In 1995, Kondo and Asai 7 described that the stripe pattern in angelfish could change due to the growth of the fish, which invoked interest in pattern formation in growing domains. Pattern formation in a growing domain has been analyzed theoretically in one-dimensional systems 29 . In particular, they begin by numerically implementing the Turing model on a growing domain and observe the increase of activator peaks due to the domain growth, which is called mode doubling. The steadystate solution of the Turing pattern on a growing domain was analyzed using a piecewise linear approximation of the reaction term, and the condition for the distinction of these three dynamics turned out to be correlated to the reverse symmetry of the system 29 . Intriguingly the reverse symmetry of the system also is important for stripe-spot selection (in this study, we use the term 'selection' to refer to the choice between stripe and spot patterns in the context of the Turing system but not to natural selection in evolutionary theory) 4 . We showed the relationship between stripe-spot selection and mode doubling (Fig. 2), which has not been fully elucidated analytically.
Previous studies have simulated tooth development using a morphodynamic mechanism in which inductive and morphogenetic mechanisms interact dynamically with each other 12,33,34 . While tooth development has been considered highly self-regulated, a recent study showed that external factors, such as physical interaction with the jawbone, is relevant for the regulation of tooth morphogenesis 35 . In this study, we simplified the model by not considering this potential external factor. Our model is implemented by two independent mechanisms: (a) global size regulation based on empirical data, and (b) local cusp patterning based on a Turing model. We thus call it a semi-morphodynamic model. While the relative position of the cusps is altered moderately, the sequence of cusp formation remains essentially unaffected with or without the external factor. While our model does not allow We could observe the formation of new peaks restricted to the boundary. (b) The distribution of u (blue) and v (orange) at the timepoint of new peak formation. Splitting (arrowhead) and insertion (arrow) are observed simultaneously. (c) Pattern formation when growth is restricted to the right edge and the reaction term is symmetric ( q = 0 . Detailed explanations of q is in "Materials and methods" section). A new peak is always generated at the site of growth. It is not clear whether splitting or insertion occurs (circle). (d) Two dimensional pattern formation of the system (q = 0). The stripe pattern is formed. (e) Pattern formation with a nonzero q term (q = 0.45) when growth is restricted to the right edge. Insertion becomes predominant (arrows). (f) 2D pattern formation of the system with a nonzero q term (q = 0.45). The spot pattern is formed. We tested how the overall size regulation of the tooth germ affected the morphogenesis of EDJ. The swapping experiments in silico between the growth function in the UM and the reaction-diffusion parameter set for LM (and vice versa) show that the spatiotemporal cusp patterns are not changed drastically (S3 Video and S4 Video). One implication of this swapping simulation is that tooth morphogenesis might be more sensitive to local 2D patterning than to global size regulation. This indicates that the increased body size and associated change of the growth pattern along the course of evolution have relatively minor effects on the pattern of cusp formation. This is consistent with the notion that rats and mice exhibit similar tooth morphologies despite more than two-fold   www.nature.com/scientificreports/ differences in body mass. It should be noted, however, that slight modifications can occur with such swapping. For example, the UM model that grew by LM growth rates yielded additional distal cusps, which might also be consistent with the patterning cascade model, where variation in small-sized cusps should be cumulative, and random variation should appear easily in later-developing cusps 36 . The related question is to identify the morphogens that act in the system described here. Actual candidates for activators and inhibitors of Turing systems have been proposed in various biological systems. In most cases, the molecules that fall into some of the most important categories of the signaling pathway for development, such as fibroblast growth factors (FGF), bone morphogenetic proteins (BMP), sonic hedgehog (SHH), and Wnt 37 , correspond to it. As we mentioned previously, in the case of tooth development, it has been proposed that BMPs act as activators, and FGF and SHH act as inhibitors 12 . These proteins also should be relevant for determining the morphology of EDJ. It remains unknown, however, exactly how each potential morphogen is regulated during the EDJ morphogenesis. Our model consisted of only two morphogens, which is clearly an oversimplification. Thus, it is difficult to apply this hypothetical model directly to certain molecules that function during real odontogenesis. Our simple model only describes the phenomenon of interest, but it may provide additional insight into morphological evolution, and it provides a better understanding of underlying developmental mechanisms 38,39 .
We may extend the framework presented here to human molars. Fig. S1 shows the EDJ of the first UM and LM of humans. The sequence of cusp formation is shared in both UM and LM in humans: mesiobuccal → mesiolingual → distobuccal → distolingual 40 . The UMs of humans exhibit a ridge, called the oblique crest, which connects the mesiolingual and distobuccal cusps. The grooves that separate these cusps are only weakly expressed or not present at all. In contrast, the LMs of humans have cusps that are delimited clearly from each other by deep grooves. Given these morphological differences between UM and LM in humans, we suggest that stripe-spot selection could also apply to human molars. While this hypothesis remains to be tested, our computational model could provide new insights into the morphogenesis of human molars. For example, recent studies on genetic disorders, such as ectodermal dysplasia, have identified genes involved in the regulation of cusp number and shape 41,42 . Currently, it is unknown whether tooth malformation is associated with changes in the sequence of cusp formation and/or it is related to the conversion between stripe-and spot-like patterns of cusp formation. Most genetic studies in mice, however, have shown reductions in the size and shape of teeth, and it is difficult to increase tooth complexity without modifying multiple signaling pathways 14,43 . However, changes in a particular signaling molecule, such as overexpression of Edar, can result in an increased number of spiky cusps 44 , which might be correlated with parameter changes in the Turing model. Thus, simulation-aided approaches have the potential to link experimental studies using model species, such as mice, with clinical research in humans, which might aid in the prevention and treatment of tooth malformation.
Our modeling suggests the difference between stripe-like UM and spot-like LM in cusp patterning. We hypothesize that stripe vs. spot patterning holds as a general rule for the evolution of mammalian teeth. For example, consider the proboscidean molar morphology. Between the Miocene and Holocene periods, the molars of ancestral 'gomphotheres' possessed spot-like conical cusps arranged in transverse rows on the crown, while recent Elephantidae exhibited more stripe-like incised ridges for shearing and cutting 45 . Thus, it is possible to detect evolutionary changes in molar morphology change utilizing 2D patterning in a Turing model, as implemented in this study. This hypothesis, however, should be tested using phylogenetic analyses based on fossils, computational modeling, and experimental reproduction from spot-like to stripe-like teeth.
Distinct pattern formation may be associated with dietary habits. Figure 5 shows that herbivores tend to have stripe-like teeth that are characterized by well-developed intercusp ridges, including the lophodont, loxodont, and selenodont. On the other hand, carnivores have spot-like cuspidate teeth, such as the carnassial (secodont) or denticulate. Omnivores are in between these two extremes. Although some omnivorous species may have different patterns between their UM and LM, such as mice (and perhaps humans), a diverse evolutionary pattern appears across mammals that roughly corresponds to their dietary adaptations.
Several characteristics of dentition are of special relevance for highly diversified mammal dietary habitats, such as relative molar size 18 and crown surface complexity 46,47 . It has been proposed that the taxonomic diversification of mammals was associated with the gain of herbivory, which subsequently led to changes in molar tooth morphology. Specifically, the increase in the number of lophs led to the "complexity" of tooth morphology that resulted in the ability to consume fibrous plant foods 48 . Such "complexity" includes high-crowned teeth with multiple lophs and suggests that this resulted in the loss of intermediate crown types. We propose that an increase in lophedness through evolution could correspond to the switch from spots to stripes in cusp patterning. Along with such a switch of the developmental program, the acquisition of a hypocone could have played a major role. Hypocones are considered to be a key innovation associated with the taxonomic diversification of herbivorous mammals 27 . In the framework of our Turing model, this change is associated with mode doubling.
The acquisition of the hypocone has evolved more than 20 times convergently among various lineages even though there are several options for adding one cusp during odonotogenesis 27 . In terms of mode doubling, the subsequent cusp appears in mode tripling under the stripe pattern cusp formation, which suggests that the hypocone would not necessarily be derived from a certain cusp or structure, but it was acquired independently in the various species in a various approach 15,27,[49][50][51][52] . Our hypothesis of switching from spot to stripe in cusp patterning is consistent with an increase in lophedness and the adaptive radiation of mammals with hypocones during the Cenozoic era.
Although the UM and LM of an individual are under identical genetic control with a common developmental architecture, most mammals have distinct morphological patterns, which suggests that the algorithm applied to UM and LM is different. Our computational modeling implies that final molar morphology could be linked to 2D pattern formation attained only by slight changes in model parameters. This suggests that the evolution of disparate morphologies may not require extensive modification of the developmental process and permits the diversification of molar morphology in mammals.

Materials and methods
Micro-CT data and 3D reconstructions. Unworn molars were obtained from the ICR mice. The first UM and LM were extracted at a postnatal age of two weeks. A total of 20 molars (n = 10 for each of UM and LM) were scanned using a µCT scanner (ELE SCAN, Nittetsu Elex, Japan; housed at Niigata University) with the following parameters: tube voltage: 72 kV and tube current: 11 µA. This resulted in an isotropic voxel resolution of 5 µm. Since the enamel-dentine junction (EDJ) in the fully formed tooth can be used as a proxy for the final configuration of the inner enamel epithelium that resulted from the patterning phase of development 53 , we selectively reconstructed EDJ for the analysis. Tooth segmentations were made between enamel and dentine, and reconstruction of the 3D model was performed with Amira (FEI Visualization Science Group).
2D pattern extraction and analyses. We evaluated the difference in cusp shape between UM and LM by extracting their 2D patterns. To normalize the different shapes of UM and LM, we projected the 3D surface to a normalized circular image. A 3D surface model of EDJ is transformed into a 2D circular image by means of morphometric mapping 54,55 . Before mapping, each tooth was aligned horizontally at its cervical line and centered using the centroid of the cervical line. The height from the cervical plane was used as a morphometric dataset for the analysis, which was sampled over the entire EDJ surface. To extract the 2D patterning, the morphometric map was binarized by using different values of thresholds. We used thresholds ranging from 40 to 60% of total EDJ height from the cervical plane because this threshold should be sensitive to pattern identification. For example, higher thresholds would only represent information on cusp tips, while lower thresholds represent the shape of the basement. Using the binarized image, we then calculated a dimensionless index of the area divided by the squared perimeter for each object, which allowed us to detect the 2D pattern difference between stripes and spots. If the pattern is spot-like, this variable should be close to that of a circle (1/4π). On the other hand, the stripe-like pattern should have smaller values because each object exhibits a complicated shape whose perimeter is relatively larger than its area.

Measurements of tooth germ growth rates.
All animal experiments were conducted in compliance with ARRIVE guidelines. The protocol for experimentation was approved by the Institutional Animal Care and Use Committee (Approval no. 27-044) of Iwate Medical University, and all methods were performed in accordance with relevant guidelines and regulations. The development of UM and LM was observed from embryonic day 14.5 (E14.5), when it starts the invagination of the epithelium and condensing mesenchyme to form its cap, until stage E18.5, when the crown cusp pattern is settled 21 . Although most cusps develop during the early bell stage, tooth mineralization has not begun. At noon of the day, when the vaginal plug was observed was considered E0.5. The UM and LM were excised from ddy mice (Japan SLC, Shizuoka, Japan) and fixed in 4% paraformaldehyde. Assuming the tooth germ is an ellipse, the long and short axis was measured every other day (Fig. 3). The Gompertz double exponential model was fitted to mean values of growth along the long and short axis from E14.5 as follows: Numerical simulation of Turing pattern in a 1D growing region. One-dimensional mockup models were implemented using Mathematica. The reaction-diffusion equations were a modified FitzHugh-Nagumo system 56,57 to allow the mathematically simple switch between stripe and spot patterning with quadratic and cubic terms 28 . The equations are as follows: where u and v are standardized activator and inhibitor concentrations, respectively. The meanings of the linear parameters are described in Fig. 6a. f u represents the positive feedback of the activator. f v represents the inhibition of activator production by the inhibitor. g u represents the promotion of inhibitor production by the activator. g v represents the decay of the inhibitor. d u and d v are the diffusion coefficients of the activator and the inhibitor, respectively. We confirmed that the parameter set gives vault-like dispersion relation (k) (Fig. 6b), indicating that the system has the pattern formation capacity. c represents the saturation of the activator. Since the existence of positive feedback, activator concentration tends to go to infinity ( c = 0 in Fig. 6c). Therefore, we set the upper and lower limit of the activator concentration using this term ( c = 1 in Fig. 6c). q represents the asymmetry of the u reaction term. When q is positive, the equilibrium point is closer to the lower limit than the upper limit. This corresponds to the lower shift of the equilibrium point, which can be implemented by several situations: 1. Constant removal of the activator (for example, by an extracellular molecule that captures the activator). 2. Additional expression of pseudoreceptor of the activator.
Similarly, when q is negative, the equilibrium point is closer to the upper limit than the lower limit. This corresponds to the upper shift of the equilibrium point. There are several biologically plausible situations to implement this effect.