Significance Experimental research on mice has yielded tremendous biological insight. However, the ∼140 million y of evolution that separate mice from humans pose a hurdle to direct application of this knowledge to humans. We report here that considerable progress for identifying genetically patterned skeletal phenotypes beyond the mouse model is possible through transdisciplinary approaches that include the anatomical sciences. Indeed, anatomy and paleontology offer unique opportunities through which to develop and test hypotheses about the underlying genetic mechanisms of the skeleton for taxa that are not well suited to experimental manipulation, such as ourselves.

Abstract Developmental genetics research on mice provides a relatively sound understanding of the genes necessary and sufficient to make mammalian teeth. However, mouse dentitions are highly derived compared with human dentitions, complicating the application of these insights to human biology. We used quantitative genetic analyses of data from living nonhuman primates and extensive osteological and paleontological collections to refine our assessment of dental phenotypes so that they better represent how the underlying genetic mechanisms actually influence anatomical variation. We identify ratios that better characterize the output of two dental genetic patterning mechanisms for primate dentitions. These two newly defined phenotypes are heritable with no measurable pleiotropic effects. When we consider how these two phenotypes vary across neontological and paleontological datasets, we find that the major Middle Miocene taxonomic shift in primate diversity is characterized by a shift in these two genetic outputs. Our results build on the mouse model by combining quantitative genetics and paleontology, and thereby elucidate how genetic mechanisms likely underlie major events in primate evolution.

The relationship between genotype and phenotype is critical to evolutionary biology, because it influences how phenotypes respond to selective pressures and evolve (1, 2). Paleontologists have long sought to incorporate the etiology of the dental phenotype to inform on questions of environmental, dietary, and adaptive change over time (3⇓⇓–6). This research has been advanced significantly by the revolution in developmental genetics over the past few decades (7). Experimental research on mice has yielded tremendous biological insight (8). However, for human phenotypes, ranging from inflammation (9) to placentation (10), the limitations of the mouse model due to the ∼140 million years ago of evolution that have occurred since our last common ancestor ∼70 Ma (11) are starting to be recognized. Here, we demonstrate how to overcome the limitations of the mouse model’s application to the primate dentition by integrating research from quantitative genetics, neontology, and paleontology. The insights gained from this transdisciplinary approach have implications for all of these seemingly disparate subdisciplines of biology.

There was considerable excitement when an inhibitory cascade (IC) mechanism, inferred from experimental manipulation of mouse molar development, showed tremendous explanatory power for the taxonomic variation observed across murines (12). The concept was soon applied to other tripartite skeletal phenotypes (13). However, the extent of the predictive power of the IC for taxa with less similar (14) and more heterodont (15⇓⇓–18) dentitions suggests that either the IC mechanism does not characterize these other organisms’ dental patterning as specifically or the characterization of the mechanistic output (the phenotype) is not accurately assessed. Despite empirical evidence for how primate dental patterning differs from dental patterning of the mouse (19, 20), proponents of the mouse model continue to adhere tightly to it. Most recently, Evans et al. (21) report that the IC model from mice applied uncritically to hominids reveals that the genus Homo is distinct from other hominid genera in having smaller third molars because of a “simple rule.”

Although the developmental mechanisms underlying tooth organogenesis are evolutionarily conserved across mammals (19, 22⇓–24), given the highly derived dentitions of mice (which lack premolars, canines, and replacement teeth) and the vastly different life histories, the finer details of dental patterning are not the same (19). Previous quantitative genetic analyses indicate that the IC model may not work as well for primates as it does for mice because the former’s phenotype can comingle the effect of a genetic patterning mechanism(s) with systemically modulated somatic factors. For example, baboon molar buccolingual width is significantly influenced by the same genes that influence the animal’s crown–rump length (20). Given that crown–rump length is highly sexually dimorphic, tooth size measurements that incorporate buccolingual width are capturing the signal of the systemic effects of body size and sex in addition to the dental patterning mechanism(s) more so than would mesiodistal length alone. Furthermore, we find that the lengths of the fourth premolar and the first through third molars have higher genetic correlations with one another than do the width measurements of these same teeth (19). This pattern of genetic correlation suggests that among primates, variation in mandibular mesiodistal tooth measurements may more directly reflect the genetic patterning mechanisms determining relative tooth sizes than do width measurements. Although Evans et al. (21) interpret their results as indicative of simply the IC, our research (19, 20) shows that the phenotype they have captured for primates includes significant pleiotropic effects with body size. The genus Homo is not only characterized by smaller third molars compared with other hominid genera (21) but also by an ∼30% increase in body size (25). Is the genus Homo different from other hominid genera because of a simple developmental rule or because of all the various factors that lead to a dramatic increase in body size? Evans et al. (21) cannot and do not distinguish between these two very different evolutionary phenomena because their phenotype definition fails to use empirical evidence from primate genetics to modify what we have learned from mouse models.

