The human EDAR 370V/A polymorphism affects tooth root morphology potentially through the modification of a reaction–diffusion system

Morphological variations in human teeth have long been recognized and, in particular, the spatial and temporal distribution of two patterns of dental features in Asia, i.e., Sinodonty and Sundadonty, have contributed to our understanding of the human migration history. However, the molecular mechanisms underlying such dental variations have not yet been completely elucidated. Recent studies have clarified that a nonsynonymous variant in the ectodysplasin A receptor gene (EDAR 370V/A; rs3827760) contributes to crown traits related to Sinodonty. In this study, we examined the association between the EDAR polymorphism and tooth root traits by using computed tomography images and identified that the effects of the EDAR variant on the number and shape of roots differed depending on the tooth type. In addition, to better understand tooth root morphogenesis, a computational analysis for patterns of tooth roots was performed, assuming a reaction–diffusion system. The computational study suggested that the complicated effects of the EDAR polymorphism could be explained when it is considered that EDAR modifies the syntheses of multiple related molecules working in the reaction–diffusion dynamics. In this study, we shed light on the molecular mechanisms of tooth root morphogenesis, which are less understood in comparison to those of tooth crown morphogenesis.


Results
The association of EDAR V370A with tooth root morphology. In the observation of maxillary tooth roots, we found single and two roots in UP1s and UP2s, and single, two, three, and four roots in UM1s and UM2s ( Fig. 1 and Supplementary Table S2). In some individuals, the UP2s on both sides were congenitally absent (2.6%). The majority of individuals had single-rooted UP2s (91.2%), three-rooted UM1s (97.8%), and three-rooted UM2s (86.4%) on both sides. For UP1s, the major phenotype was single roots on both sides (65.1%), but individuals with two or more roots were also commonly observed. Regarding the mandibular tooth roots, we observed individuals with two-rooted (75.3%) and three-rooted (19.2%) LM1s on both sides, and some had a heterogeneous (two/three) phenotype (5.4%) ( Fig. 1 and Supplementary Table S3). The most frequent phenotype for LM2s was two roots on both sides (56.5%), but individuals with C-shaped roots were also commonly observed.
In the genotyping of 255 individuals for EDAR V370A, we found 25V/V homozygotes (9.0%), 120V/A heterozygotes (47.1%), and 112 A/A homozygotes (43.9%). Using Fisher's exact test, we observed significant associations between the EDAR genotype and the root shape of UP1s and LM2s (Table 1). Logistic regression analysis including age, sex, and region as covariates showed that the EDAR genotype was associated with the root shape of LM1s as well as UP1s and LM2s (Table 2). Interestingly, the 370A allele was associated with a decreased and increased number of UP1 and LM1 roots, respectively, while the allele was associated with the C-shaped root of LM2s. We also confirmed that there was no association between the region and tooth root morphology, but there was a significant association between age and the root number of UP1s.
We also observed the tooth crown traits such as shoveling in UI1s, Carabelli cusp in UM1s, and the cusp numbers of UM2s and LM2s (Supplementary Table S4), and calculated the Spearman's rank correlation coefficients between dental traits ( www.nature.com/scientificreports/   (Table 3). Among the root traits that were significantly associated with the EDAR genotype, the LM2 root trait was significantly correlated with UI1 shoveling, but the UP1 and LM1 root traits were not. It is worth noting that the correlations between the crown and root traits of the same tooth (UM1, UM2, and LM2) were not significant. Additionally, regarding the root numbers of two adjacent teeth, a positive correlation was observed in UP1 vs UP2 and UM1 vs UM2, whereas a negative, but not significant, correlation was observed in UP2 vs UM1 and LM1 vs LM2.
Computer simulation assuming a reaction-diffusion model. The effects of EDAR V370A on tooth root morphology cannot be interpreted in a simple way since it varies depending on tooth type. For further understanding of the effects of EDAR V370A, we performed a computer simulation for two-dimensional cell patterns assuming a reaction-diffusion model. In this simulation, we hypothesized that (1) the EDAR polymor- www.nature.com/scientificreports/ phism positively modifies the induction of both activator and inhibitor in a pure reaction-diffusion system that determines the cell patterns in root development ( Fig. 2A) and that (2) the areas that become future roots or future furcation areas in the apical surface of dental papilla are determined by the cell pattern at the beginning of root formation (Fig. 3). Our simplistic model can produce various spotted or reverse spotted patterns that can be interpreted as various root types, including one root, two roots, three roots, four roots, and C-shaped root ( Fig. 2B-E and Supplementary Figs. S1-S4). The root shape pattern depends on parameters α s and γ, which represent the strengths of induction for activator and inhibitor, respectively ( Fig. 2B-E). Interestingly, α s and γ have opposite effects; the root number increases as α s decreases or γ increases, and C-shaped root patterns appear when α s is high or γ is low. As the size of the area (i.e. total cell number) of the apical surface increases over time, the root number increases. In addition, the cell pattern can be flipped depending on u max , the saturation value of the activator. When the equilibrium value of the activator u 0 is closer to 0 than to u max (i.e. u max = 10u 0 ; Fig. 2B,D), areas where the activator is concentrated form spots. However, the pattern is reversed on the condition that u 0 is closer to u max than to 0 (i.e. u max = 1.1u 0 ; Fig. 2C,E). In case of spotted patterns, it can be interpreted that the activator in the reaction-diffusion system also functions as a signaling molecule that induces mesenchymal cell proliferations to form roots. In contrast, reverse-spotted patterns indicate that the activator functions as a molecule that causes the inhibition of mesenchymal cell proliferation and the promotion of furcation.

