 Research Article
 Open Access
 Published:
The evolution of anthropoid molar proportions
BMC Evolutionary Biologyvolume 16, Article number: 110 (2016)
Abstract
Background
Developmental processes that underpin morphological variation have become a focus of interest when attempting to interpret macroevolutionary patterns. Recently, the Dental Inhibitory Cascade (dic) model has been suggested to explain much of the variation in mammalian molar size proportions. We tested the macroevolutionary implications of this model using anthropoid primate species (n=100), focusing on overall morphological patterns, as well as predictions made about molar size variability, direct developmental control, and diet.
Results
Of the species sampled, 56 % had centroids that fell within regions of molar proportion morphospace consistent with the dic model. We also found that the third molar had greater variation in size than either the first or second molars, as expected by the model. Some dic model predictions were not supported, however, such as the expected proportion of M _{2}/M _{1} when the third molar is absent. Furthermore, we found that some variability in third molar size could not be explained by the influence of the inhibitory cascade. Overall, we found considerable cladespecific differences in relative molar sizes among anthropoid primates, with hominoids and cercopithecins strongly divergent from dic model predictions, and platyrrhines, colobines, and papionins more consistent with the inhibitory cascade. Finally, we investigated reasons why some clades deviated from dic model expectations. Adaptations for frugivory (e.g., bunodont cusp relief) appeared to be one driver of relatively larger second molars and have evolved independently in multiple lineages of anthropoids.
Conclusions
The dic model explains some of the variation in anthropoid primate molar proportions. However, there are interesting deviations away from this broad mammalian pattern, particularly in hominoids and cercopithecins, which suggest the model is only one of multiple mechanisms determining morphological variability in mammalian teeth.
Background
Mammalian dentition is complex and variable, adapting to the dietary, environmental, and social demands of a species’ niche. Despite the considerable range of variability in dental forms, certain morphological and evolutionary patterns emerge repeatedly in phylogeneticallydistant taxa. For example, the hypocone (distolingual cusp on upper molars) evolved independently more than 20 times [1], hypsodonty (highcrowned molar teeth) evolved in multiple species of ungulates [2], lagomorphs [3], and South American rodents [4, 5], and tribosphenic molars may have evolved twice in mammals [6, 7] (but see [8]). During evolution, teeth are more frequently lost than gained [9–11] and the last developing teeth within a particular dental field are almost always lost first [12, 13]. While many of these patterns likely evolved in response to similar selective pressures [1, 14], there is evidence that developmental constraints limit molar shape [15]. Thus it is reasonable to assume that some of the dental variability of mammalian teeth is constrained by developmental processes.
In the past decade, several studies have explored developmental mechanisms potentially responsible for dental morphological variation [16–19]. Kavanagh and colleagues [16] developed the dental inhibitory cascade (DIC) model to explain the evolution and development of molar proportions. Through experimental alteration of developing mouse molars, the authors suggested that much of mammalian molar proportion diversity could be attributed to a simple, highly conserved pattern of differences in the timing and concentration of molecules that activate (a) or inhibit (i) molar initiation and proliferation. When the enamel knot of the earliest developing mandibular molar (M _{1}) is initiated, it expresses inhibitory signals (including ectodin, wise, Bmp3, and follistatin) that inhibit, or greatly delay, the subsequent formation of M _{2} and M _{3}. A change in the expression or timing of these inhibitory signals was proposed to determine the eventual size of M _{2} and M _{3}, by altering the relative timing of molar initiation.
According to Kavanagh et al. [16], therefore, the patterning of molar sizes works as an inhibitory cascade, so that the timing of M _{1} development cumulatively affects the timing of M _{2} and M _{3} development, as the DIC mechanism operates along the developing tooth row. Decreased M _{1} inhibition allows M _{2} and M _{3} to form earlier and ultimately grow larger than M _{1}, while increased inhibition restricts the size of subsequently developing molars and may eventually lead to their loss. The DIC model can be characterized by the following equation:
where y is molar area, as determined by position along the tooth row (x), and by the relative strength of activators (a) and inhibitors (i), defined by the a/i ratio. Molar areas are derived from Eq. 1 as M _{1}=1,M _{2}=a/i,M _{3}=2(a/i)−1. Thus, balance between activation and inhibition leads to equal sized molars, while decreases in inhibition result in larger distal molar areas with a cumulative effect from M _{2} to M _{3}. Since inhibition is hypothesized to cascade distally along the molar row, it is possible to predict the relative size of each molar using the following relationship:
As a test of the DIC model’s predictions in a macroevolutionary context, Kavanagh et al. [16] demonstrated that most of the molar patterns seen in murid rodents can be explained by alteration of the ratio of activation and inhibition molecules expressed during early tooth development.
One extension of the model is that in morphospace regions of high inhibition (M _{1}≫M _{2}≫M _{3}), we would predict minimal temporal overlap between M _{1} and M _{2} calcification. As the level of inhibition decreases, however, the amount of overlap should increase, so that in morphospace regions with the lowest levels of inhibition (M _{1}≪M _{2}≪M _{3}) temporal overlap would be substantial.
The DIC model highlights the importance of developmental timing in determining dental metrics, by establishing the pattern of hierarchical dependence within a particular developmental system and by asserting broad applicability across mammals (indicating stability under even the most extreme fluctuations in absolute developmental time).
While the model does not predict variation in cusp height, enamel thickness, or eruption schedule, one implication of the model is that molar size proportions are independent from these and other traits determined later in development. For example, our interpretation is that under this model, fruit bats, golden moles, and tree shrews not only exhibit similar molar size proportions, but proportions that were determined by fundamentally the same developmental interactions. The DIC model is innovative in suggesting that development may also drive molar size variation in addition to other variables, such as dietary adaptation, phylogenetic inertia, and allometry [20–22].
Researchers have begun to test the predictions of the DIC model in taxa, such as rodents and carnivores, which have shown rapid expansion of dental morphospace over their recent evolutionary history. A study of Mesozoic and Cenozoic mammals suggested that the DIC model could be plesiomorphic for this clade [23]. Yet, while some studies support the original findings of the DIC model [22, 24], others have presented evidence that calls into question its broad phylogenetic applicability [25–29]. For example, voles show an expansion of M _{1} size that drives molar proportions into morphospace (M _{1}>M _{2}<M _{3}) previously thought unoccupiable [26], while canids exhibit a reduced major axis regression slope (0.45±0.07) between M _{3}/M _{1} and M _{2}/M _{1} that is much smaller than predicted by the DIC model [28]. Though the authors of both the above studies explain deviations from the DIC model as the result of changes in differential evolvability of the M _{1}, the inability of the DIC model to explain all relationships between molar size proportions across Mammalia is worth further study. In particular, the model may not generalize to species with longer life histories, as the decay rates of activation and inhibition molecules may be sufficient enough to lessen the effect of the model [18, 30, 31].
Bernal et al. [29] assessed molar size proportion variation among platyrrhines, an extant radiation of South American primates, and found limited support for the DIC model. While the results of a phylogenydependent regression analysis corresponded with the predictions of the DIC model, an examination of platyrrhine molar proportions using ordinary least squares and reduced major axis regression, as well as an assessment of intraspecific variability in molar proportions, deviated from the DIC model’s predictions. The authors attribute these findings to the highly derived dentition of some platyrrhine clades, including the loss of the third molar in marmosets and tamarins, and reduction of the third molar in Cebinae [32].
We suggest that a broader consideration of the anthropoid dentition reveals multiple other morphological observations that are incongruent with the predictions of the DIC model. For example, great apes exhibit a M _{1}<M _{2}>M _{3} molar pattern, in contrast to the M _{1}<M _{2}<M _{3} pattern of cercopithecoids [33–35]. The M _{1}<M _{2}>M _{3} pattern is difficult to explain under the DIC model. Moreover, M _{3} agenesis in humans, which unlike in dwarfed marmosets and tamarins is polymorphic, has been reported to occur without large changes in relative M _{2} size [36–39].
Our aim in this paper is to test three predictions, and two extensions, of the DIC model using a broad sample of extant anthropoid primates. If anthropoid molar proportions follow the DIC model, this represents further evidence of the model’s robustness across mammalian taxa. If, however, anthropoids deviate from the DIC model’s predictions, it is possible that other factors play an equally important role in determining their relative molar sizes. Previous studies that have sought to test the DIC model’s developmental predictions have relied solely on measurements of molar crown size (e.g., [22, 29]). Anthropoid primates offer a unique opportunity to test the developmental predictions of the DIC model, as data are available on the timing of calcification in successive molars, allowing ontogenetic patterns to be evaluated. To test the model’s predictions we employ an integrative approach that considers dental metrics, development, diet, and function.
In the following section we report the results of testing three macroevolutionary predictions (1–3), and two extensions (4–5), of the DIC model:

1.
Molar area relationships: The relationship among molar areas is: M _{1}=1,M _{2}=a/i,M _{3}=2(a/i)−1. Balance between a/i leads to equal area molars, while increasing inhibition enlarges distal molar areas with a cumulative effect from M _{2} to M _{3}. The relationship between M _{3}/M _{1} and M _{2}/M _{1} molar proportions can best be described by the linear equation: M _{3}/M _{1}=2(M _{2}/M _{1})−1. Three regions of molar proportion morphospace are consistent with the dic model: a) a high a/i region where M _{1}<M _{2}<M _{3}, b) a low a/i region where M _{1}>M _{2}>M _{3}, and c) a region where M _{3}/M _{1}=2(M _{2}/M _{1})−1. For mammals with three molars, M _{2} accounts for onethird of total lower molar area.