To extend the mouse model to our own lineage, gene-forward approaches such as are used in mice and other model organisms are not practical. Therefore, the next phase of knowledge will come from phenotype-back methods. To demonstrate this approach, we combined the strength of three datasets—quantitative genetics, neontology, and paleontology—to test hypotheses about how to better define the output of dental patterning mechanisms in primates. If a phenotype definition corresponds directly to a patterning mechanism, it should meet four expectations: (i) the genetically patterned (GP) trait will be significantly heritable, with little to no covariate effects; (ii) it will capture variation that differentiates taxa, representing genetic variance exploited by selection; (iii) different GP traits will be uncorrelated within a species (but may be phenotypically correlated across closely related taxa because of their phylogenetic proximity); and (iv) the fossil record will reveal differing evolutionary histories for GP traits because the genetically independent traits have responded to selection and/or drift relatively independently.

We identify two dental phenotypes that meet these four criteria. These phenotypes characterize the output of two (as-yet-unknown) patterning mechanisms that influence the relative proportions of the various teeth within the two genetic submodules of the primate permanent mandibular postcanine dentition (26). The first GP phenotype is the molar module component (MMC). The MMC follows Kavanagh et al.’s proposed developmental IC model (12), but by using only mesiodistal length, we capture the GP output more directly by removing the same modulating, pleiotropic effects that confounded Evans et al.’s analysis (21): MMC = M 3 l / M 1 l . The second GP phenotype is the premolar/molar module (PMM), defined as PMM = M 2 l / P 4 l , where l is the mesiodistal length of the respective molar crown, and subscripts indicate the position of the tooth crown (M, molar; P, premolar). We use only the length of the second molar in the PMM because it is a stable predictor of the overall length of the molar series [for murines (12) and for Old World monkeys (OWMs) (SI Appendix, Fig. S1)].

Discussion The dramatic evolution of these two GP dental traits across the Old World higher primates shows a remarkable correlation with global climate and vegetation change as reflected in decreasing δ18O and δ13C values following the Middle Miocene Climatic Optimum (34, 35). Paleobotanical and faunal evidence suggests that the largely mosaic forests and woodlands of Oligocene Africa transitioned in the Middle Miocene to more heterogeneous landscapes (36, 37). This climate change further correlates with the long-recognized shift in species diversity between apes and OWMs (Fig. 3D). Before the Middle Miocene, apes were abundant and diverse across Africa and Eurasia (38), whereas OWMs were still rare. However, by the Pliocene, monkeys had experienced a taxonomic and geographic expansion across the Old World (39) and ape diversity had dramatically reduced to a handful of relict genera that persist today. The origins of the human lineage are embedded within this evolutionary transition (31, 40). Evolutionary changes in the two GP phenotypes introduced here appear to reflect how their two underlying genetic mechanisms responded to selective pressures acting on the primate postcanine dentition. In turn, these mechanisms may provide a broad and deep evolutionary explanation for primate adaptation during the Neogene. Although Miocene apes occupied the same PMM space that extant apes and all of the hominids demonstrate (Fig. 2C) during the Miocene, the ape MMC phenotype experienced a downward shift as the first and third molars became more equivalent in size. During the Plio-Pleistocene, the papionins, and even some of the larger fossil colobines, evolved MMC and PMM phenotypes that occupy the morphospace of the now-extinct Miocene apes. From the perspective of the patterning of the postcanine dentition, baboon morphology replaced the Miocene ape strategy, which may well underlie and/or reflect the demographic shifts hypothesized for this evolutionary competition (41) (i.e., during the later Miocene, hominoids increasingly invested in parental care of fewer offspring compared to the papionin reproductive strategy). Given that teeth develop and erupt at different times over a long ontogenetic period and that premolars replace deciduous teeth and molars do not (and the mechanisms relating them are unknown), the MMC and PMM may represent two different responses to selective pressures that vary, in part, on the ontogenetic timing of the pressure. Although the earliest written record of human anatomical studies is 3,000 y old, when cadaver dissection became acceptable in the second century A.D., anatomical sciences became the foundation of medical research (42). However, since the discovery of DNA’s structure and the synthesis of organic compounds from inorganic molecules (43, 44), biology’s focus has progressively shifted from the macroscopic toward the microscopic. The wealth of new insight gained about the genetic and cellular underpinnings of biology might be argued as justification for this relegation of the comparative anatomical sciences, and their associated museum collections, to the disciplinary sidelines. However, because developmental genetics today aims for a genotype/phenotype map that can inform us about our own species’ biology, the many millions of years of evolution that stand between H. sapiens and mice still pose a formidable hurdle (8, 9, 45). We report here that considerable progress for identifying GP skeletal phenotypes beyond the mouse model is possible through transdisciplinary approaches that include the anatomical sciences. Indeed, anatomy and paleontology offer unique opportunities through which to develop and test hypotheses about the underlying genetic mechanisms of the skeleton for taxa that are not well suited to experimental manipulation, such as ourselves.