Discussion
Previous studies have shown an association between the EDAR 370A allele and crown characteristics, such as incisor shoveling and double shoveling, an increased number of LM2 cusps, and an increased crown size, which are typical phenotypes in Sinodonty 7-9 . However, the association between the EDAR 370A allele and root characteristics has not yet been examined. In this study, we observed the root morphology using CT images and found that the EDAR 370A allele was significantly associated with the root phenotypes related to Sinodonty, i.e., single-rooted UP1s, three-rooted LM1s, and LM2s with a C-shaped single root. Worthy of attention here is that the effects of the EDAR 370A allele differed between the three tooth types: the 370A allele was associated with a decreased and increased number of roots in UP1s and LM1s, respectively, and an unseparated root in LM2s. We confirmed that UI1 shoveling is strongly correlated with the EDAR genotype 7,8 . Although the LM2 cusp number had been reported to be associated with the EDAR genotype 8 , their correlation was not significant in the present study. This study showed that the associations of the EDAR genotype with UM2, LM1, and LM2 root traits are weaker than that with UI1 shoveling. In general, correlations among the EDAR-associated crown and root traits were lower than their correlations with the EDAR genotype (Table 3). In addition, no significant correlation was observed even between the crown and root traits of the same tooth. These results suggest how complicated the dental morphogenesis is.
In root formation, HERS is thought to play important roles [51][52][53] . HERS is derived from the epithelium at the cuff of enamel formation, which is referred to as the cervical loop. The first step of root formation occurs through the elongation of HERS and interactions between HERS and mesenchymal cells. Preodontoblasts are co-localized and make contact with HERS, and then differentiate into odontoblasts that secrete radicular dentin. In the multi-rooted teeth of mice, HERS shows tongue-shaped inward protrusions that join to form furcation [54][55][56][57][58] . However, the mechanisms of the HERS invagination into furcation areas have not fully elucidated. A recent study has suggested that tongue-shaped epithelial protrusions in the cervical loop, termed cervical tongues, appear before HERS is formed, and that the positions of the cervical tongues are partly dependent on the shape of the crown and the positions of cusps 59 . In addition, some have argued that the shape of the cervical loop and the direction of HERS elongation are passively determined by mesenchymal cell proliferation in the dental papilla. Mesenchymal cell proliferation is active in root-forming areas, but not in furcation areas where HERS can be elongated horizontally at the border between the dental papilla and the dental follicle [59][60][61] . In human and rat molars, unlike murine molars, the formation of dentin islands appears at the center of the furcation area in the dental papilla independently of the surrounding dentin [62][63][64] . Then, the furcation area is formed by fusion of dentin islands and the surrounding dentin. www.nature.com/scientificreports/ In this study, we assumed that the root morphology is determined by self-organizing cell distributions in the apical surface of the dental papilla at the beginning of root formation (Fig. 3). Several patterns of root morphology, including the C-shaped root, could be reproduced by a pure reaction-diffusion model, in which we implemented the growth of the area of the apical surface but not the invagination of epithelial cells (Fig. 2). We recognize that our model is simplistic and the real system must be more complicated. To pursue realistic models, we need to consider growths, migrations and interactions of epithelial and mesenchymal cells and various interactions between multiple numbers of molecules. However, even a simple model can be expected to be useful for understanding the essential molecular dynamics in the tooth root morphogenesis. The inferences drawn from the simulation study were as follows: (1) spotted and reverse-spotted patterns can be produced depending on the saturation value of the activator (u max ); (2) an increased size of the apical surface at the beginning of root formation is associated with an increased number of roots; (3) the strength of activator synthesis (α s ) is negatively associated with the root number; and (4) the strength of inhibitor synthesis (γ) is positively associated with the root number.
So far, several signaling pathways, including Fgfs, Bmps, Shh, and Wnts, have been proposed to be involved in the reaction-diffusion system in studies of hair, feathers, and tooth crown 40,44,65 . In particular, the Wnt/βcatenin signaling pathway works in the reaction-diffusion system in hair follicle spacing 41 , and is thus a strong candidate in tooth root formation as well. In the initiation of root formation, Wnt10a and β-catenin are expressed in preodontoblasts, odontoblasts and HERS, which are distributed similar to the reverse-spotted pattern in our simulation [66][67][68][69] . Dkk1, which likely acts as the inhibitor in the reaction-diffusion system, is also expressed in these cells 68,70,71 . Therefore, the localizations of both Wnt10a and Dkk1 are consistent with our pure reaction-diffusion model under the condition of a low saturation value of the activator (u max ) (Fig. 2C,E). Other supportive evidences are that Wnt10a is involved in dentinogenesis and mineralization through dentin sialophosphoprotein 66 , and that disruption of the Wnt/β-catenin pathway or overexpression of Dkk1 arrests root odontogenesis during tooth root development 69,71 . In addition, the Wnt/β-catenin pathway has been shown to inhibit dental pulp stem cell differentiation 72 . Another candidate is the system involving BMPs and FGFs 73 . Bmp4 is expressed in preodontoblasts around HERS 74 and it regulates HERS elongation 75 . The real system is likely to be more complicated and to involve multiple signaling pathways and their dynamic interactions 40,48,58,[76][77][78] . Further experimental evidence is required to support that tooth root morphogenesis is controlled through reaction-diffusion dynamics.
At the root formation stage, Edar expression appears in HERS 50 . Although we did not explicitly implement the HERS elongation in our simulation study, it is hypothesized that HERS invagenates depending on the pattern of signaling molecules and EDAR modifies the cell distribution through direct and indirect effects on these signaling molecules [79][80][81][82][83][84][85][86] . Our simulation results provide a possible explanation for the discrepancy in the effect of the EDAR 370A allele between tooth types. Given that the EDAR 370A allele positively modifies both activator and inhibitor syntheses, the net effect can depend on tooth type and determine whether the number of roots increases or decreases: it is possible that the EDAR 370A allele has a stronger influence on activator synthesis than on inhibitor synthesis in UP1s and LM2s, whereas the opposite happens in LM1s.
Additionally, the opposite EDAR 370A effects between LM1s and LM2s may be attributed to the partitioning of the areas for the two adjacent teeth. The first stage of tooth development is the formation of the dental lamina. Then, the dental region is partitioned into compartments for incisor, canine, premolar, and molar regions. Finally, these compartments are further partitioned to form each tooth 43 . In the development of molars, the dental placodes appear sequentially from anterior to posterior and, therefore, an increased size of LM1 may cause a decreased size of LM2 within a restricted space of the compartment. EDAR also participates in the development of the signaling centers of the placodes in molars 45 . In this study, we observed a negative correlation between LM1 and LM2 root numbers. However, the correlation is not significant and weaker than their correlations with the EDAR genotype. These results suggest that the area partitioning between LM1s and LM2s is not a major factor that affects their root morphology.
In conclusion, this study showed that the EDAR 370V/A polymorphism is associated with the root morphology of UP1s, LM1s, and LM2s, but its effects on the root number are not simple. The EDAR genotype is a strong genetic factor that is a determinant of the tooth morphology dimorphism observed in Asian populations, but there must be other genetic factors involved as well. Previous studies have demonstrated that common genetic variants of WNT10A and PAX9 are associated with tooth crown size 87,88 , and these variants may thus also be associated with root morphology. The present study also demonstrated that the examination of the modulation of morphology by a common variant is useful for understanding the mechanism of morphogenesis. The molecular interactions that occur during tooth morphogenesis are not yet fully understood. However, the accumulation Figure 2. Computational analysis of tooth root morphogenesis assuming a reaction-diffusion system. (A) Schema of the activator-inhibitor system. Interactions between two diffusible substances, namely activator and inhibitor, generate a self-organized spatial pattern. In this simulation, the EDAR variant is assumed to be associated with the self activation of activator and the activation of inhibitor by activator. (B-E) A representative result (left) and the summary (right) of 20 independent simulations. The variation in the cell pattern is caused depending on α s (strength of activator synthesis) (B, C) and γ (strength of inhibitor synthesis) (D, E). In addition, "spotted" and "reverse-spotted" patterns are emerged under the conditions that the maximum concentration of the activator is high (u max = 10 u 0 ; B, D) and low (u max = 1.1 u 0 ; C, E), respectively. The blue color indicates a high concentration of the activator. Detailed conditions for computer simulations are described in "Materials and methods", and parameter values used are as follows: α s = 1.5- 1.8 and γ = 1.0 (B, C); α s = 1.  www.nature.com/scientificreports/ of data from continued efforts will enable the complete picture of the molecular mechanisms of tooth morphogenesis to be elucidated.