2.
Molar area variability: M _{3} is predicted to have the greatest size variability of any lower molar.

3.
M _{3} agenesis: M _{3} agenesis occurs when M _{2} is less than half the size of M _{1}.

4.
Developmental timing: In regions of molar proportion morphospace subject to high inhibition (M _{1}≫M _{2}≫M _{3}) there is minimal temporal overlap between M _{1} and M _{2} calcification, while in morphospace regions of low inhibition (M _{1}≪M _{2}≪M _{3}) temporal overlap is substantial.

5.
Primary dietary category and molar area proportions: Molar proportions can be a measure of diet [16].
Results and discussion
To test macroevolutionary predictions of the DIC model we used previously published lower molar occlusal area data from 100 extant anthropoid primates (Table 1; Additional file 1: Table S1; [40]). We analyzed these data using Bayesian phylogenetic generalized linear mixed models (PGLMM), and employed a Markov chain Monte Carlo (MCMC) approach to sample the posterior distribution of parameter space. Phylogenetic signal for each model was quantified using Pagel’s λ parameter [41]. All models achieved convergence to a stationary posterior distribution, as determined by Hiedelberger and Welch’s convergence diagnostic [42] and low levels (<0.1) of autocorrelation. Our MCMC sampling strategies resulted in effective sample sizes of 10,000 for each parameter. Detailed descriptions of the anthropoid sample, phylogenetic tree, morphometric variables, and regression models used are provided in the Methods section.
To evaluate if DIC model expectations were credible, given our anthropoid sample, we calculated 95 % highest density intervals (HDI) of posterior distribution parameter values. The DIC model specifies point predictions, rather than interval predictions, of parameter values. For any empirical dataset, the probability is therefore extremely small that molar proportions would exactly match DIC model expectations. To provide a more useful way of calculating the probability that DIC model predictions were correct, given our data, we created a region of practical equivalence (ROPE) around each DIC model parameter expectation. A ROPE is an interval that encloses values of the parameter that are, for practical purposes, negligibly different from the point value [43, 44]. We used the ROPE as a decision tool for determining whether DIC model predicted parameter values were credible and/or probable for the sampled anthropoid taxa.
To determine a range of slope/intercept values that might be deemed practically equivalent to the DIC mathematical model’s predictions, we used, as a starting point, experimental evidence reported in Kavanagh et al. [16]. This experimental evidence yielded a slope of 2.024 and intercept of 0.997, thus deviating slightly from the DIC mathematical model’s predictions of 2 and −1. For each PGLMM, we calculated posterior probabilities for ROPEs of several sizes. The inclusion of ROPEs of any size in our analyses was a conservative measure, increasing the chance of DIC model predictions being corroborated, compared with using only the point predictions of the strict mathematical model.
Molar area corrected estimates
The product of maximum linear breadth and length measurements often overestimates molar occlusal area (Additional file 1: Figure S1). We employed a correction (see Methods for details) to make these rectangular areas more comparable to areas calculated by tracing outlines around molar occlusal perimeters. We found that the average difference between outline areas and corrected areas was much smaller (o a−c a: \(\hat {\beta } = 0.079\), 95 % HDI −0.73, 0.99; in mm^{2}) than that between outline areas and rectangular areas (o a−r a: \(\hat {\beta } = 5.73\), 95 % HDI 4.87, 6.61). This indicated that corrected areas provided a better approximation to molar outline occlusal area (Additional file 1: Figure S2).
Molar area relationships
Molar area proportions
The DIC model predicts that the relationship between M _{3}/M _{1} and M _{2}/M _{1} is best characterized by a line with a slope of 2 and an intercept of −1. In addition, two other regions of molar proportion morphospace are also consistent with the DIC model: a) a high a/i region where M _{1}<M _{2}<M _{3}, and b) a low a/i region where M _{1}>M _{2}>M _{3}. We calculated the posterior probability that anthropoid lower molar proportions can best be described by a line with parameters contained within ROPE intervals for the interspecific slope [1.90,2.10] and intercept [−1.10,−0.90]. We also determined whether these ROPEs contained credible parameter values for our anthropoid sample by estimating 95 % HDIs for the interspecific slope and intercept.
We found that 56 % of the sampled anthropoid species’ centroids fell inside regions of molar proportion morphospace consistent with the DIC model (Fig. 1, Additional file 1: Figure S3). However, cercopithecins, hominids, and hylobatids occupied the M _{1}<M _{2}>M _{3} region of morphospace, indicating that their dental phenotypes cannot have originated from alterations in the a/i ratio. As in other mammalian clades [24, 26], we found that the M _{1}>M _{2}<M _{3} region of morphospace was unoccupied, suggesting that some dental morphologies result from developmental mechanisms that rarely occur in mammalian evolution.
For the entire anthropoid sample, we found that the posterior probability (\(\mathbb {P}\)) of the interspecific slope, \(\mathbb {P}(\beta _{1\mathrm {B}} \in [1.90, 2.10]) = 0.10\), or intercept, \(\mathbb {P}(\beta _{0} \in [1.10, 0.90]) = 0.37\), being within the ROPE was low. The ROPE was partly encompassed by the 95 % HDI (1.42, 1.99), while the posterior mean was 1.72 (Table 2). For clades within Anthropoidea there was marked heterogeneity in slope and intercept estimates. Some taxa (e.g., platyrrhines and papionins) had posterior mean slope values closely mirroring DIC model predictions, while others (e.g., hominoids and cercopithecins) had much shallower slopes (Table 2, Fig. 2). Hominoids and cercopithecins both had 95 % HDIs that excluded the ROPE, indicating that DIC model predicted molar proportions are not credible for these taxa. All of the clades analyzed, however, exhibited low probabilities of slope or intercept values being within the DIC molar proportion ROPE (Additional file 1: Table S2, Figures S4–5).
These results indicate that anthropoids likely exhibit a shallower slope for molar proportion relationships than the DIC model predicts, suggesting that on average M _{3} increases in size less than twice as rapidly as M _{2}. Just over onehalf of sampled anthropoid species’ centroids fall within regions of molar proportion morphospace consistent with the DIC model, but several lower ranked clades (e.g., Hominidae, Hylobatidae, Cercopithecini) fall into the M _{1}<M _{2}>M _{3} region of morphospace and display molar proportion slopes that are substantially shallower than the DIC predicted slope. Catarrhines as a whole deviate more from the DIC model, and occupy regions of molar proportion morphospace with lower inhibition, than platyrrhines. Within the New World monkeys, only alouattines exhibited larger distal molars as a result of low inhibition. This variability highlights considerable cladespecific differences in molar proportion relationships among anthropoid primates. While the DIC model may provide a plausible explanation for molar proportions in platyrrhines, colobines, and papionins, the probability is low that the inhibitory cascade is a substantial factor driving relative molar size in hominoids and cercopithecins.
In the above analysis, we tested macroevolutionary predictions of the DIC model related to molar proportion morphospace. However, this analysis is potentially sensitive to the inclusion of ratios as variables in the regression model [45]. Therefore, as a robustness check, we tested two mathematically equivalent predictions of the DIC model that relate to the morphospace of raw molar areas, using regression models that included only raw molar area variables.
M _{2} relative area
For mammals with three molars, the DIC model predicts that M _{2} accounts for onethird of total lower molar area (M _{T}). We defined a ROPE of [0.323, 0.343] around this onethird point value. For anthropoids, we found that the probability of relative M _{2} area being onethird of M _{T} was extremely low: \(\mathbb {P}(M_{2}/M_{\mathrm {T}} \in [0.323, 0.343]) < 0.001\). The 95 % HDI (0.349,0.360) did not overlap the ROPE. This indicates that, for anthropoids, the proportion of lower molar total area attributable to M _{2} is likely larger than 0.333 (Table 3). Several lower ranked clades also followed this pattern of having relatively large M _{2}s. Both hominoids (95 % HDI 0.347,0.384) and cercopithecins (95 % HDI 0.374,0.403) exhibited relatively larger M _{2}s than anthropoids as a whole. In contrast, colobines (95 % HDI 0.327,0.362) and papionins (95 % HDI 0.328,0.356) deviated from this pattern, with 95 % HDIs that encompassed most of the ROPE (Additional file 1: Table S3, Figures S6–7). Overall, these results echo those of the molar proportions analysis. Hominoids and cercopithecins deviate most from DIC model predictions, while papionins, colobines, and platyrrhines are most consistent with the inhibitory cascade. Lucas et al. [46] have previously reported that M _{2} accounts for onethird of M _{T} area in anthropoids. In contrast, our results indicate that M _{2} is relatively larger in Anthropoidea, accounting for slightly more than onethird of M _{T} area, though we note the strong phylogenetic component driving variability in relative M _{2} size.
Molar areas
The DIC model predicts that the relationship among molar areas is: M _{1}=1,M _{2}=a/i,M _{3}=2(a/i)−1. Equal area molars result when the a/i ratio is balanced, while decreasing inhibition results in larger distal molars with a cumulative effect from M _{2} to M _{3}. Since the effect of relative a/i levels cascades distally along the molar field, one implication of this model is that the effect of M _{1} to M _{2} is repeated from M _{2} to M _{3}.
We modeled the tripartite relationship of lower molar areas using path analysis to estimate how much the relationship between M _{1} and M _{3} area was mediated by M _{2} area. The ‘direct’ path (c ^{′}; Fig. 3) was the simple relationship between M _{1} and M _{3}, while the ‘indirect’ path (ab) was the product of a proximal path (a) between M _{1} and M _{2} and a distal path (b) between M _{2} and M _{3} [47]. Our effect size measure was the proportion of total lower molar size variability accounted for by M _{2} acting as mediator: p r _{ m }=a b/a b+c ^{′}. Under a strict interpretation of the DIC model, M _{1} size influences M _{3} size only through the size of M _{2} (i.e., p r _{ m }=1). We defined a ROPE of [0.9, 1] for the DIC p r _{ m } expectation.
For anthropoids, the probability was extremely low that p r _{ m } was within the ROPE: \(\mathbb {P}(\mathit {pr_{m}} \in [0.9, 1]) < 0.001\). The p r _{ m } posterior mean (0.73) was just over twothirds, with a 95 % HDI (0.68, 0.78) that did not encompass the ROPE (Table 4). There was some heterogeneity in p r _{ m } among clades within Anthropoidea. Posterior means for all taxa fell well below DIC model predictions. The 95 % HDIs for platyrrhines and catarrhines excluded the ROPE, though some upper bound highest density limits approached 1 for catarrhine subclades. For cercopithecoids the p r _{ m } was much higher than in other anthropoids (\(\widehat {\mathit {pr_{m}}} = 0.86\), 95 % HDI: 0.75, 0.98) and partially overlapped the ROPE. Overall, the path analysis suggests, therefore, that Old World monkeys are more likely than other anthropoids to have molar sizes determined by an inhibitory cascade mechanism. However, even for cercopithecoids, the posterior probabilities of parameter estimates being inside the ROPE were moderate: \(\mathbb {P}(\widehat {\mathit {pr_{m}}} \in [0.9, 1]) <= 0.42\) (Additional file 1: Table S4, Figure S8).
Only about twothirds of the total effect of M _{1} size on M _{3} size was mediated through M _{2} size in anthropoids; much less than expected under the DIC model. While some of the intertaxon heterogeneity in p r _{ m } may be driven by body size (e.g., it is possible that constraints of the a/i model may be relaxed in largerbodied taxa), both the largest and smallest taxa in our study seem to deviate most from the DIC model expectation of p r _{ m }. Surprisingly, those taxa for whom M _{3} is the largest lower molar (i.e., taxa hypothesized to have the least amount of inhibition) display the greatest mediation through M _{2} (e.g., compare Papionini to Platyrrhini). Overall, the path model either implies that there is a mechanism by which M _{1} size directly affects M _{3} size, or in some taxa a secondary developmental mechanism exists that solely affects M _{2} size.
Renvoise et al. [26] conducted mediation analysis on arvicoline rodents, finding that M _{2} size predicted M _{3} size (b path; Fig. 3) much more reliably than M _{1} size predicted M _{3} size (c ^{′} path; Fig. 3). This is consistent with the DIC model, but Renvoise et al. did not calculate the overall effect of the indirect path through M _{2} (the product of a and b paths), which is a more appropriate summary of the overall mediation effect and a better comparator to the direct path (c ^{′}) [48]. It is therefore difficult to compare Renvoise et al.’s results with our own. However, both analyses suggest that while M _{2} size plays a large role in predicting M _{3} size, factors external to the DIC model are also involved in determining distal molar size.
Body mass and fit/deviation from the DIC model
While not a prediction of the DIC model, the relationship between body size (as it influences absolute developmental timing) and fit/deviation from the DIC molar proportion line is of interest because long life histories (associated with larger body sizes) may provide a mechanism through which developmental constraints of the a/i pathway can be relaxed in a broad mammalian context.
We tested the association between species’ absolute perpendicular distance from the DIC molar proportion line and log body mass. For our entire sample of anthropoids, the 95 % HDI (−4.35,−0.39) for the specificlevel slope was negative, indicating that increased deviation from the DIC model line is associated with a reduction in body size (Additional file 1: Figure S9a). Since callitrichins have extremely small body mass and deviate markedly from the DIC model line, we also fitted this model including only noncallitrichin anthropoids, to see if marmosets and tamarins were driving this negative relationship. We found that removing callitrichins reduced the likelihood of a negative association between body size and DIC model deviation, as the resulting 95 % HDI (−3.50,0.47) overlapped zero (Additional file 1: Figure S9b). In addition, among anthropoids, cercopithecins exhibit the greatest deviation from DIC model expectations and yet are relatively smallbodied. We therefore fitted a third model, including only noncercopithecin and noncallitrichin anthropoids, to see how the smallbodied and extremely deviant cercopithecins were influencing the relationship between body mass and DIC deviation. We found that removing both of these clades resulted in a positive association between body size and DIC model deviation, as the posterior mean of the interspecific slope was positive (4.2), though the 95 % HDI (−0.70,9.04) overlapped zero (Additional file 1: Figure S9c).
Molar area variability
The DIC model predicts that M _{3} will have the greatest size variability of any lower molar. We modeled the difference in average molar area variance (using the coefficient of variation for small sample sizes) between each lower molar type across the anthropoid sample. We found that M _{3} size was indeed more variable than either M _{1} size (posterior mean difference 0.015, 95 % HDI 0.009, 0.022) or M _{2} size (posterior mean difference 0.010, 95 % HDI 0.004, 0.017; λ=0.35, 95 % HDI 0.08, 0.63; Additional file 1: Figure S10). This is consistent with both the DIC model and previously theorized patterns of dental development in mammals (e.g., Butler’s [49] morphogenetic field theory).
M _{3} agenesis
The DIC model makes predictions about expected M _{2}/M _{1} proportions when the third molar is absent. Specifically, M _{3} agenesis is predicted to occur when M _{2} is less than half the size of M _{1}. In anthropoids, we instead found that M _{3} agenesis occurred when the M _{2}/M _{1} ratio was much larger (posterior mean: 0.87, 95 % HDI: 0.64, 1.09; λ=0.82, 95 % HDI: 0.75, 0.89) and that the expected M _{2}/M _{1} ratio under the DIC model had very low probability: \(\mathbb {P}(M_{2}/M_{1} < 0.5 \mid M_{3}^{\mathit {ag}}) = 0.0003\). Within a sample of modern humans (n=66) polymorphic for M _{3} agenesis (Fig. 4, Additional file 1: Table S5), we also found that relative M _{2} size was much larger than predicted under the DIC model in individuals with congenitally missing M _{3}s (posterior mean: 0.92, 95 % HDI: 0.83, 1.01) and that the DIC expected M _{2}/M _{1} ratio was improbable: \(\mathbb {P}(M_{2}/M_{1} < 0.5 \mid M_{3}^{\mathit {ag}}) = 0.0003\). In addition, the mean M _{2}/M _{1} proportion for individuals not exhibiting M _{3} agenesis (posterior mean 0.95, 95 % HDI: 0.88, 1.03) was similar to those with congenitally missing M _{3}s (posterior mean difference 0.03, 95 % HDI: −0.03, 0.09).
Third molar variations are seen with some frequency among primates. Modern humans, for example, are polymorphic in third molar number [50]. Evans et al. [51] have recently proposed that the DIC is a mechanism that explains both the relative and absolute sizes of permanent molars in modern humans, including the relative size of M _{2} when the third molar is absent. Our results demonstrate, in contrast, that modern humans have second molars subequal in size to first molars, whether or not the M _{3} is congenitally absent [52].
In marmosets and tamarins, M _{3} agenesis is likely to have occurred independently three times [32] and there is some reduction of M _{2} relative to other platyrrhines. However, our results agree with those of Bernal et al. [29] that M _{2} is still larger than predicted by the DIC model. One other platyrrhine, the extinct genus Xenothrix, has third molar agenesis, though it too has an M _{2}/M _{1} proportion (0.74) much larger than the DIC model expectation [53].
The departure of callitrichins and modern humans from the expectations of the DIC developmental model suggests that mechanisms other than alterations in the a/i ratio may regulate variation in molar number and relative size. With the exception of certain syndromes, the developmental mechanisms leading to third molar agenesis are unknown, even in humans [18, 50]. There is some developmental evidence that supports a lack of association between relative M _{2} size and M _{3} agenesis, as Cai et al. [54] found that tooth size and tooth number were largely independent in rodents. Overall, anthropoids with third molar agenesis have M _{2}/M _{1} proportions that are much larger than predicted by the DIC model, which contrasts with Carnivora, whose M _{2}/M _{1} proportions range from larger (raccoons) to much smaller (Canidae) than expected under the model, with many forms being congruent with the model [24, 28].
Developmental timing
One extension of the DIC model concerns the relative timing of molar development, specifically that temporal overlap in molar calcification initiation will increase with higher a/i. An assessment of nine primate species for which molar calcification data were available suggests the opposite trend. We predicted that temporal overlap would be highest in species with small M _{1}s. Instead we found that genera with relatively large M _{1}s, such as Varecia, had the greatest amount of temporal overlap between the development of M _{1} and M _{2}, while genera like Macaca, with the smallest M _{1}s, had the least amount of temporal overlap (Fig. 5). Though other factors, such as growth rates of subsequent molars, may also contribute to relative molar proportions, the evidence suggests that species with large posterior molars initiate them long after their smaller anterior molar crowns are complete.
A large M _{3} was associated with greater temporal separation between M _{1} and M _{2} initiation, even when controlling for M _{1} size. Though this pattern could also be consistent with a weak correlation between early stage (e.g., enamel knot formation and cusp differentiation) and late stage (e.g., calcification) developmental events, there is no histological evidence to support marked deviation in the relative schedule of these events. This suggests that, for anthropoids, the DIC model cannot fully explain variation in molar proportions.
Primary dietary category and molar area proportions
Kavanagh et al. [16] suggested that molar proportions could be associated with dietary preference. We tested this hypothesis by estimating whether anthropoid species that primarily eat fruit are more likely to occupy the M _{1}<M _{2}>M _{3} region of molar proportion morphospace than anthropoids that primarily eat other food items. We found that frugivorous anthropoids had much higher odds of occupying the M _{1}<M _{2}>M _{3} region of molar proportion morphospace (Odds Ratio: 9.51, 95 % HDI: 2.73, 32.3; λ=0.78, 95 % HDI: 0.68, 0.87; Fig. 6, Additional file 1: Figure S11). This equates to a 851 % (95 % HDI: 173 %, 3131 %) increase in the odds of an anthropoid species being frugivorous if M _{2} is its largest lower molar.
Four taxonomic groups (Hominidae, Hylobatidae, Cercopithecini, Atelinae) that predominantly occupy the M _{1}<M _{2}>M _{3} molar proportion morphospace have diets that are primarily frugivorous, with lowcrowned bunodont molar cusps. Though some researchers have argued that molar proportions represent adaptations to processing food with different material properties [20, 55, 56], there is minimal evidence to support the M _{1}<M _{2}>M _{3} pattern as an adaptation for frugivory. There is, however, evidence to suggest that different cusp shapes and relief profiles have large effects on processing food with different mechanical properties [57, 58]. Broadly, folivores tend to have highcrowned or lophed teeth that aid the shearing of leafy plant material, while frugivores have lower, more bulbous, cusps that aid in mastication of soft pulpy fruit and the emission of juices [59]. If the strong selective pressure on reducing cusp height acts on genes that are involved in both the establishment of cusp relief and the regulation of molar sizes, then selection for lower cusp height would constrain the developmental pathways through which molars could evolve.
Conclusions
Our results provide only limited support for the hypothesis that anthropoid molar proportions are partially governed by the dental inhibitory cascade model. We found a large amount of variability in relative and absolute molar sizes across anthropoid groups. Some clades (platyrrhines, colobines, and papionins) showed patterns of relative molar size consistent with changes in a/i concentrations, while others (hominoids and cercopithecins) diverged markedly from the expectations of the inhibitory cascade.
While the DIC model may have been the plesiomorphic model for mammalian molar proportions, in anthropoid primates and other taxa it is likely that secondary developmental pathways have influenced relative molar sizes [60]. We argue that the molar proportions of anthropoids result from the influence of two constraints in different directions: 1) cooption of the same signaling molecules (e.g., Shh, PDGF) for different developmental processes throughout dental evolution constrains the relationship between cusp height and relative molar proportion; and, 2) relaxation of constraints between activation and inhibition molecules caused by lengthening of absolute developmental time.
Cooption of genes occurs throughout dental evolution, including the development of toothed jaws [61, 62], regionalization in the dentition [63] and dental complexity [64, 65]. Dental features are interdependent, both within a tooth [66] and along a tooth row [67, 68], and it is likely that a similar relationship exists throughout developmental time between molar proportions and cusp relief. For example, Shh is known to regulate Sostdc1 (a probable inhibitor of posterior molars) and alter the size of the enamel knot, influencing the ultimate size of cusps [17, 69, 70]. Similarly, altering concentrations of PDGF can change both molar size by 17 % and cusp height by 40 % [71]. Thus it is probable that selection for change in cusp relief has a related effect on molar proportions that drives species into regions of morphospace that are unexpected under the DIC model.
We propose that a relaxation of the constraints of the a/i pathway brought on by long life histories may allow the evolution of molar proportions that are inconsistent with the DIC model. Activation molecules are expressed only for a short period of time and these upregulate inhibitory molecules. While largerbodied individuals will typically have longer and more intense expression of both types of genes, there is an absolute limit on how long they can be expressed because of their molecular decay rate [72]. Smallerbodied mammals that have slowed rates of embryonic development, such as callitrichins [73], may similarly be more prone to molar loss.
Thus, while smallbodied and quick developing species such as mice may be greatly constrained by a/i genes, largerbodied taxa may evolve in ways that are not consistent with the DIC model. This possibility reiterates longknown concerns about using a smallbodied species with fast life history to generalize aspects of evodevo across Mammalia [74]. Heterochronic change in the eruption timing of mammalian molars may also relax these developmental constraints, though this has not been tested in the developmental literature. This is potentially the case for many of the subfossil and living lemurids, who have a fast dental development schedule with early eruption of the molars [75], in addition to slowdeveloping callitrichins [29].
In this study, our aim was to investigate macroevolutionary predictions of the DIC model in anthropoid primates, specifically testing whether the developmental mechanisms proposed by the model matched observed morphological variation in relative molar areas. We found some congruence between the results of our analysis and the DIC model. Over onehalf of the sampled species’ centroids fell within DIC model expected regions of molar proportion morphospace. In particular, platyrrhines, colobines, and papionins showed patterns of relative molar size consistent with changes in a/i concentrations. In anthropoids, however, it seems relatively easy to diverge from the strict predictions of the inhibitory cascade, which was particularly the case for hominoids and cercopithecins who occupied the M _{1}<M _{2}>M _{1} region of molar proportion morphospace. Diverse groups of nonprimate taxa, for example canids [28], Notoungulata and Astrapotheria [22], and arvicoline rodents [26], also show divergence from DIC model expectations. Though we found only weak evidence of a positive association between deviation from the DIC model and body size, this could potentially be a result of the majority of the sampled anthropoid taxa having medium body size. A more robust test of this hypothesis will require a broader mammalian sample. While the DIC model explains some of the variation in mammalian molar proportions, there are likely other important developmental pathways being coopted to produce molar rows with larger M _{2}s than expected or those missing M _{3}s. Further studies investigating the ontogenetic and molecular processes underlying different species will shed more light on the specific developmental pathways that determine molar proportions across mammals.
Methods
Materials
We obtained previously published data on lower molar areas of anthropoid primates from Plavcan [76]. This source reports maximum mesiodistal length and buccolingual breadth (separately for the trigonid and talonid) of each molar crown, to the nearest 1/100 mm, for 3283 specimens of 106 species. A total of 126 specimens were excluded from analysis due to either high levels of wear or absence of reporting wear condition. A further 24 specimens of indeterminate sex, 1 zoo specimen, and 221 specimens with incomplete lower molar dentitions were also excluded. To insure reasonably accurate estimates of species mean molar proportions and coefficients of variation, we included only species with at least four specimens that sampled both sexes. These exclusion criteria winnowed the original sample to 2,895 specimens from 100 taxa. Two subspecies of Pan troglodytes (P. t. schweinfurthii and P. t. troglodytes) were retained in the analysis. These anthropoid taxa represent a broad phylogenetic sample from Cercopithecoidea (n=61), Platyrrhini (n=25), and Hominoidea (n=14) (Table 1, Additional file 1: Table S1). We updated genus and species names to reflect the taxonomy of Wilson and Reeder [77].
For the analysis of M _{2}/M _{1} proportions when M _{3} exhibits agenesis, we collected lower molar area data from 66 modern human specimens across 6 populations (Additional file 1: Table S5) from the Museum of Comparative Zoology, Harvard. To confirm M _{3} agenesis, radiographs were taken using a SKYSEA™ dental portable radiography machine and Ergonom selfdeveloping Xray film. Each quadrant was radiographed with the occlusal plane at M _{2} perpendicular to the Xray. In addition, solely for the analysis of developmental timing, we obtained species mean molar proportion values of one strepsirrhine primate (Varecia variegata) from Swindler [59] (Additional file 1: Table S6).
To account for phylogenetic dependence among molar data during analysis, we used a majorityrule consensus tree of the 100 sampled anthropoid taxa, which was derived from the 10K Trees website and other sources (Additional file 1: Figure S12; [78, 79]). This phylogeny was generated from several sources of molecular data, including nuclear DNA, mitochondrial DNA, and Ychromosome sequences. We used a single phylogenetic tree in our analyses, rather than a Bayesian posterior distribution of trees, as multiple trees were not available for a proportion of the sampled taxa. A potential limitation of the study, therefore, is that phylogenetic uncertainty was not accounted for in our models [80].
To examine the relationships between primate molar proportions and various traits of interest, we compiled data on: a) the calcification of successive molar crowns (temporal overlap of M _{1} and M _{2}) from histological and radiographic sources (Additional file 1: Table S7), b) primary dietary category (fruit, leaves, insects, seeds, animals, omnivore) from the primatological literature (Additional file 1: Table S1), and c) body mass (Additional file 1: Table S1). Morphometric data, phylogenetic data (NEXUS format), and R code are archived in an online Zenodo data repository [40].
To test the efficacy of our molar area correction method, we collected molar measurements from three anthropoid species (Alouatta seniculus, Cercopithecus mitis, Homo sapiens) from the Museum of Comparative Zoology, Harvard. For each species, 10 wildshot specimens (5 male, 5 female) were selected. Only adult specimens were used, as determined by complete fusion of sphenooccipital synchrondrosis and epiphyses on the clavicle and pelvis [81, 82].
Data acquisition and preparation
We calculated molar crown areas for each anthropoid specimen using the product of buccolingual breadth (averaging trigonid and talonid buccolingual measurements) and mesiodistal length taken from Plavcan [76]. Since sample sizes for each sex were not always balanced, we calculated speciesspecific weighted means of molar crown areas by first aggregating to sexspecific species mean areas and then sexpooled species mean areas (Additional file 1: Table S1). We also calculated a smallsample coefficient of variation (c v) for each molar type and species combination.
For the polymorphic modern human sample (n=66), we extracted outline areas from scaled superior occlusal photographs using the Polygon tool in Image J (v. 1.45r) [83]. Specimens were photographed using a Canon Digital Rebel XS 10 megapixel digital SLR camera fitted with an EFS 60mm macro lens, following protocols established by GómezRobles and colleagues [84].
The product of linear breadth and length measurements often overestimates molar occlusal area (Additional file 1: Figure S1; [26]). We therefore used a correction to make these rectangular areas more comparable to areas calculated by tracing outlines around molar occlusal perimeters. Outline and rectangular areas from three species in different anthropoid clades, were extracted from scaled superior occlusal photographs using the Polygon and StraightSegment tools in Image J (v. 1.45r) [83]. Specimens were photographed using the same equipment and protocols described above.
We calculated areas from length (l) and breadth (b) variables to estimate each molar’s area as both a rectangle (r a=b l) and ellipse (e a=π(b l)/4) and then compared these estimates with the molar’s outline area (o a) to solve for a coefficient (x) of molar shape: x=o a−e a/r a−e a. We then created corrected areas (c a) by applying these coefficients, solving for area: c a=r a(x)+e a(1−x). We tested the usefulness of our correction method by comparing the difference between outline areas and corrected areas to the difference between outline areas and rectangular areas using a GLMM.
Selection of methods
Several studies investigating DIC model predictions with mammalian molar areas have used a reduced (standardized) major axis estimator (RMA; e.g., [16, 22, 26, 28, 29]), which assumes both Y and X are random variables measured with error and thus seeks to minimize error orthogonal to the model. We chose not to use RMA, since this method is known to produce large overestimates of slope when the ratio of error variances between Y and X is not unity [85, 86] or when the ratio of error variances does not equal the ratio of variances [87]. With empirical data, both these assumptions are unlikely to be met [88]. In addition, specieslevel observations are correlated with phylogeny [89] and therefore violate the assumption of statistical independence that is fundamental to many regression models, including RMA. Accounting for phylogenetic dependence in the data is possible, but complicated in a RMA framework. Currently the only software implementation is for the simple case of bivariate Gaussian data. Notably, none of the above studies accounted for phylogenetic dependence in their RMA models.
One alternative regression method some studies have used (e.g., [29]) is based on a generalized least squares estimator (PGLS) with an additional parameter (λ) that scales offdiagonal elements of the variancecovariance matrix in the error term, according to an hypothesis of phylogenetic relationships. This method, in common with the ordinary least squares estimator, treats explanatory variables as fixed and therefore minimizes error only in the Y variable. In addition, it is currently not possible to fit models that incorporate individuallevel data, and thus measurement error, using PGLS.
Due to the shortcomings of the above methods, we elected to use Bayesian phylogenetic generalized linear mixed models, estimated using a Markov chain Monte Carlo approach. It is straightforward to incorporate phylogenetic information into PGLMMs to account for interspecies dependence [90]. In addition, the PGLMM estimator has the advantage of being able to model individuallevel data, thus allowing intraspecific variance (both polymorphism and measurement error) to be accounted for in estimates of amongspecies relationships [80, 91]. Variation below the species level is considered a source of uncertainty (error) in specificlevel variables and, if ignored, can produce biased estimates of amongspecies relationships and inflated type I error rates [92–96]. Incorporating specificlevel error in models often decreases the precision with which parameters can be estimated, but provides more realistic interval estimates for those parameters. We fitted all models using the MCMC glmm() function from package MCMC glmm v. 2.21 [91] in R v. 3.1.3 [97].
Phylogenetic generalized linear mixed models
General model components and priors
The basic components of the Bayesian PGLMM are:
The random component (Eq. 3a) comprises a vector of observed trait values y and a probability distribution \(\mathcal {D}\) with mean μ and variance ϕ. For the logistic distribution ϕ is a constant (ϕ=π ^{2}/3), while for the Gaussian distribution it has to be estimated (ϕ=σ ^{2}). The systematic component (Eq. 3d) contains the linear part of the model, with the linear predictor η composed of a design matrix W that relates the explanatory variables (w _{1}…w _{ n }) to the data and a vector of location effects θ (‘fixed’ and ‘random’ effects). In between these two parts, a hierarchical layer is added to the model (Eq. 3c), where l is a hypothetical latent variable composed of the linear predictor η and a vector of residuals e (or the effect due to additive dispersion for nonGaussian models). The final component (Eq. 3b) joins the random and systematic parts of the model together using a link function g(·), or more usually its inverse g^{−1}(·). This function is used to relate the latent variable l to the observed data y, by transforming l into a quantity g^{−1}(l) that is the expectation of the distribution of y. For the Bernoulli distribution this would be the logistic function μ=logit^{−1}(l)=logistic(l) while for the Gaussian distribution it is the identity function μ=l.
In our analyses, when the response variable was Gaussian and an identity link function specified, then y was assumed to be normally distributed (\(\mathcal {N}\)) with expectation equal to the latent variable:
When the response variable was binary and a logit link function specified, then y was assumed to be Bernoulli distributed (\(\mathcal {B}\)) with expectation equal to the logistically transformed latent variable:
In a Bayesian analysis there is no distinction between fixed and random effects, with all effects treated as random [98]. However, since these terms are prevalent in the literature we note that W can be further deconstructed into design matrices for ‘fixed’ (X) and ‘random’ (Z) effects, and θ can be deconstructed into vectors of ‘fixed’ (β) and ‘random’ (γ) effects parameters:
The ‘fixed’ effects were assumed a priori to be independently distributed with zero mean and specified variance (\({\sigma ^{2}_{b}}\)) or scale (γ _{ b }):
where I is an identity matrix. To represent diffuse prior knowledge of the ‘fixed’ effects, for Gaussian response models (Formula 6a) we used a normal distribution and set \({\sigma ^{2}_{b}}\) to be 1×10^{8}, while for categorical responses (Formula 6b) we used a Cauchy distribution (\(\mathcal {C}\)) with γ _{ b } equal to π ^{2}/3+v, where v is the total variance of the ‘random’ and residual effects. This Cauchy prior is approximately flat on the probability scale [99].
In our models, the ‘random’ effects design matrix was composed of:
where p is a phylogenetic (correlated) random variable, s is a speciesspecific (i.i.d.) random variable, and γ _{ i } are vectors of ‘random’ effects parameters. The above terms accounted for betweenspecies variation in intercepts due to both phylogenetic and multiple measurement effects [100]. It is also possible to model amongspecies variation in slopes, but this requires a large intraspecific sample size for all species [101]. Since interspecies slope heterogeneity in our sample of anthropoids was low (Additional file 1: Figure S13), and sample sizes for some species were small (Additional file 1: Table S1), we did not increase the complexity of our models by adding terms for random slopes.
The ‘random’ effects were also assumed to be normally distributed about zero, but with variances (\({\sigma ^{2}_{p}}\) and \({\sigma ^{2}_{s}}\)) inferred a posteriori:
where Σ is a scaled matrix of phylogenetic covariances calculated from an ultrametric majorityrule consensus tree of the sampled anthropoid taxa [80]. Including the p term in the model thus makes the assumption that phylogenetic effects are correlated according to Σ. To embody our defuse prior knowledge of the ‘random’ effects we specified parameter expanded priors for p and s, with scale = 1, degree of belief = 1, mean = 0, and variance =1×10^{3} [91]. It has been suggested that parameter expanded priors are preferable when weakly informative priors are sought [102].
The error term (e) was also assumed to be normally distributed:
where the residual variance \({\sigma ^{2}_{e}}\) was estimated for Gaussian responses and fixed at π ^{2}/3 for Bernoulli responses.
Phylogenetic signal
Phylogenetic signal for each model was calculated using Pagel’s [41] λ parameter:
where \({\sigma ^{2}_{e}}\) was estimated for Gaussian responses and fixed at π ^{2}/3 for categorical responses. The λ parameter is measured on the interval [0, 1] and multiplies offdiagonal elements of Σ to reflect the pattern of observed covariance in trait values under a Brownian motion model.
Withingroup centering for specieslevel effects
When multiple measurements for each species are included in regression models the resulting β parameters measure a combination of between and withinspecies effects. To disentangle inter from intraspecific effects we therefore used withingroup centering [100, 103]. This method partitioned each predictor x into two components, one containing the specific mean of x, the other containing a measure of withinspecies variability (deviation of individual measurements from the specific mean). For individual j belonging to species i the generic model was thus:
where:
with J _{ i } being the number of individuals in species i. For each predictor x, there were thus two slopes: β _{B} the betweenspecies slope and β _{W} the (common) slope within each species [100, 103].
MCMC sampling and model convergence criteria
For Gaussian response models we sampled from the posterior distribution using MCMC for 1.1×10^{7} iterations, with a burnin period of 1×10^{6} and thinning interval of 1000. For Bernoulli response models the posterior MCMC sample was run for 1.5×10^{7} iterations, with a burnin period of 3×10^{6} and thinning interval of 1200. We determined that models had converged to a stationary posterior distribution when Hiedelberger and Welch’s diagnostic test [42] was passed and autocorrelation levels dipped below 0.1.
Data analysis
Molar area proportion relationships
To test the first prediction of the DIC model we fitted a regression model of molar proportions. The relationship between M _{3}/M _{1} and M _{2}/M _{1} occlusal area proportions can be described by the following model:
where the parameters β _{0} and β _{1B} are the betweenspecies intercept and slope of a line representing the population average of M _{3}/M _{1} conditional on a given value of M _{2}/M _{1}, and β _{1W} is the (pooled) withinspecies slope.
M _{2} relative area
To determine if M _{2} area accounts for onethird of lower molar total area, we regressed the sum of lower molar areas on M _{2} area using the following model:
where M _{2} is the occlusal surface area of the second lower molar, M _{T} is the sum of all three lower molar areas, and β _{1B} and β _{1W} are between and (pooled) withinspecies estimates of the proportion of total lower molar area accounted for by M _{2}.
Molar area relationships
We assessed how much the relationship between M _{1} and M _{3} area was mediated by M _{2} area, using two models (Eqs. 14a, 14b) to estimate parameters of the three pathways representing molar area relationships:
where M _{ i } is the occlusal surface area and β _{ iB} and β _{ iW} are the between and (pooled) withinspecies slopes of the i ^{th} lower molar. The indirect effect was calculated using Eq. 14c, where \(\beta _{2\mathrm {B}}^{(t)}\) and \(\beta _{1\mathrm {B}}^{\prime (t)}\) are the t ^{th} parameter estimates for t=1,…,T MCMC draws from the posterior distribution. The direct effect was calculated in the same manner using Eq. 14d. Point estimates of the mediated (\(\widehat {ab}\)) and direct (\(\lvert \widehat {c^{\prime }} \rvert \)) effects are the mean of these individual draws, while the 95 % HDI is given by the (.025,.975) sample quantiles of the posterior draws [48]. We calculated the proportion mediated \(\widehat {\mathit {pr_{m}}}\) as the ratio of indirect effect to total effect, with the latter being the sum of the indirect effect and the absolute value of the direct effect (Eq. 14e) [48, 104].
Body mass and fit/deviation from the DIC model
We used the following model to test the association between fit/deviation from the DIC model line and body mass:
where \(\text {ln}(\texttt {mass}_{i}^{\text {kg}})\) is the natural log of species mean body mass and \(\lvert \texttt {dic}_{\mathit {ij}}^{\texttt {dev}} \rvert \) is the absolute perpendicular distance from the DIC model line to each individual’s location in molar proportion space.
Molar area variability
We tested whether average third molar area variability (\(M_{3}^{\mathit {cv}}\)) across all anthropoid species sampled was larger than either average first (\(M_{1}^{\mathit {cv}}\)) or second (\(M_{2}^{\mathit {cv}}\)) molar area variability using the following model:
where:
and
and where \(\widehat {\mathit {cv}}\) is the sample estimate of the coefficient of variation for small sample sizes:
We used \(M_{3}^{\mathit {cv}}\) as the reference level for the indicator variables.
M _{3} agenesis
We fitted the following model to determine whether M _{3} agenesis across our anthropoid sample is associated with values of M _{2}/M _{1}<0.5:
where:
and
We also fitted a model to determine whether M _{3} agenesis within a polymorphic modern human sample is associated with values of M _{2}/M _{1}<0.5:
where r is a populationspecific (i.i.d.) random variable with i levels, distributed as:
and where:
and
Primary dietary category and molar area proportions
We used the following logistic regression model to test the association between primary dietary category and the relative size of M _{2}:
where:
This equation models the probability of being in the M _{1}<M _{2}>M _{3} region of molar proportion morphospace as a function of a species’ primary dietary category.
Ethics and consent to participate
Not applicable.
Consent to publish
Not applicable.
Availability of data and materials
The morphometric and phylogenetic data supporting the results of this article, together with an R script that replicates the analyses, are available in the Zenodo repository (http://dx.doi.org/10.5281/zenodo.50108).
Abbreviations
 c v :

coefficient of variation for small sample sizes
 dic :

dental inhibitory cascade
 hdi :

highest density interval
 M _{1} :

first lower molar
 M _{2} :

second lower molar
 M _{3} :

third lower molar
 \(\mathbb {P}\) :

posterior probability
 rope :

region of practical equivalence
References
 1
Hunter JP, Jernvall J. The hypocone as a key innovation in mammalian evolution. Proc Natl Acad Sci USA. 1995; 92:10718–10722.
 2
Mendoza M, Palmqvist P. Hypsodonty in ungulates: an adaptation for grass consumption or for foraging in open habitat?J Zool. 2008; 274:134–42.
 3
Jardine PE, Janis CM, Sahney S, Benton MJ. Grit not grass: concordant patterns of early origin of hypsodonty in Great Plains ungulates and Glires. Palaeogeogr Palaeoclimatol Palaeoecol. 2012; 365:1–10.
 4
von Koenigswald W. Evolutionary trends in the enamel of rodent incisors In: Luckett WP, Hartenberger JL, editors. Evolutionary Relationships Among Rodents: Comments and Conclusions. New York: Springer: 1985. p. 403–22.
 5
Williams SH, Kay RF. A comparative test of adaptive explanations for hypsodonty in ungulates and rodents. J Mammal Evol. 2001; 8:207–29.
 6
Luo ZX, Cifelli RL, KielanJaworowska Z. Dual origin of tribosphenic mammals. Nature. 2001; 409:53–7.
 7
Luo ZX, Qiang J, Yuan CX. Convergent dental adaptations in pseudotribosphenic and tribosphenic mammals. Nature. 2007; 450:93–7.
 8
Woodburne MO, Rich TH, Springer MS. The evolution of tribospheny and the antiquity of mammalian clades. Mol Phylo Evol. 2003; 28:360–85.
 9
Werdeline L. Supernumerary teeth in Lynx lynx and the irreversibility of evolution. J Zool. 1987; 2:259–66.
 10
Wiens JJ. Reevolution of lost mandibular teeth in frogs after more than 200 million years and reevaluating dollo’s law. Evol. 2011; 65:1283–1296.
 11
Godinot M. Primate origins: a reappraisal of historical data favoring tupaiid affinities In: Ravosa MJ, Dagosto M, editors. Primate Origins: Adaptations and Evolution. New York: Springer: 2007. p. 83–142.
 12
Gingerich PD. Homologies of the anterior teeth in Indriidae and a functional basis for dental reduction in primates. Am J Phys Anthropol. 1977; 47:387–93.
 13
Rose KD. The Beginning of the Age of Mammals. Baltimore: JHU Press; 2006.
 14
Jernvall J, Fortelius M. Common mammals drive the evolutionary increase of hypsodonty in the Neogene. Nature. 2002; 417:538–40.
 15
Renaud S, Auffray JC, Michaux J. Conserved phenotypic variation patterns, evolution along lines of least resistance, and departure due to selection in fossil rodents. Evol. 2006; 60(8):1701–1717.
 16
Kavanagh KD, Evans AR, Jernvall J. Predicting evolutionary patterns of mammalian teeth from development. Nature. 2007; 449:427–33.
 17
Munne PM, Felszeghy S, Jussila M, Suomalainen M, Thesleff I, Jernvall J. Splitting placodes: effects of bone morphogenetic protein and Activin on the pattering and identity of mouse incisors. Evol Dev. 2009; 12:383–92.
 18
Nieminen P. Genetic basis of tooth agenesis. J Exp Zool B. 2009; 312:320–42.
 19
Prochazka J, Pantalacci S, Churava S, Rothova M, Lambert A, Lesot H, Klein O, Peterka M, Laudet V, Pterkova R. Patterning by heritage in mouse molar row development. Proc Natl Acad Sci USA. 2010; 107:15497–502.
 20
Lucas PW, Corlett RT, Luke DA. Postcanine tooth size and diet in anthropoid primates. Z Morphol Anthropol. 1986; 76:253–76.
 21
Lucas PW, Teaford MF. Functional morphology of colobine teeth In: Davies G, Oates J, editors. Colobine Monkeys: Their Ecology, Behaviour and Evolution. Cambridge: Cambridge University Press: 1994. p. 173–203.
 22
Wilson LAB, SánchezVillagra MR, Madden RH, Kay RF. Testing a developmental model in the fossil record: molar proportions in South American ungulates. Paleobiol. 2012; 38:308–21.
 23
Halliday TJD, Goswami A. Testing the inhibitory cascade model in Mesozoic and Cenozoic mammaliaforms. BMC Evol Biol. 2013; 13:79–90.
 24
Polly PD. Evolutionary biology: development with a bite. Nature. 2007; 449:413–5.
 25
Laffont R, Renvoise E, Navarro N, Alibert P, Montuire S. Morphological modularity and assessment of developmental processes within the vole dental row (Microtus arvalis, Arvicolinae, Rodentia). Evol Dev. 2009; 11:302–11.
 26
Renvoisé E, Evans AR, Jebrane A, Labruere C, Laffont R, Montuire S. Evolution of mammal tooth patterns: new insights from a developmental prediction model. Evol. 2009; 63:1327–40.
 27
Labonne G, Laffont R, Renvoisé E, Jebrane A, Labruère C, ChateauSmith C, Navarro N, Montuire S. When less means more: evolutionary and developmental hypotheses in rodent molars. J Evol Biol. 2012; 25:2102–111.
 28
Asahara M. Unique inhibitory cascade pattern of molars in canids contributing to their potential to evolutionary plasticity of diet. Ecol Evol. 2013; 3(2):278–85.
 29
Bernal V, Gonzalez PN, Perez SI. Developmental processes, evolvability and dental diversification of New World monkeys. Evol Biol. 2013; 40:532–41.
 30
Pouladi MA, Morton AJ, Hayden MR. Choosing an animal model for the study of Huntington’s disease. Nat Rev Neurosci. 2013; 14(10):708–21.
 31
Peterkova R, Hovorakova M, Peterka M, Lesot H. Threedimensional analysis of the early development of the dentition. Australian Dent J. 2014; 59(S1):55–80.
 32
Scott JE. Lost and found: the third molars of Callimico goeldii and the evolution of the callitrichine postcanine dentition. J Hum Evol. 2015; 83:65–73.
 33
Simons EL. Parapithecus grangeri of the African Oligocene: an archaic catarrhine without lower incisors. J Hum Evol. 1986; 15:205–13.
 34
Andrews P, Martin L. Cladistic relationships of extant and fossil hominoids. J Hum Evol. 1987; 16:101–18.
 35
Miller ER, Benefit BR, McCrossin ML, Plavcan JM, Leakey MG, ElBarkooky AN, Hamdan MA, Abdel Gawad MK, Hassan SM, Simons EL. Systematics of early and middle Miocene Old World monkeys. J Hum Evol. 2009; 57:195–211.
 36
Garn SM, Lewis AB, Kerewsky RS. Relative molar size and fossil taxonomy. Am Anthropol. 1964; 66:587–92.
 37
Garn SM, Lewis AB. The gradient and pattern of crownsize reduction in simple hypodontia. Angle Orthod. 1970; 40:51–8.
 38
Sofaer JA. A model relating developmental interaction and differential evolutionary reduction of tooth size. Evol. 1973; 27:427–34.
 39
Ford SM. Callitrichids as phyletic dwarfs, and the place of the Callitrichidae in Platyrrhini. Primates. 1980; 21:31–43.
 40
Worthington S. Anthropoid morphometric and phylogenetic data, with R replication code. 2016. Zenodo Repository, http://dx.doi.org/10.5281/zenodo.50108.
 41
Pagel M. Inferring the historical patterns of biological evolution. Nature. 1999; 401:877–84.
 42
Hiedelberger P, Welch PD. Simulation run length control in the presence of an initial transient. Opns Res. 1983; 31:1109–1144.
 43
Kruschke JK. Bayesian assessment of null values via parameter estimation and model comparison. Perspect Psychol Sci. 2011; 6(3):299–312.
 44
Kruschke JK. Bayesian estimation supersedes the t test. J Exp Psychol Gen. 2013; 142(2):573–603.
 45
Jasieński M, Bazzaz FA. The fallacy of ratios and the testability of models in biology. Oikos. 1999; 84(2):321–6.
 46
Lucas PW, Corlett RT, Luke DA. Sexual dimorphism of tooth size in anthropoids. Hum Evol. 1986; 1:23–9.
 47
Sobel ME. Asymptotic confidence intervals for direct effects in structural equation models In: Leinhardt S, editor. Sociological Methodology. Washington DC: American Sociological Association: 1982. p. 290–312.
 48
Yuan Y, MacKinnon DP. Bayesian mediation analysis. Psychol Methods. 2009; 14(4):301–22.
 49
Butler PM. Studies of the mammalian dentition. Differentiation of the postcanine teeth. Proc Zool Soc Lond B. 1939; B109:1–36.
 50
Carter KE, Worthington S. Morphologic and demographic predictors of third molar agenesis: a systematic review and metaanalysis. J Dent Res. 2015; 94(7):886–94.
 51
Evans AR, Daly ES, Catlett KK, Paul KS, King SJ, Skinner MM, Nesse HP, Hublin JJ, Townsend GC, Schwartz GT, Jernvall J. A simple rule governs the evolution and development of hominin tooth size. Nature. 2016; 530:477–81.
 52
Carter KE, Worthington S, Smith TM. The evolution of third molar agenesis in humans. Am J Phys Anthropol. 2013; 150(Suppl. 56):95–5.
 53
MacPhee RDE, Horovitz I. New craniodental remains of the Quaternary Jamaican monkey Xenothrix mcgregori (Xenotrichini, Callicebinae, Pitheciidae), with a reconsideration of the Aotus hypothesis. Amer Mus Nov. 2004; 3434:1–51.
 54
Cai J, Cho SW, Kim JY, Lee MJ, Cha YG, Jung HS. Patterning the size and number of tooth and its cusps. Dev Biol. 2007; 304(2):499–507.
 55
Gingerich PD, Smith BH, Rosenberg KR. Allometric scaling in the dentition of primates and prediction of body weight from tooth size in fossils. Am J Phys Anthropol. 1982; 58:81–100.
 56
Koyabu DB, Endo H. Craniodental mechanics and diet in Asian colobines: morphological evidence of mature seed predation and sclerocarpy. Am J Phys Anthropol. 2010; 142:137–48.
 57
Walker P, Murray P. An assessment of masticatory efficiency in a series of anthropoid primates with special reference to the Colobinae and Cercopithecinae In: Tuttle RH, editor. Primate Functional Morphology and Evolution. The Hague: Mouton: 1975. p. 135–50.
 58
Jheon AH, Seidel K, Biehs B, Klein OD. From molecules to mastication: the development and evolution of teeth. Dev Biol. 2012; 2:165–82.
 59
Swindler DR. Primate Dentition: An Introduction to the Teeth of Nonhuman Primates. Cambridge: Cambridge University Press; 2002.
 60
Tummers M, Thesleff I. The importance of signal pathway modulation in all aspects of tooth development. J Exp Zool B Mol Dev Evol. 2009; 312(4):309–19.
 61
Fraser GJ, Hulsey CD, Bloomquist RF, Uyesugi K, Manley NR, Streelman JT. An ancient gene network is coopted for teeth on old and new jaws. PLoS Biol. 2009; 7:1000031.
 62
Fraser GJ, Smith MM. Evolution of developmental pattern for vertebrate dentitions: an oropharyngeal specific mechanism. J Exp Zool B. 2011; 316:99–112.
 63
Buchtová M, Handrigan GR, Tucker AS, Lozanoff S, Town L, Fu K, Diewert VM, Wicking C, Richman JM. Initiation and patterning of the snake dentition are dependent on Sonic hedgehog signaling. Dev Biol. 2008; 319(1):132–45.
 64
Jernvall J. Linking development with generation of novelty in mammalian teeth. Proc Natl Acad Sci USA. 2000; 97:2641–645.
 65
Streelman JT, Albertson RC. Evolution of novelty in the cichlid dentition. J Exp Zool B Mol Dev Evol. 2006; 306:216–26.
 66
Kangas AT, Evans AR, Thesleff I, Jernvall J. Nonindependence of mammalian dental characters. Nature. 2004; 432:211–4.
 67
Hlusko LJ, Mahaney MC. Quantitative genetics, pleiotropy, and morphological integration in the dentition of Papio hamadryas. Evol Biol. 2009; 36(1):5–18.
 68
Hlusko LJ, Sage RD, Mahaney MC. Modularity in the mammalian dentition: mice and monkeys share a common dental genetic architecture. J Exp Zool B. 2011; 316:21–49.
 69
Ahn Y, Sanderson BW, Klein OD, Krumlauf R. Inhibition of Wnt signaling by Wise (Sostdc1) and negative feedback from Shh controls tooth number and patterning. Development. 2010; 137:3221–231.
 70
Cho SW, Kwak S, Woolley TE, Lee MJ, Kim EJ, Baker RE, Kim HJ, Shin JS, Tickle C, Maini PK, Jung HS. Interactions between Shh, Sostdc1 and Wnt signaling and a new feedback loop for spatial patterning of the teeth. Development. 2011; 138(9):1807–16.
 71
Wu N, Iwamoto T, Sugawara Y, Futaki M, Yoshizaki K, Yamamoto S, Yamada A, Nakamura T, Nonaka K, Fukumoto S. PDGFs regulate tooth germ proliferation and ameloblast differentiation. Arch Oral Biol. 2010; 55(6):426–34.
 72
Gilbert SF. Evolution through developmental changes: how alterations in development cause evolutionary changes in anatomy In: Auletta G, Leclerk M, Martínez RA, editors. Biological Evolution: Facts and Theories. Roma: Gregorian And Biblical Press: 2011. p. 153–67.
 73
Oerke AK, Hesitermann M, Kuderlin U, Martin RD, Hodges JK. Monitoring reproduction in Callitrichidae by means of ultrasonography. Evol Anthropol. 2002; 11(Suppl. 1):183–5.
 74
Harjunmaa E, Seidel K, Hakkinen T, Renvoise E, Corfe IJ, Kallonen A, Zhang ZQ, Evans AR, Mikkola ML, SalazarCiudad I, Klein OD, Jernvall J. Replaying evolutionary transitions from the dental fossil record. Nature. 2014; 512(7512):44–8.
 75
Schwartz GT, Samonds KE, Godfrey LR, Jungers WL, Simons EL. Dental microstructure and life history in subfossil Malagasy lemurs. Proc Natl Acad Sci USA. 2002; 99:6124–129.
 76
Plavcan JM. 1990. Sexual Dimorphism in the Dentition of Extant Anthropoid Primates. Ph.D. Dissertation. University of Michigan, Ann Arbor, MI.
 77
Wilson DE, Reeder DM, editors. Mammal Species of the World. A Taxonomic and Geographic Reference, 3rd edn. Baltimore, MD: John Hopkins University Press; 2005.
 78
Arnold C, Matthews LJ, Nunn CL. The 10kTrees website: a new online resource for primate phylogeny. Evol Anthropol. 2010; 19(3):114–8.
 79
Meyer D, Rinaldi ID, Ramlee H, PerwitasariFarajallah D, Hodges JK, Roos C. Mitochondrial phylogeny of leaf monkeys (genus Presbytis, Eschscholtz, 1821) with implications for taxonomy and conservation. Mol Phylo Evol. 2011; 59:311–9.
 80
de Villemereuil P, Wells JA, Edwards RD, Blomberg SP. Bayesian models for comparative analysis integrating phylogenetic uncertainty. BMC Evol Biol. 2012; 12:102–18.
 81
Molleson TI. The Biology of Human Ageing In: Bittles AH, Collins KJ, editors. Cambridge: Cambridge University Press: 1986. p. 95–118.
 82
Kreitner KF, Schweden FJ, Riepert T, Nafe B, Thelen M. Bone age determination based on the study of the medial extremity of the clavicle. Eur Radiol. 1998; 8:1116–22.
 83
Abràmoff DMD, Magalhães DPJ, Ram DSJ. Image processing with ImageJ. Biophoton Int. 2004; 11(7):36–42.
 84
GómezRobles A, MartinónTorres M, Bermúdez de Castro JM, Margvelashvili A, Bastir M, Arsuaga JL, PérezPérez A, Estebaranz F, Martínez LM. A geometric morphometric analysis of hominin upper first molar shape. J Hum Evol. 2007; 53:272–85.
 85
Carroll RJ, Ruppert D. The use and misuse of orthogonal regression in linear errorsinvariables models. J Amer Stat. 1996; 50(1):1–6.
 86
Smith RJ. Use and misuse of the reduced major axis for linefitting. Am J Phys Anthropol. 2009; 140:476–86.
 87
Warton DI, Wright IJ, Falster DS, Westoby M. Bivariate linefitting methods for allometry. Biol Rev. 2006; 81(2):259–91.
 88
Carroll RJ. Measurement Error in Nonlinear Models: a Modern Perspective. Monographs on Statistics and Applied Probability. Boca Raton: Chapman & Hall/CRC; 2006.
 89
Felsenstein J. Phylogenies and the comparative method. Am Nat. 1985; 125(1):1–15.
 90
Hadfield JD, Nakagawa S. General quantitative genetic methods for comparative biology: phylogenies, taxonomies, metaanalysis and multitrait models for continuous and categorical characters. J Evol Biol. 2010; 23(3):494–508.
 91
Hadfield JD. mcmc methods for multiresponse generalised linear mixed models: the MCMCglmm R package. J Stat Soft. 2010; 33(2):1–22.
 92
Harmon LJ, Losos JB. The effect of intraspecific sample size on type I and type II error rates in comparative studies. Evol. 2005; 59:2705–710.
 93
Ives AR, Midford PE, Garland T. Withinspecies variation and measurement error in phylogenetic comparative methods. Syst Biol. 2007; 56(2):252–70.
 94
Felsenstein J. Comparative methods with sampling error and withinspecies variation: contrasts revisited and revised. Am Nat. 2008; 171(6):713–25.
 95
Garamszegi LZ, Møller AP. Effects of sample size and intraspecific variation in phylogenetic comparative studies: a metaanalytic review. Biol Rev. 2010; 85:797–805.
 96
Garamszegi LZ. Uncertainties due to withinspecies variation in comparative studies: measurement errors and statistical weights In: Garamszegi LZ, editor. Modern Phylogenetic Comparative Methods and Their Application in Evolutionary Biology. Berlin: Springer: 2014. p. 157–99.
 97
R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2015. R Foundation for Statistical Computing. http://www.Rproject.org/.
 98
Gelman A, Carlin JB, Stern HS, Rubin DB. Bayesian Data Analysis, 3rd edn. Boca Raton: CRC Press; 2013.
 99
Gelman A, Jakulin A, Pittau MG, Su YS. A weakly informative default prior distribution for logistic and other regression models. Ann Appl Stat. 2008; 2(4):1360–83.
 100
de Villemereuil P, Nakagawa S. General quantitative genetic methods for comparative biology In: Garamszegi LZ, editor. Modern Phylogenetic Comparative Methods and Their Application in Evolutionary Biology. Berlin: Springer: 2014. p. 287–303.
 101
Phillimore AB, Hadfield JD, Jones OR, Smithers RJ. Differences in spawning date between populations of common frog reveal local adaptation. Proc Natl Acad Sci USA. 2010; 107(18):8292–297.
 102
Gelman A. Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper). Bayesian Anal. 2006; 1(3):515–34.
 103