Materials and Methods Quantitative Genetics Analyses. We collected linear measurements of the mandibular postcanine teeth for 632 baboons that are part of a captive, pedigreed breeding colony of Papio hamadryas (46) housed at the Southwest National Primate Research Center (SNPRC) in San Antonio, Texas. Genetic management of the colony allows for data collection from animals that are not inbred. All nonfounder animals in this study resulted from matings that were random with respect to dental, skeletal, and developmental phenotype. The female-to-male sex ratio is ∼2:1. The animals from which linear tooth size measurements were collected [described in detail elsewhere (47)] are distributed across 11 extended pedigrees that are three to five generations deep. The mean number of animals with data per pedigree was 44, and these individuals typically occupied the lower two or three generations of each pedigree. All pedigree data management and preparation was facilitated through use of the computer package PEDSYS (48). The SNPRC Institutional Animal Care and Use Committee, in accordance with the established guidelines (49), approved all procedures related to the treatment of the baboons during the conduct of this study. Measurements used in these analyses are from the left mandible, except for the measurements of the MMC, which are taken from the right due to sample size limitations on the left. Tooth area is defined as the mesiodistal length of the tooth multiplied by the mesial buccolingual width, where applicable. All trait values were inverse-normalized before analysis to account for residual kurtosis. Statistical Genetic Analyses. We conducted our genetic analyses within a maximum likelihood framework, the standard approach in quantitative genetics for finding parameter estimates that maximize the likelihood of the genetic model, conditional on all available data. We used a variance decomposition approach, which, as implemented in the computer package SOLAR (27, 50), makes possible analyses of data from pedigrees of arbitrary size and complexity. Per classical quantitative genetics theory and method, this approach models the phenotypic covariance for each trait of interest within a pedigree as Ω = 2 Φ σ G 2 + I σ E 2 , where Φ is a matrix of kinship coefficients for all relative pairs, σ G 2 is the additive genetic variance, I is an identity matrix, and σ E 2 is the environmental variance. These components of the phenotypic variance are additive, such that σ P 2 = σ G 2 + σ E 2 . Heritability, or the proportion of the phenotypic variance attributable to additive genetic effects, was estimated as h 2 = σ G 2 / σ P 2 (with 1 − h2 = e2 being variance due to nonadditive genetic factors). We modeled a phenotype of an individual as a linear function of the measurements on the trait, the means of these traits in the population, the covariates and their regression coefficients, plus additive genetic values and random environmental deviations. (The mean effects of sex and age, which served as a proxy for dental wear, were included as covariates when/if they were found to be significant.) We assessed the significance of the maximum likelihood estimates for heritability and other parameters by means of likelihood ratio tests. Twice the difference of the maximum likelihoods of a general model (in which all parameters are estimated) and a restricted model (in which the value of a parameter to be tested is held constant at some value, usually 0) are compared. This difference is distributed asymptotically approximately as either a 1/2:1/2 mixture of c2 and a point mass at 0 for tests of parameters like h2 for which a value of 0 in a restricted model is at a boundary of the parameter space, or as a c2 variate for tests of covariates for which 0 is not a boundary value (51). In both cases, degrees of freedom are equal to the difference in the number of estimated parameters in the two models (52). However, in tests of parameters like h2, whose values may be fixed at a boundary of their parameter space in the null model, the appropriate significance level is obtained by halving the P value (51). Extant Phenotypic Data. We collected standard linear size measurements [following Swindler (53)] from 723 extant OWMs (representing seven genera and nine species) and 199 extant apes (representing four genera and five species) (SI Appendix, Tables S1 and S6). We also collected standard linear size measurements for 56 fossil OWMs (seven genera) and 165 fossil apes (15 genera and 26 species) (SI Appendix, Tables S4 and S7). Data included in our phenotypic analyses are published in Dataset S1. Descriptive statistics are presented in SI Appendix, Tables S2, S5, and S8. Phylogenetic Analyses. We used a consensus molecular chronogram from a Bayesian phylogenetic analysis of genetic data downloaded from the web site 10kTrees, version 3 (54). The 10kTrees dataset is based on sequences of 11 mitochondrial and six autosomal genes sampled from the GenBank. Details on tree inference methodology, including fossil calibration and inference of divergence dates, for the consensus chronogram are available through the 10kTrees documentation. One taxon missing from the 10kTrees sample was Presbytis rubicunda, which was added as a sister taxon to Presbytis melalophos in the consensus tree using the R package ape, version 3.0-11 (55). The divergence date for P. rubicunda and P. melalophos has been contentious, in part, because mitochondrial data suggest that populations historically called P. melalophos are paraphyletic with P. rubicunda nested within the historical P. melalophos clade, presumably due to hybridization or incomplete lineage sorting (56). To address this uncertainty, we used two divergence dates for P. rubicunda: 1.31 Ma based on mitochondrial DNA (56) and 2.5 Ma based on nuclear DNA (57, 58). Results were nearly identical between the two analyses, so we chose the date of 1.31 Ma (56) for final analyses. Divergence dates for other taxa as inferred by 10kTrees and referring to the nodes in SI Appendix are provided in SI Appendix, Fig. S3 and Table S11. Phylogenetic Least Squares Regression. To understand better how variation in the MMC and PMM is distributed across primates while taking evolutionary relationships into account, we fit each factor within a phylogenetic mixed model (59) using a Bayesian framework in the R package MCMCglmm, version 2.17 (60, 61). The Bayesian framework was in keeping with the tree estimation methodology used to estimate tree structure used by 10kTrees and enable us to incorporate phylogenetic uncertainty into our models. To control for possible sex and body size effects, we included the mesial width of the left M 2 (a noted proxy for body size) and sex as random effects along with the consensus molecular chronogram from 10kTrees described above. We used noninformative priors corresponding to an inverse-gamma distribution with shape and scale parameters equal to 0.01. For both the MMC and PMM, we allowed the Markov chains a burn-in period of 4,000 iterations, after which we ran 140,000 iterations while sampling every 20th iteration for the posterior distribution. Pagel’s λ was then estimated from the posterior distribution, along with the 95% credibility intervals. The significance of the deviation of Pagel’s λ from 0 and 1 was tested with mean values for all extant taxa using the gls function in the R package nlme, version 3.1-113 (59). Results are presented in SI Appendix, Tables S9 and S10. Morphological Disparity Through Time. Morphological disparity through time was calculated for extant catarrhine taxa measured in this study following the protocol described by Harmon et al. (62) and using the R package geiger, version 2.0.3 (63). For each taxon, the mean value of the MMC and PMM for all extant taxa was used for analysis. Disparity is calculated from average pairwise Euclidian distance between species. Relative disparities for each subclade are represented by the subclade’s (or preceding node’s) disparity divided by the disparity of the entire clade. Values near 0 imply that a particular clade contains little of the overall variation and that variation in the trait is partitioned between subclades rather than within them, whereas values near 1 suggest that a clade contains a large amount of that variation and that clades may overlap in trait space. Results are presented in Fig. 3.