Materials and methods
Subjects. The subjects were 255 Japanese patients (98 males and 157 females; 20-69 years of age) who underwent cone-beam CT scanning (3D Accuitomo F17D, J. Morita Mfg. Corp., Kyoto, Japan) for the purpose of dental/oral surgery and treatments at the University of the Ryukyus Hospital. In a questionnaire, we asked for the birthplaces of their grandparents, i.e., either mainland Japan or the Ryukyu Islands (Supplementary Table S1). We excluded individuals who had any grandparent(s) originating from countries other than Japan. Plaster casts of permanent dentition were obtained from each subject. Saliva (2.0 ml) was collected and stored using Oragene DNA (DNA Genotek, Ottawa, Canada). We obtained written informed consent from the patients for their participation in this study. This study was approved by the research ethics committee of the University of the Ryukyus and all research was performed in accordance with relevant guidelines/regulations.

Morphological analysis.
Morphological analysis of the CT images (voxel size: 250 μm) was performed using the three-dimensional imaging software Amira ver 6.0.0 (Thermo Fisher Scientific, Waltham, USA). After volume rendering of the Digital Imaging and Communications in Medicine (DICOM) data, we defined the occlusal plane using the mesial contact between both sides of the maxillary central incisors and the mesial buccal cusps of the left and right first molars (Fig. 1). Since UP1s, UP2s, UM1s, UM2s, LM1s, and LM2s have variations in the root number and shape, we observed these teeth in slice images parallel to the occlusal plane at intervals of 1.0 mm. Based on the images, we counted the number of roots and categorized the tooth root phenotypes as single (1), two (2), three (3), four (4), and C-shaped (C) root(s), or congenital missing (m); the phenotype  Tables S2 and  S3). For each tooth type, when either the left or right tooth was absent due to tooth extraction or treatment, we excluded the case from the analysis of the root morphology.
Using the plaster casts from the subjects. we also observed four tooth crown traits (shoveling of UI1s, Carabelli cusp of UM1s, the cusp number of UM1, and the cusp number of LM2) according to the the Arizona State University dental anthropology system 89 . We observed both left and right teeth. When tooth was absent or unobservable on one side, the phenotype on the other side was used. When the phenotypes are different between the left and right teeth, the less frequent phenotype was adopted. SNP genotyping. DNA was extracted from saliva that had been collected and stored in Oragene DNA collection kits using standard methods. Genotyping for EDAR 370 V/A (T1540C; rs3827760) was performed using a Taqman Genotyping Assay (Thermo Fisher Scientific, Waltham, USA).