van de Pol M, Wright J. A simple method for distinguishing within versus betweensubject effects using mixed models. Anim Behav. 2009; 77(3):753–8.
 104
Alwin DF, Hauser RM. The decomposition of effects in path analysis. Am Sociol Rev. 1975; 40:37–47.
Acknowledgements
We thank Michael Plavcan for providing dental measurement data. We also thank Tanya M. Smith, David Pilbeam, and two anonymous reviewers for helpful discussion and constructive comments on earlier versions of the manuscript. We are grateful to Judith Chupasko for providing access to the primate skeletal material in her care at the Museum of Comparative Zoology, Harvard.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
KC conceived the study, participated in the study design, collected modern human molar data from the Museum of Comparative Zoology at Harvard, and contributed to writing the manuscript. SW participated in the study design, contributed to writing the manuscript, wrote the R code, conducted the statistical analyses, and generated the figures and tables. Both authors read and approved the final manuscript.
Authors’ information
KC conducted this research as part of her Ph.D. degree and is broadly interested in the evolution and development of mammalian teeth. SW is an evolutionary biologist specializing in statistical modeling, morphometrics, comparative and phylogenetic methods.
Funding
This study was funded by a National Science Foundation GRFP grant (DGE1144152) to K. Carter.
Additional file
Additional file 1
Table S1. Anthropoid (n=100) lower molar proportions, body mass, and diet summary data. Table S2. Molar area proportion pglmm posterior probabilities. Table S3. Relative M _{2} area pglmm posterior probabilities. Table S4. Proportion mediated (p r _{ m }) pglmm posterior probabilities. Table S5. Modern human (n=66) lower molar proportions data. Table S6. Strepsirrhine (n=1) lower molar proportions summary data. Table S7. Primate molar crown formation time data. Figure S1. Molar area correction coefficients method. Figure S2. Molar area correction coefficients test. Figure S3. Anthropoid lower molar proportions of individual specimens (n=2,895). Figure S4. Posterior density of molar proportion slope (interspecific). Figure S5. Posterior density of molar proportion intercept. Figure S6. Posterior density of relative M _{2} area slope (interspecific). Figure S7. Posterior density of relative M _{2} area intercept. Figure S8. Posterior density of proportion mediated (interspecific). Figure S9. Ln body mass as a function of fit/deviation from the dic molar proportion model line. Figure S10. Anthropoid lower molar proportions intraspecific variation. Figure S11. Anthropoid lower molar proportions with respect to relative M _{2} area. Figure S12. Molecular phylogeny of anthropoid primates (n=100). Figure S13. Anthropoid lower molar proportions intraspecific error. Supplementary references. (1–43). (PDF 3891 kb).
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Molar proportions
 Dental inhibitory cascade
 Anthropoid primates
 Evodevo
 Bayesian phylogenetic generalized linear mixed models