Acknowledgments We thank C. Owen Lovejoy, P. David Polly, Gen Suwa, Tim D. White, and three anonymous reviewers for their feedback on earlier versions of this manuscript. We thank the following people for assistance with data collection and/or project development: Julia Addiss, Stephen Akerson, Sarah Amugongo, Liz Bates, Josh Carlson, Selene Clay, Josh Cohen, Theresa Grieco, Anne Holden, Michaela Huffman, Daniel Lopez, Kurtis Morrish, Jackie Moustakas, Alicia Murua-Gonzalez, Whitney Reiner, Oliver Rizk, Antoine Souron, Risa Takenaka, Kara Timmins, Mallory Watkins, Madsen H. White, Jeffrey Yoshihara, Sunwoo Yu, and Arta Zowghi. We thank the Southwest National Primate Research, supported by NIH National Center for Research Resources Grant P51 RR013986, for access and assistance with the quantitative genetic analysis of baboon dental variation. The following museums provided access to their skeletal collections: American Museum of Natural History, Cleveland Museum of Natural History, Phoebe A. Hearst Museum of Anthropology, University of California Museum of Vertebrate Zoology, and the Smithsonian Institution’s National Museum of Natural History. This work was supported by National Science Foundation Division of Behavioral and Cognitive Sciences Grants 0500179, 0616308, and 0130277 (to L.J.H.). The analyses presented herein relied specifically on the data provided in SI Appendix.

Footnotes Author contributions: L.J.H. designed research; L.J.H. and M.C.M. performed research; C.A.S., T.A.M., M.F.B., and M.C.M. contributed new reagents/analytic tools; L.J.H., C.A.S., T.A.M., M.F.B., and M.C.M. analyzed data; and L.J.H. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Data deposition: Full Old World monkey phenotypic data are archived online at the Dryad Digital Repository (doi:10.5061/dryad.693j8).

See Commentary on page 9142.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1605901113/-/DCSupplemental.