Statistical analysis.
To examine the association between the EDAR 370V/A genotype and tooth morphology, Fisher's exact test and logistic regression analysis were performed using JMP (SAS Institute Japan Ltd., Tokyo, Japan). In these analyses, the phenotypes for each tooth were dichotomized into two classes by merging or excluding some phenotypes. For instance, when the tooth root shape differed between the left and right sides, the case was classified into the less common phenotype. As covariates in the logistic regression analysis, we included sex (male = 1; female = 2) and region. Region was defined as the number of grandparents originating from the Ryukyu Islands (0, 1, 2, 3, or 4; Supplementary Table S1). The EDAR 370V/A genotype is denoted by the number of the A allele (V/V = 0; V/A = 1; A/A = 2). Spearman's rank correlation and P values were calculated among the EDAR genotype, tooth root traits, and tooth crown traits.
Computational analysis. The numerical calculations were implemented in C, and the graphics of cell networks were made in Mathematica (Wolfram Research Inc., Champaign, USA). In the numerical calculations, two-dimensional cell patterns were simulated based on repeated cycles of three steps: cell division, cell network dynamics, and reaction-diffusion dynamics 90 .
Throughout the calculation, the number of cells increased from 10 to 1000 cells. The cell network system used in this study was similar to that of Prusinkiewicz and Lindenmayer 91 , and has been described elsewhere 90 . Cells proliferate and move by outward turgor pressures as follows: polygonal cells are tightly arranged in a twodimensional space and are separated from each other by a straight cell wall. At a regular time interval T = 50.0, the cell with the largest area divides by the cell wall that passes through the gravity center of the cell in a random direction. The position of vertex i, u i = (x i , y i ) , moves depending on the force acting on vertex i ( F i ), which consists of elastic forces of cell walls connecting to the vertex ( F Sj ) and turgor pressures of surrounding cells ( F Pj ) ( Supplementary Fig. S5). The vertex position changes according to the following equations: where F i is the total force acting on vertex i; j indicates the cell wall connecting to vertex i; F Sj is the elastic force of cell wall j acting on vertex i; l j is the length of cell wall j; l 0 is the constant that corresponds to the rest length of cell wall; F Pj is the sum of P j1 and P j2 , which are effects of turgor pressures acting on cell wall j from neighboring cell j1 and j2, respectively; the strengths of P j1 and P j2 are inverse proportional to V j1 and V j2 (i.e. areas of cell j1 and j2), respectively; and k T = 0.5, k S = 1.0, l 0 = 0.3, and k P = 1.0 are constants. Vertex position (i.e. cellular position) was calculated to the equilibrium state.
The reaction-diffusion dynamics of the activator (u i ) and inhibitor (v i ) in the ith cell is described by the following form of equations: F Sj = k s l j − l 0 F Pj = P j1 + P j2 P jm = k P /V jm (m = 1, 2) www.nature.com/scientificreports/ where α s , α d , α m, β, γ, δ, ε, D u , D v , and u max are constants. Equations (1a) or (1b) includes terms for the synthesis, degradation, and diffusion of the activator or inhibitor, respectively. The activator is induced by itself in the strength α s and repressed by the inhibitor in the intensity β. In contrast, the inhibitor is induced by the activator in the strength γ. The activator or inhibitor decays at the rate α d or δ, respectively. ψ modifies the rate of degradation of the activator in marginal cells, which prevents spot patterns from appearing in the edge of cell networks. The activator or inhibitor diffuses between adjacent cells with the diffusion coefficient D u or D v , respectively.
The equilibrium values of u i and v i in non-marginal cells are described by u 0 = δε/(βγ − (α s − α d )δ) and v 0 = γ ε/(βγ − (α s − α d )δ) , respectively, and u max (> u 0 ) is the constant that corresponds to the saturation value of u i (i.e. 0 ≤ u i ≤ u max ). In each of the reaction-diffusion dynamics steps, numerical calculations using the Euler method were carried out until reaching an almost steady state when the total time Td = 50.0 and the time step dt = 0.02. The parameter values in the numerical calculations are described in Fig. 2. For each parameter set, we performed 20 independent calculations in which the initial values of variables were given as their equilibrium with a random fluctuation of 1.0%.