Evolution of Craniodental Correlates of Diet in African Bovidae

J Mammal Evol DOI 10.1007/s10914-016-9329-x ORIGINAL PAPER Evolution of Craniodental Correlates of Diet in African Bovidae Ignacio A. Lazagabaster 1,2 & John Rowan 1,2 & Jason M. Kamilar 2,3,4 & Kaye E. Reed 1,2 # Springer Science+Business Media New York 2016 Abstract Establishing the relationship between craniodental morphology and dietary ecology in extant species permits inferences to be made about the ecology and biology of fossil species and the habitats they inhabited. Previous work linking diet and craniodental morphology has historically relied upon categorical classifications of diet and has not considered the phylogenetic signal (i.e., non-independence) of morphology due to shared evolutionary history. Here we use phylogenetic comparative methods to analyze the relationship between diet and eight craniodental indices for 40 species of African Bovidae using both categorical and continuous (stable carbon isotopes of enamel, δ13C) classifications of diet. In addition, we examine three modes of evolution that best explain interspecific variation in each of these indices, including: Brownian Motion (BM), Early Burst (EB), and Ornstein-Uhlenbeck (OU). Our results show that the hypsodonty index (HI), the length of the masseteric fossa relative to facial depth (MAS-F), and the length of the diastema relative to the total toothrow length (DIAS-TR) are the best predictors of diet among African bovids. These indices are best explained by either a BM or OU mode of evolution. Our findings have important implications * Ignacio A. Lazagabaster [email protected] 1 Institute of Human Origins, Arizona State University, Tempe, AZ 85282, USA 2 School of Human Evolution and Social Change, Arizona State University, Tempe, AZ 85282, USA 3 Department of Anthropology, University of Massachusetts, Amherst, MA 01003, USA 4 Graduate Program in Organismic and Evolutionary Biology, University of Massachusetts, Amherst, MA 01003, USA for understanding the evolution of craniodental traits and reconstructing the diet of fossil mammals, especially bovids. Keywords Stable isotope ecology . Phylogenetic comparative methods . Africa . Bovidae . Hypsodonty Introduction Extant bovids (~ 143 species) comprise one of the most diverse and abundant clades of mammalian herbivores and perform critical roles in the functioning of modern terrestrial ecosystems (Janis 2007; Heller et al. 2013; Hopcraft et al. 2015; Shorrocks and Bates 2015). Bovids are most diverse in Africa (~ 86 living species), where several endemic antilopine tribes (e.g., Aepycerotini, Alcelaphini, Cephalophini, Hippotragini, Reduncini) dominate the herbivore faunas, especially in the savannas of eastern and southern Africa (Kingdon 2013; Groves 2014). Fossils from the middle Miocene onwards document the gradual emergence of bovid-dominated mammal communities in Africa, with major tribal radiations beginning during the latest Miocene and culminating in extraordinary species diversity by the Plio-Pleistocene and leading up to the present day (Gentry 1990, 2010; Vrba 1995; Bibi et al. 2009). Due to their abundance and habitat and dietary specificity, bovids have been extensively used for paleoenvironmental reconstruction (Vrba 1975, 1985, 1988; Bobe and Eck 2001; Bobe and Behrensmeyer 2004; Reed 2008). Most paleoenvironmental reconstructions are based on ecological analogy between fossil and living species (i.e., taxonomic uniformitarianism), although others have relied on linking morphology with habitat and dietary ecology (i.e., ecomorphology). Establishing morphological correlates of habitat and dietary ecology in African bovids permits inferences about the ecology and biology of fossil species, and the J Mammal Evol environments they inhabited (Spencer 1995, 1997; Kappelman et al. 1997). Previous work has investigated the relationship between diet and craniodental morphology for ungulates in general (Jernvall et al. 1996; Pérez-Barbería and Gordon 1999, 2001; Mendoza et al. 2002; Janis 2008; Damuth and Janis 2011) and African bovids in particular (Spencer 1995, 1997; Sponheimer et al. 1999). Most studies have used categorical classifications of diet (e.g., a grazer/browser dichotomy) that may potentially ‘bin’ several distinct feeding ecologies. Relatively recent measures of diet on continuous scales (e.g., carbon stable isotopes (δ13C) of tooth enamel) have become widely available for many extant species (e.g., Cerling et al. 2003, 2015; Sponheimer et al. 2003) and permit investigation of more subtle patterns between diet and morphology. In addition, most research has not considered the effects of phylogenetic signal (i.e., non-independence) of craniodental traits among species (but see Pérez-Barbería and Gordon 1999; Perez-Barberia and Gordon 2001 for exceptions). As defined by Blomberg and Garland, phylogenetic signal is the Btendency for related species to resemble each other more than they resemble species drawn at random from the tree^ (Blomberg and Garland 2002: 905). Phylogenetic signal may possibly influence the relationship between craniodental morphology and diet, as bovid species are likely to exhibit morphologies more similar to their closest relatives (e.g., species within the same tribe) than with more distantly related species to which they have more dietary similarity. Examining craniodental traits in a phylogenetic context allows us to test different modes of evolution by which they may have evolved. There are three major modes of evolutionary change that can be modeled: Brownian Motion (BM), Ornstein-Uhlenbeck (OU), and Early Burst (EB). There are other modes of evolution but these are the three best understood and most commonly used in the literature (Nunn 2011). In a BM model, trait variation increases in proportion to time and is proportional to the sum of the branch lengths from the root to the tips (Freckleton and Harvey 2006; Cooper and Purvis 2010; Harmon et al. 2010; Kamilar and Cooper 2013). Alternatively, in an OU model, a trait evolves towards an optimal phenotype, which can result from natural selection toward an optimum trait value (Lande 1976; Felsenstein 1988; Hansen 1997; Cooper and Purvis 2010; Harmon et al. 2010). As argued by Harmon et al. (2010), an OU model may also imply neutral evolution taking place in a tightly constrained part of morphospace such that trait values are limited to a certain range. Lastly, an EB model implies an early rapid diversification of the trait among lineages, with the rate of change slowing as time progresses (Harmon et al. 2010). This mode of trait evolution is consistent with adaptive radiations (Cooper and Purvis 2010), understood here as Bmore or less simultaneous divergences of numerous lines from more or less the same ancestral adaptive type,^ as originally described by Simpson (1955). Here we use phylogenetic comparative methods to analyze the relationship between craniodental morphology and diet in African bovids using both categorical and continuous classifications of diet. We then analyze the mode of evolution of each morphological correlate of diet. Understanding the evolution of bovid craniodental traits is important, as Cantalapiedra et al. (2014) have shown that the diversification of ruminants was tied to dietary evolution over the Cenozoic. Thus, morphological traits tightly linked to dietary ecology are likely to have been important in facilitating the divergence and subsequent diversification of several ruminant groups, including bovids. Materials and Methods Data Collection We compiled craniodental data for 40 species of African bovids from Mendoza et al. (2002) and generated eight morphological indices that have been previously linked with diets in ungulates (Table 1; Fig. 1). We assigned each species a categorical dietary class (browser, grazer, mixed-feeder) based on literature sources (Kingdon 2013; Cantalapiedra et al. 2014). In addition, we collected stable carbon isotope (δ13C) data for 35 species from various sources (Cerling et al. 2003, 2015; Sponheimer et al. 2003; Codron et al. 2005, 2007; Codron and Brink 2007; Louys and Faith 2015). Because the δ13C values were calculated from different organic sources (e.g., feces, bone, teeth), they were standardized by converting them to dietary δ13C values and averaged for each species, following the procedure of Louys and Faith (2015). The eight morphological indices included: 1) premolarmolar row length (PR-MR), the length of the lower premolar row divided by the lower molar row; 2) muzzle width-palatal width (MW-PW), the width of the muzzle at the premaxillarymaxillary junction divided by the width of the palate between the maxillary second molars; 3) hypsodonty index (HI), the unworn height of the lower third molar divided by its width; 4) diastema-toothrow length (DIA-TR), the length of the lower diastema divided by the length of the entire lower toothrow; 5) incisor index (IW1-IW2), the width of the lower first incisor divided by the lower second incisor; 6) masseteric fossa-facial depth (MAS-F), the depth of the masseteric fossa relative to facial depth under the orbit; 7) the length of the paraoccipital process divided by the length of the base of the posterior part of the skull (PP-BS); and 8) the basicranial angle (CA). Dietary Correlates We used phylogenetic comparative methods to analyze the relationship between morphology and diet. We used a J Mammal Evol Table 1 African bovid species with body mass (BM), diet category (B = browser, M = mixed feeder, G = grazer) and carbon stable isotopes (δ13C). Body mass values come from Mendoza et al. (2002). Dietary categories come from Kingdon (2013) and Cantalapiedra et al. (2014) δ13C Ref. (δ13C) Species (n = 40) Tribe BM Diet Addax nasomaculatus Aepyceros melampus Hippotragini Aepycerontini 111 53 G M – −18.83 – 1,2 Alcelaphus buselaphus Alcelaphini 136 G −12.16 1,2 Ammotragus lervia Antidorcas marsupialis Caprini Antilopini 86 31 M M – −24.20 – 2 Cephalophus dorsalis Cephalophini 20 B −28.73 1,2 Cephalophus silvicultor Cephalophus spadix Connochaetes gnou Cephalophini Cephalophini Alcelaphini 61 57 136 B B G −28.35 −27.50 −15.70 1,2 1 3 Connochaetes taurinus Alcelaphini 216 G −13.01 1,2 Damaliscus pygargus Alcelaphini 69 G −13.60 1,2 Dorcatragus megalotis Antilopini 9 M – – Eudorcas thomsonii Gazella dorcas Hippotragus equinus Antilopini Antilopini Hippotragini 20 23 270 M M G −17.10 – −12.83 1,2 – 1,2 Hippotragus niger Hippotragini 226 G −11.91 1,2 Kobus ellipsiprymnus Kobus leche Litocranius walleri Reduncini Reduncini Antilopini 205 87 42 G G B −13.52 −18.20 −25.73 1,2 2 1,2 Madoqua kirkii Nanger granti Neotragus moschatus Oreotragus oreotragus Antilopini Antilopini Neotragini Oreotragini 5 62 5 14 B B B B −25.92 −24.02 −27.39 −25.10 1 1,2 1 1,2 Oryx gazella Ourebia ourebi Pelea capreolus Philantomba monticola Raphicerus campestris Redunca arundinum Redunca fulvorufula Sylvicapra grimmia Syncerus caffer Tragelaphus angasii Tragelaphus buxtoni Tragelaphus eurycerus Tragelaphus imberbis Tragelaphus oryx Tragelaphus scriptus Tragelaphus spekii Tragelaphus strepsiceros Hippotragini Antilopini Reduncini 170 18 32 G M B −15.30 −16.03 – Cephalophini Antilopini Reduncini Reduncini Cephalophini Bovini Tragelaphini Tragelaphini Tragelaphini Tragelaphini Tragelaphini Tragelaphini Tragelaphini Tragelaphini 6 14 62 31 13 620 91 183 205 77 511 58 74 215 B B G G B G B B B B B B B B −27.68 −26.33 −13.55 −12.56 −26.06 −14.39 −24.50 −26.92 −27.30 −24.61 −24.79 −27.65 −26.35 −25.47 4 1,2 – 1 1,5 1,2 1,2 1,2 1,2 1,5 1 1 1 1,2 1,2 1,2 1,2 1 5 Cerling et al. (2015); 2 Louys and Faith (2015); 3 Codron and Brink (2007); 4 Codron et al. (2005); Codron et al. (2007) phylogeny of all species in our dataset based on Hassanin et al. (2012). All indices were log-transformed to ensure normality. Dietary classes (ANOVA) We performed ANOVAs controlling for phylogenetic non-independence using the function aov.phylo in the package geiger (Pennell et al. 2014; Harmon et al. 2015). ANOVAs were used to predict categorical classifications of diet from craniodental indices for 40 species. δ13C values (PGLS) We used δ13C values of 35 species as a continuous measure of diet and analyzed the relationship J Mammal Evol a PP MAS BS MZ b PW c DIA d Fig. 1 Craniodental measurements used in this study. The length and width of the third molar (hypsodonty index, HI), and the angle of the braincase (CA) are not shown. F = depth of the face under the orbit; MAS = length of the masseteric fossa; PP = length of the paraoccipital process; MZ = muzzle width; PW = palatal width; DIA = length of the diastema; PR = premolar row length; MR = molar row length; IW1 = width of the first lower incisor; IW2 = width of the second lower incisor. The cranium pictured belongs to Alcelaphus buselaphus (USNM 161953 from the National Museum of Natural History, Washington, D.C.) between δ13C values and craniodental indices using regression. First, we used phylogenetic generalized least-squares (PGLS) regressions to predict δ13C values from each of the craniodental indices in a series of bivariate models. Second, we analyzed the relationship between δ13C and all of the craniodental indices using PGLS multiple regression. We used corrected Akaike Information Criterion (AICc) to determine what indices were the best predictors of diet by generating the sum of AICc weights for each craniodental index. AICc was also used to weigh various combinations of indices and we consider all models with a Δ AICc value less than 2 to be ‘equally good’ following Burnham and Anderson (2002). We used the function pgls in the caper package (Orme et al. 2013) for the PGLS analyses. AICc multimodel inference was performed using the dredge function from the MuMIn package (Barton 2014). Following Burnham and Anderson (2002), we averaged the models within the top 95 % of model weight. Phylogenetic PCA We used phylogenetic principal components analysis (PCA) to compress the craniodental indices into a series of components while controlling for phylogeny using the function phyl.pca from the phytools package (Revell 2009, 2012). The eigenvalues of principal component 1 (PC1) against PC2 and PC1 against PC3 were then plotted in bivariate plots to visualize the discrimination of diets by craniodental indices. We plotted PC3 as it explained a similar amount of variance to PC2 (see results). Evolutionary Mode All of the craniodental indices were fitted to three types of evolutionary models: Brownian Motion (BM), Ornstein-Uhlenbeck (OU), and Early Burst (EB). We considered the model with the lowest AICc value the best fit for each index. However, models with a Δ AICc value less than 2 were considered to be ‘equally J Mammal Evol good’ following Burnham and Anderson (2002). To model the evolution of craniodental indices, we used the function fitcontinuous from the geiger package in R (Harmon et al. 2015). Table 3 Summary statistics of phylogenetic generalized least squares (PGLS) bivariate regressions predicting carbon stable isotopes (δ13C) from craniodental indices Phylogenetic bivariate regression Results Dietary Correlates ANOVA (dietary classes) Phylogenetic ANOVAs demonstrated that HI is the strongest predictor of diet among African bovids using categorical dietary classes but that MZ-PW, MAS-F, and PR-MR are also significant predictors (p < 0.05) (Table 2). CA and PP-BS were non-significant (p > 0.05), while DIA-TR and IW1-IW2 had p-values approaching significance (p = 0.05). PGLS (δ13C) The bivariate PGLS regressions indicated that all indices are strong predictors of δ13C isotopes (p < 0.05) with the exception of PP-BS, which showed very low correlation with diet (p = 0.03, r2 = 0.11) (Table 3; Fig. 2). The PGLS multiple regression analyses, however, showed that only HI, DIA-TR, MAS-F, MZ-PW, and IW1-IW2 are significant predictors of δ13C values (overall model r2 = 0.557, p < 0.001). These variables have sum AICc weights over 0.50 (Table 4A). Sum AICc weights for HI and DIA-TR are 1, indicating that they were included in all models and are the most important predictors of diet. CA, PR-MR, and PP-BS have sum AICc weights under 0.15 and are poor predictors. HI, DIA-TR, MAS-F, and IW1-IW2 were included in the top model, but the second model also included CA, while many of the subsequent models had MZ-PW but not IW1-IW2 (Table 4B). Phylogenetic PCA The results of the PCA are shown in Table 5 and Fig. 3. PC1 explains 28.3 % of the variance, while PC2 and PC3 explain 18.9 % and 16.5 % of the variance, respectively. Together the first three components explain 63.7 % of the variance overall. The indices that are loading most strongly on PC1 are MAS-F, PR-MR, and CA, but also Table 2 Summary statistics of phylogenetic one-way analysis of variance (phylogenetic ANOVA) predicting categorical diets (B, M, G) from craniodental indices Phylogenetic ANOVA Variable HI MAS-F MZ-PW CA DIA-TR IW1-IW2 PR-MR PP-BS F-statistic 17.485 11.292 12.597 7.234 9.118 8.418 12.664 3.658 p-value 0.007 0.030 0.019 0.098 0.050 0.050 0.015 0.221 Variable Estimate t-statistic p-value Adj.r2 HI MAS-F 11.378 −21.289 5.520 −5.006 <0.001 <0.001 0.464 0.414 MZ-PW CA DIA-TR IW1-IW2 17.950 −100.321 15.154 −7.005 5.946 −3.792 3.720 −3.332 <0.001 0.001 0.001 0.002 0.503 0.282 0.274 0.229 PR-MR PP-BS −14.726 −20.669 −3.498 −2.249 0.001 0.031 0.248 0.107 HI, IW1-IW2, and DIA-TR. The indices that are loading more strongly on PC2 are PP-BS, MZ-PW, and CA, while those that load more strongly on PC3 are DIA-TR, MZ-PW, and HI (Table 5). Evolutionary Mode The AICc values obtained by fitting the three modes of evolution to each of the craniodental indices are illustrated in Fig. 4. There is a lot of variability in the mode of evolution the craniodental indices followed. Two indices (DIA-TR and PR-MR) follow a Brownian Motion (BM) model, while two others (IW1-IW2, and MZ-PW) follow and Early Burst (EB) model (Figs. 4 and 5). The remaining four indices (HI, MASF, PP-BS, and CA) follow an Ornstein-Uhlenbeck (OU) model (Table 6). It is important to note that some of the indices have AICc values that are less than two units different from the AICc value of another model. In these cases, we consider that the second model could also be valid for explaining the evolution of that index (see footnotes in Table 6). Discussion This is the first study to our knowledge that has explored the relationship between craniodental morphology and diet in African bovids accounting for phylogenetic effects and investigating the mode of trait evolution. In addition, we have used δ13C isotopes as a continuous proxy for diet, together with a more traditional tripartite categorical classification of diet. Integration of the δ13C values is important, as stable carbon isotopes are an increasingly common proxy in paleoecological studies of fossil bovids and other ungulates in Africa (e.g., Cerling et al. 2011, 2015; Bedaso et al. 2013; Bibi et al. 2013). Our results support previous studies that have found a strong relationship between diet and hypsodonty, diastema length, muzzle width, masseteric fossa length, the ratio of the J Mammal Evol MAS-F =0.41 MZ-PW HI =0.50 =0.46 13C 13C 13C =0.28 Browser CA DIA-TR Mixed feeder Grazer =0.27 13C 13C =0.11 PP-BS IW1-IW2 =0.23 =0.25 13C 13C 13C PR-MR Fig. 2 Bivariate plots showing the relationship between log-transformed craniodental indices and carbon stable isotopes (δ13C) in African bovids. Points are colored by diet category: Browser – red circles; Mixed feeder – blue squares; Grazer – green triangles lower first and second incisor widths, basicranial angle, the lower premolar-molar row length ratio, and the length of the paraoccipital process (Janis and Ehrhardt 1988; Janis and Fortelius 1988; Janis and Thomason 1995; Williams and Kay 2001; Perez-Barberia and Gordon 2001; Mendoza et al. 2002; Kaiser et al. 2013). However, the importance of craniodental indices varies across our statistical analyses. Thus, we focus our discussion on HI, MAS-F, and DIASTR, as these three indices were consistently shown to be highly significant predictors of diet. The three indices we found to be the best predictors of diet (HI, MAS-F, DIAS-TR) have clear functional or biomechanical relationships with the physical properties of food items. The HI is higher in grazing species because their teeth are exposed to higher wear rates during the lifespan of an individual than in browsers or mixed-feeders (Janis and Fortelius 1988; Damuth and Janis 2011). Higher wear rates are either due to the presence of hard silica particles in grass or due to the increased intake of grit accompanying the food (Damuth and Janis 2011; Hummel et al. 2011; Kaiser et al. 2013; Lucas J Mammal Evol weights; B) Best models. Only models with ΔAICc <2 are shown. [SE = standard error; df = degrees of freedom; log-lik = logarithmic likelihood] Table 4 Phylogenetic generalized least squares (PGLS) multiple regressions predicting carbon stable isotopes (δ13C) from craniodental indices. A) Summary statistics for each craniodental index, including AICc 4 A) 4B) Variable Estimate SE z-statistic p-value Weight HI 3.890 1.450 2.684 0.007 1 DIA-TR MAS-F 7.158 −7.728 2.511 3.142 2.850 2.460 0.004 0.014 1 0.85 MZ-PW 5.388 3.144 1.713 0.087 0.66 IW1-IW2 CA −2.894 −24.790 1.606 12.010 1.802 2.064 0.071 0.039 0.50 0.15 PR-MR 5.456 3.364 1.622 0.105 0.13 PP-BS −6.778 4.246 1.596 0.110 0.12 Model 1 HI DIA-TR MAS-F Variables IW1-IW2 2 3 HI HI DIA-TR DIA-TR MAS-F IW1-IW2 IW1-IW2 MZ-PW MZ-PW 4 HI DIA-TR MAS-F MZ-PW 5 6 7 HI HI HI DIA-TR DIA-TR DIA-TR MAS-F MAS-F MAS-F MZ-PW MZ-PW 8 9 10 HI HI HI DIA-TR DIA-TR DIA-TR MAS-F MZ-PW MZ-PW MAS-F IW1-IW2 CA PR-MR PR-MR df 5 log-lik −81.46 ΔAICc 0.00 Weight 0.06 6 6 −80.40 −80.47 0.81 0.95 0.04 0.04 6 −80.51 1.04 0.04 PP-BS 6 5 4 −80.56 −82.03 −83.44 1.13 1.15 1.23 0.04 0.04 0.04 PP-BS 7 5 6 −79.05 −82.33 −80.88 1.28 1.75 1.77 0.03 0.03 0.03 CA PP-BS Total model: F = 6.34; SE = 6.831; Adj. r2 = 0.557; p < 0.001 et al. 2014). Likewise, grazers also tend to have wider masseteric fossae for the insertion of larger masseter muscles, which increases the strength and efficiency of the masticatory apparatus when processing relatively abrasive grasses and associated grit (Janis and Thomason 1995; Mendoza et al. 2002). Grazers also have longer diastemata than browsers and mixedfeeders. The function of the diastema, however, is still Table 5 Summary of phylogenetic principal components analysis (phylogenetic PCA) on craniodental indices of African bovids. The loadings for the first three principal components for each index and the variance explained by each principal component are shown Variable PC1 PC2 PC3 PR-MR IW1-IW2 HI MZ-PW DIA-TR PP-BS CA MAS-F Eigenvalue Variance (%) Cumulative Variance (%) −0.420 −0.370 0.372 0.208 0.044 −0.014 −0.418 −0.573 1.505 28.3 28.3 −0.302 0.271 0.256 0.449 0.107 −0.603 0.432 0.084 1.230 18.9 47.2 0.162 −0.231 −0.384 0.559 0.646 0.205 0.029 0.007 1.147 16.5 63.7 debated. According to the ‘fracture scaling hypothesis,’ the length of the diastema varies as a consequence of oversized jaws in grazers (Lucas 2004). Thus, the diastema may be sizedependent and highly correlated to the overall size of the mandible. As grazers are generally larger than browsers, the diastema may be seen as simply a consequence of their larger body mass. However, it should be noted that there are some African bovid species with relatively small body mass, such as Litocranius walleri or Nanger granti, that present very long diastemata (Fig. 5). Greaves (1978) postulated that the presence of the diastema was related to size but that the length of the posterior molar row was constrained to improve the effectiveness of the bite force. The mode of evolution that best explains interspecific variation in craniodental traits is very heterogeneous, including among those that most strongly predicted diet. Some of the indices, such as IW1-IW2 and MZ-PW, follow an EB mode of evolution. This implies that these indices evolved early in the radiation of the clade and then their evolution subsequently slowed down through time. It is possible that these indices evolved rapidly during the adaptive radiation of bovid tribes in Africa and for biomechanical or adaptive reasons they remained virtually unaltered in the subsequent evolutionary history of each tribe. If this is indeed true, fossil species, especially those at the base of the clade, should show similar trait values to those of their modern relatives. Testing this J Mammal Evol J Mammal Evol Phylogenetic principal components analysis (PCA) plot of all craniodental indices. Top, plot of PC1 against PC2; bottom, plot of PC1 against PC3. Taxa are colored by dietary category: Browser – red circles; Mixed feeder – blue squares; Grazer – green triangles ƒFig. 3 hypothesis is out of the scope of this paper but should be the focus of future studies incorporating the fossil record of African Bovidae. Cantalapiedra et al. (2014) suggested that early bovids were most likely mixed feeders and it can therefore be expected that the morphology of their skull and dentition was similar to the morphology of mixed-feeders today. Thus it appears that IW1-IW2 and MZ-PW, both of which follow an EB mode of evolution, evolved very early in tribal radiations and are examples of cases where expanding morphological and ecological disparity was coupled with increasing taxonomic diversity. Other traits are best explained by an OU model of evolution, suggesting that there is an optimal adaptive plateau (a mean value) towards which the trait is evolving. There are many reasons for a trait to evolve towards an optimum and causation is difficult to decipher. The most plausible explanations involve either some type of functional constraint or a biomechanical optimization of that trait. The best predictor of diet according to our analyses is HI and this trait is better explained by an OU mode of evolution (although see footnote in Table 6). In some cases, the values of this index are highly variable among closely related species. For example, Litocranius walleri is a very brachydont browsing species (HI = 1.32), while its sister taxon, the mixed-feeder Antidorcas marsupialis, is comparatively hypsodont (HI = 4.89) (Fig. 5). These differences are higher than what is expected according to branch length, which is equal to time of evolutionary divergence assuming a constant rate of change. We hypothesize that the primary cause of this divergence is diet. This argument may still be sustained if we prefer to consider that the traits were better suited to a BM mode of evolution, because the within-clade variability of this trait is still comparatively high (and thus the phylogenetic signal is low). The craniodental variables MAS-F, PP-BS, and CA also follow an OU model and the same argument explained above Conclusions Establishing craniodental correlates of dietary ecology in extant species permits inferences to be made about the ecology and biology of fossil species and the habitats they inhabited. Previous work linking diet and craniodental morphology has historically relied upon categorical classifications of diet and has not considered the phylogenetic signal (i.e., non-independence) of morphology due to shared evolutionary history. Here we show using phylogenetic comparative methods that the best predictors of diet in African bovids are the hypsodonty index, the length of the diastema relative to that of the toothrow, and the length of the masseteric fossa in relation to the height of the face. These indices are best suited to either OU or BM evolutionary mode. We hypothesize that diet is a key factor in explaining their evolution. As habitats changed during the Cenozoic, bovids adapted to new dietary niches by modifying their skull and their teeth. Some of these evolutionary adaptations resulted in convergences in different clades, suggesting that some traits were selected towards an optimal value (i.e., those following an OU mode of evolution). Other indices studied in this work were not powerful discriminants of diet in African bovids and were better explained by EB. If dietary diversification was responsible in part for the evolution of these traits, then the most important changes happened early during the radiation of the different tribes and then stayed similar because no substantial changes were necessary to optimize intake and processing of food. Overall, this study provides a framework 0.37 Weighted AICc values Fig. 4 Relative AICc weights of the three evolutionary models (BM, OU, EB) calculated for each craniodental index. Lower relative AICc weights indicate better fit can be followed. However, it is reasonable to expect that the less correlation there is between a craniodental index and diet, the greater the likelihood of other forces acting on it. For example, in the case of PP-BS and CA, these traits are either in the posterior part of the skull or occur in relation to the braincase and it is reasonable to expect that there is more biomechanical constraint associated with morphology of the highly integrated vertebral trunk-skull complex. 0.36 0.35 0.34 BM 0.33 OU 0.32 EB 0.31 0.30 HI MAS-FAC CA DIA-TR IW1-IW2 PR-MR Craniodental variables PP-BS MZ-PW J Mammal Evol Mode of evolution OU/BM BM/OU OU/BM EB EB BM OU OU Caprini Hippotragini Alcelaphini Cephalophini Oreotragini Reduncini Bovidae Antilopini Neotragini Aepycerontini Bovini Tragelaphini 0 2 4 6 HI 0 1 2 3 0 1 2 3 0 .3 .6 .9 0 2 4 6 50 150 0 .3 .6 .9 CA DIA-TR MAS-F MZ-PW IW1-IW2 Browser Mixed feeder Grazer 0 .3 .6 .9 PR-MR PP-BS Craniodental indices Fig. 5 Phylogenetic tree of African Bovidae based on Hassanin et al. (2012). The values for each craniodental index (raw values) are shown in barplots colored by diet category: Browser – red; Mixed feeder – blue; Grazer – green to understand the diet of fossil species while controlling for phylogeny using an independent and continuous proxy for diet (stable carbon isotopes of enamel, δ13C). It will also aid Table 6 Best models fitted to craniodental indices of African bovids. The AICc values for each model are shown. The summary includes evolutionary parameters (sigma-square, zeta-zero and alpha). Only the parameters of the best model are shown [loglik = logarithmic likelihood] in improving our knowledge of the evolutionary history of African bovids, the most diverse group of large herbivores in Africa. σ2 α Variable Best model AICc (BM) AICc (OU) AICc (EB) CA DIA-TR HI IW1-IW2 MAS-F OU BM/OU1 OU/BM1 EB/BM1 OU/BM1 −136.959 −13.916 24.167 21.962 −22.657 −142.284 −12.628 23.423 24.104 −23.188 −134.617 −11.573 26.509 21.02 −20.314 74.475 9.120 −8.878 −7.177 14.427 0.007 0.065 0.266 0.524 0.080 5.033 0.679 1.166 1.122 0.692 2.237 NA 0.953 −1.853 0.865 MZ-PW PP-BS PR-MR EB OU BM −52.556 37.985 −26.701 −50.213 34.055 −24.358 −59.95 40.327 −24.617 33.308 −13.694 15.512 0.176 0.604 0.047 0.795 3.356 −0.592 −3.281 2.357 NA log-lik Z0 1 Best models with AICc values separated from other model from less than 2 units. In the case of DIA-TR, the second best model is OU. In the case of HI, IW1-IW2 and MAS-FAC, the second best model is BM σ2 (sigma-square): rate of evolution, also known as Brownian motion rate parameter α (alpha): in cases where an OU model is selected, alpha quantifies the strength of attraction to optima, also known as the rubber band parameter of the OU process. Where the best fitting model is EB, alpha is the exponential rate change parameter, with positive values indicating an acceleration of rates through time and negative values indicating a deceleration Z0 (zeta-zero): ancestral state or root node value J Mammal Evol Acknowledgments We thank the editors and reviewers of Journal of Mammalian Evolution for their most helpful comments. We thank J. Cantalapiedra for helpful discussions and suggestions. I. A. Lazagabaster was supported by Fundacion Obra Social La Caixa Graduate Fellowship. J. Rowan was supported by a National Science Foundation Graduate Research Fellowship. References Barton K (2014) Package BMuMIn^: multi-model inference. R Package Version 1.9.13. Available at: https://cran.r-project.org/web/ packages/MuMIn/index.html Bedaso ZK, Wynn JG, Alemseged Z, Geraads D (2013) Dietary and paleoenvironmental reconstruction using stable isotopes of herbivore tooth enamel from middle Pliocene Dikika, Ethiopia: implication for Australopithecus afarensis habitat and food resources. J Hum Evol 64:21–38 Bibi F, Bukhsianidze M, Gentry AW, Geraads D, Kostopoulos DS, Vrba ES (2009) The fossil record and evolution of Bovidae: state of the field. Palaeontol Electron 12:1–11 Bibi F, Souron A, Bocherens H, Uno K, Boisserie JR (2013) Ecological change in the lower Omo Valley around 2.8 Ma. Biol Lett 9: 20120890 Blomberg SP, Garland T (2002) Tempo and mode in evolution: phylogenetic inertia, adaptation and comparative methods. J Evol Biol 15: 899–910 Bobe R, Behrensmeyer AK (2004) The expansion of grassland ecosystems in Africa in relation to mammalian evolution and the origin of the genus Homo. Palaeogeogr Palaeoclimatol Palaeoecol 207:399– 420. doi:10.1016/j.palaeo.2003.09.033 Bobe R, Eck GG (2001) Responses of African bovids to Pliocene climatic change. Paleobiology 27:1–48. doi:10.1666/0094-8373(2001)027 Burnham KP, Anderson DR (2002) Model Selection and Multimodel Inference: A Practical Information-Theoretic approach. Springer, New York Cantalapiedra JL, FitzJohn RG, Kuhn TS, Fernández MH, DeMiguel D, Azanza B, Morales J, Mooers AØ (2014) Dietary innovations spurred the diversification of ruminants during the Caenozoic. Proc R Soc B Biol Sci 281:20132746. doi:10.1098/rspb.2013.2746 Cerling TE, Andanje SA, Blumenthal SA, Brown FH, Chritz KL, Harris JM, Hart JA, Kirera FM, Kaleme P, Leakey LN (2015) Dietary changes of large herbivores in the Turkana Basin, Kenya from 4 to 1 Ma. Proc Natl Acad Sci USA 112:11467–11472 Cerling TE, Harris JM, Passey BH (2003) Diets of East African Bovidae based on stable isotope analysis. J Mammal 84:456–470 Cerling TE, Levin NE, Passey BH (2011) Stable isotope ecology in the Omo-Turkana Basin. Evol Anthropol 20:228–237. doi:10.1002/ evan.20326 Codron D, Brink JS (2007) Trophic ecology of two savanna grazers, blue wildebeest Connochaetes taurinus and black wildebeest Connochaetes gnou. Eur J Wildl Res 53:90–99 Codron D, Codron J, Lee-Thorp JA, Sponheimer M, De Ruiter D, Sealy J, Grant R, Fourie N (2007) Diets of savanna ungulates from stable carbon isotope composition of faeces. J Zool 273:21–29 Codron D, Codron J, Sponheimer M, Lee-Thorp JA, Robinson T, Grant CC, De Ruiter D (2005) Assessing diet in savanna herbivores using stable carbon isotope ratios of faeces. Koedoe 48:115–124 Cooper N, Purvis A (2010) Body size evolution in mammals: complexity in tempo and mode. Am Nat 175:727–738 Damuth J, Janis CM (2011) On the relationship between hypsodonty and feeding ecology in ungulate mammals, and its utility in palaeoecology. Biol Rev 86:733–758. doi:10.1111/j.1469-185X.2011.00176.x Felsenstein J (1988) Phylogenies and quantitative characters. Annu Rev Ecol Syst 445–471 Freckleton RP, Harvey PH (2006) Detecting non-Brownian trait evolution in adaptive radiations. PLoS Biol 4:e373 Gentry AW (1990) Evolution and dispersal of African Bovidae. In: Bubenik GA, Bubenik AB (eds) Horns, Pronghorns, and Antlers: Evolution, Morphology, Physiology, and Social Significance. Springer-Verlag, New York, pp 195–227 Gentry AW (2010) Bovidae. In: Werdelin L, Sanders WJ (eds) Cenozoic Mammals of Africa. University of California Press, Berkeley, pp 741–796 Greaves WS (1978) The jaw lever system in ungulates: a new model. J Zool 184:271–285 Groves C (2014) Current taxonomy and diversity of crown ruminants above the species level. Zitteliana 5–14 Hansen TF (1997) Stabilizing selection and the comparative analysis of adaptation. Evolution 51:1341–1351 Harmon LJ, Losos JB, Davies JT, Gillespie RG, Gittleman JL, Bryan Jennings W, Kozak KH, McPeek MA, Moreno-Roark F, Near TJ (2010) Early bursts of body size and shape evolution are rare in comparative data. Evolution 64:2385–2396 Harmon L, Weir J, Brock C, Glor R, Challenger W, Hunt G, FitzJohn R, Pennell M, Slater G, Brown J (2015) Package Bgeiger.^ Available at: https://cran.r-project.org/web/packages/geiger/index.html Hassanin A, Delsuc F, Ropiquet A, Hammer C, van Vuuren BJ, Matthee C, Ruiz-Garcia M, Catzeflis F, Areskoug V, Nguyen TT (2012) Pattern and timing of diversification of Cetartiodactyla (Mammalia, Laurasiatheria), as revealed by a comprehensive analysis of mitochondrial genomes. C R Biol 335:32–50 Heller R, Frandsen P, Lorenzen ED, Siegismund HR (2013) Are there really twice as many bovid species as we thought? Syst Biol 62: 490–493 Hopcraft JGC, Sinclair ARE, Holdo RM, Mwangomo E, Mduma S, Thirgood S, Borner M, Fryxell JM, Olff H (2015) Why are wildebeest the most abundant herbivore in the Serengeti ecosystem? In: Sinclair ARE, Metzger KL, Mduma SAR, Fryxell JM (eds) Serengeti IV: Sustaining Biodiversity in a Coupled Human-Natural System. University of Chicago Press, Chicago and London, pp 35– 72 Hummel J, Findeisen E, Südekum KH, Ruf I, Kaiser TM, Bucher M, Clauss M, Codron D (2011) Another one bites the dust: faecal silica levels in large herbivores correlate with high-crowned teeth. Proc R Soc B Biol Sci 278:1742–1747. doi:10.1098/rspb.2010.1939 Janis CM (2007) Artiodactyl paleoecology and evolutionary trends. In: Prothero DR, Foss SE (eds) The Evolution of Artiodactyls. The John Hopkins University Press, Baltimore, pp 292–302 Janis C (2008) An evolutionary history of browsing and grazing ungulates. In: Gordon IJ, Prins HHT (eds) The Ecology of Browsing and Grazing. Springer-Verlag, Berlin, pp 21–45 Janis CM, Ehrhardt D (1988) Correlation of relative muzzle width and relative incisor width with dietary preference in ungulates. Zool J Linn Soc 92:267–284 Janis CM, Fortelius M (1988) On the means whereby mammals achieve increased functional durability of their dentitions, with special reference to limiting factors. Biol Rev 63:197–230. doi:10.1111/j.1469185X.1988.tb00630.x Janis CM, Thomason JJ (1995) Correlations between craniodental morphology and feeding behavior in ungulates: reciprocal illumination between living and fossil taxa. In: Thomason JJ (ed) Functional Morphology in Vertebrate Paleontology. Cambridge University Press, New York, pp 76–98 Jernvall J, Hunter JP, Fortelius M (1996) Molar tooth diversity, disparity, and ecology in Cenozoic ungulate radiations. Science 274(5292): 1489–1492 Kaiser TM, Müller DWH, Fortelius M, Schulz E, Codron D, Clauss M (2013) Hypsodonty and tooth facet development in relation to diet J Mammal Evol and habitat in herbivorous ungulates: implications for understanding tooth wear. Mammal Rev 43:34–46. doi:10.1111/j.1365-2907.2011. 00203.x Kamilar JM, Cooper N (2013) Phylogenetic signal in primate behaviour, ecology and life history. Philos Trans R Soc B Biol Sci 368: 20120341 Kappelman J, Plummer T, Bishop L, Duncan A, Appleton S (1997) Bovids as indicators of Plio-Pleistocene paleoenvironments in East Africa. J Hum Evol 32:229–256 Kingdon J (2013) The Kingdon Field Guide to African Mammals. Bloomsbury Publishing Plc, London Lande R (1976) Natural selection and random genetic drift in phenotypic evolution. Evolution 30(2): 314–334 Louys J, Faith JT (2015) Phylogenetic topology mapped onto dietary ecospace reveals multiple pathways in the evolution of the herbivorous niche in African Bovidae. J Zool Syst Evol Res 53:140–154 Lucas PW (2004) Dental Functional Morphology: How Teeth Work. Cambridge University Press, Cambridge Lucas PW, Casteren A van, Al-Fadhalah K, Almusallam AS, Henry AG, Michael S, Watzke J, Reed DA, Diekwisch TGH, Strait DS, Atkins AG (2014) The role of dust, grit and phytoliths in tooth wear. Ann Zool Fennici 51:143–152 Mendoza M, Janis CM, Palmqvist P (2002) Characterizing complex craniodental patterns related to feeding behaviour in ungulates: a multivariate approach. J Zool 258:223–246 Nunn CL (2011) The Comparative Approach in Evolutionary Anthropology and Biology. University of Chicago Press, Chicago Orme D, Freckleton R, Thomas G, Petzoldt T, Fritz S, Isaac N, Pearse W (2013) CAPER: Comparative Analyses of Phylogenetics and Evolution in R. Version 0.5.2. Available at: https://cran.r-project. org/web/packages/caper/index.html Pennell MW, Eastman JM, Slater GJ, Brown JW, Uyeda JC, FitzJohn RG, Alfaro ME, Harmon LJ (2014) Geiger v2. 0: an expanded suite of methods for fitting macroevolutionary models to phylogenetic trees. Bioinformatics 30(15):2216–2218 Pérez-Barbería FJ, Gordon IJ (1999) The functional relationship between feeding type and jaw and cranial morphology in ungulates. Oecologia 118:157–165 Perez-Barberia FJ, Gordon IJ (2001) Relationships between oral morphology and feeding style in the Ungulata: a phylogenetically controlled evaluation. Proc R Soc B Biol Sci 268:1023–1032. doi:10. 1098/rspb.2001.1619 Reed KE (2008) Paleoecological patterns at the Hadar hominin site, Afar regional state, Ethiopia. J Hum Evol 54:743–768 Revell LJ (2009) Size-correction and principal components for interspecific comparative studies. Evolution 63:3258–3268 Revell LJ (2012) Phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol Evol 3:217–223 Shorrocks B, Bates W (2015) The Biology of African Savannahs. Oxford University Press, New York Simpson GG (1955) Major Features of Evolution. Columbia University Press, New York Spencer LM (1995) Morphological correlates of dietary resource partitioning in the African Bovidae. J Mammal 76:448–471. doi: 10.2307/1382355 Spencer LM (1997) Dietary adaptations of Plio-Pleistocene Bovidae: implications for hominid habitat use. J Hum Evol 32:201–228. doi:10.1006/jhev.1996.0102 Sponheimer M, Lee-Thorp JA, DeRuiter DJ, Smith JM, van der Merwe NJ, Reed KE, Grant CC, Ayliffe LK, Robinson TF, Heidelberger C (2003) Diets of southern African Bovidae: stable isotope evidence. J Mammal 84:471–479 Sponheimer M, Reed KE, Lee-Thorp JA (1999) Combining isotopic and ecomorphological data to refine bovid paleodietary reconstruction: a case study from the Makapansgat Limeworks hominin locality. J Hum Evol 36:705–718. doi:10.1006/jhev.1999.0300 Vrba ES (1975) Some evidence of chronology and palaeoecology of Sterkfontein, Swartkrans and Kromdraai from the fossil Bovidae. Nature 254:301–304. doi:10.1038/254301a0 Vrba ES (1985) Paleoecology of early Hominidae, with special reference to Sterkfontein, Swartkrans and Kromdraai. In: Coppens Y (ed) L’Environnement des Hominidés aux Plio-Pléistocène. Masson, Paris, pp 345–369 Vrba ES (1988) Late Pliocene climatic events and hominid evolution. In: Grine FE (ed) Evolutionary History of the BRobust^ Australopithecines. Aldine de Gruyter, New York, pp 405–426 Vrba ES (1995) The fossil record of African antelopes (Mammalia, Bovidae) in relation to human evolution and paleoclimate. In: Vrba ES, Denton G, Burckle T, Partridge T (eds) Paleoclimate and Evolution with Emphasis on Human Origins. Yale University Press, New Haven, pp 385–424 Williams SH, Kay RF (2001) A comparative test of adaptive explanations for hypsodonty in ungulates and rodents. J Mamm Evol 8:207–229