A multivariate ecogeographic analysis of macaque craniodental variation

Abstract Objectives To infer the ecogeographic conditions that underlie the evolutionary diversification of macaques, we investigated the within‐ and between‐species relationships of craniodental dimensions, geography, and environment in extant macaque species. We studied evolutionary processes by contrasting macroevolutionary patterns, phylogeny, and within‐species associations. Materials and Methods Sixty‐three linear measurements of the permanent dentition and skull along with data about climate, ecology (environment), and spatial geography were collected for 711 specimens of 12 macaque species and analyzed by a multivariate approach. Phylogenetic two‐block partial least squares was used to identify patterns of covariance between craniodental and environmental variation. Phylogenetic reduced rank regression was employed to analyze spatial clines in morphological variation. Results Between‐species associations consisted of two distinct multivariate patterns. The first represents overall craniodental size and is negatively associated with temperature and habitat, but positively with latitude. The second pattern shows an antero‐posterior tooth size contrast related to diet, rainfall, and habitat productivity. After controlling for phylogeny, however, the latter dimension was diminished. Within‐species analyses neither revealed significant association between morphology, environment, and geography, nor evidence of isolation by distance. Discussion We found evidence for environmental adaptation in macaque body and craniodental size, primarily driven by selection for thermoregulation. This pattern cannot be explained by the within‐species pattern, indicating an evolved genetic basis for the between‐species relationship. The dietary signal in relative tooth size, by contrast, can largely be explained by phylogeny. This cautions against adaptive interpretations of phenotype–environment associations when phylogeny is not explicitly modelled.

temperatures on animal body size at higher latitudes (Mayr, 1956;Meiri and Dayan, 2003), but in some taxa rainfall may better explain the relationship between size and latitude (Ashton et al., 2000;Millien et al., 2006), including some nonhuman primates (Cardini et al., 2007;Frost et al., 2003). Such contrasting ecological correlates of Bergmann's rule indicate that different selective forces may give rise to the same pattern (Meiri and Dayan, 2003).
In a taxon-wide study of Bergmann's rule in primates, a positive relationship between latitude and body mass was found among non-Malagasy primates (Harcourt and Schreier, 2009). However, at a lower taxonomic level, and after controlling for phylogeny, the pattern only persisted in macaque species living on the Asian continental shelf (Harcourt and Schreier, 2009). Furthermore, size gradients that correlate not with latitude but with longitude have been retrieved for cranial size within several African cercopithecid primates (vervet monkeys : Cardini et al., 2007;red colobus monkeys: Cardini and Elton, 2009; greater spotnosed and blue monkey: Cardini et al., 2010;and baboons: Dunn et al., 2013). By contrast, cranial shape varies more strongly along a latitudinal than a longitudinal gradient between several Neotropical species of howler and capuchin monkeys Meloro et al., 2014), and between some (but not all) macaques (Ito et al., 2014). Bergmann's rule is, however, a special case of the more general principle of ecomorphology; Bergmann's rule concerns temperature and body size, but evolutionary ecomorphology investigates multiple interactions between environmental parameters and organismal shape and size.
To date, only a limited number of primate studies have included multiple climate and ecological variables (e.g., Cardini et al., 2007;Harvati and Weaver, 2006;Kamilar et al., 2012;Meloro et al., 2014;Viguier, 2004). From these studies, a mixed pattern of the environmental correlates of morphological size and shape variation in primates emerges. Rainfall and other humidity measures, as indicators of habitat productivity, are relevant in explaining cranial variation in vervet monkeys (Cardini et al., 2007), some Malagasy sifakas (Lehman et al., 2005), and lemurs (Viguier, 2004). In New World capuchin monkeys, however, both rainfall and temperature are important climatic predictors of skull shape (C aceres et al., 2014). A recent environmental analysis of Malagasy strepsirrhine body mass revealed that diet and climate were weak predictors of body size, but that there was a strong phylogenetic effect (Kamilar et al., 2012). In modern humans, signals of population history in cranial variation have been found to be stronger than, or even drive, climatic signatures, highlighting the role of population structure and genetic drift (Betti et al., 2010;Harvati and Weaver, 2006;Roseman and Auerbach, 2015). It is becoming increasingly apparent that primate evolution, within and between species, has been characterized by a complex interplay of different selective forces and neutral processes.
Here, we carry out, to our knowledge, the first detailed multivariate analysis of craniodental dimensions and their relation to geographic distribution, climate and species' ecology in the radiation of macaques (Cercopithecidae: Macaca) in a phylogenetic framework. Macaques are an interesting taxon because they diversified widely and rapidly during times of considerable environmental change in the Pliocene and Pleistocene (Abegg and Thierry, 2002;Brandon-Jones, 1996;Delson, 1980), and because they continue to occupy a range of different habitats across southern, central, eastern, and insular Asia and North Africa today (Fooden, 1982).
The main focus of our study is to elucidate the role of the environment in phenotypic and taxonomic diversification at the macroevolutionary level by studying between-species variation within a given phylogeny. Our choice of phenotype is the dentition and associated cranial structures, as teeth are minimally plastic and are therefore likely to carry stronger evolutionary signals than other phenotypes. We compare the interspecific patterns to the intraspecific patterns in order to infer evolutionary processes. Macaques are well known for occupying a great diversity of environments in terms of climate, geographical distribution, and resource exploitation. Therefore, we also investigate here how multiple and diverse aspects of the environment may have impacted diversification of macaques, by specifically studying how they interact in their association to macaque craniodental variation. We expect the effect of climate on macaque thermoregulation and the influence of resource exploitation on the dentition to underlie the functional links between macaque craniodental morphology and the environment. Specifically, we test if the Bergmann effect is reflected in dental patterns, as expected on the basis of a Bergmannian trend in macaque body mass (Harcourt and Schreier, 2009), and explore its ecological and environmental correlates. Furthermore, we expect that dental patterns, such as relative incisor and molar size, covary with food type and thus in turn with climate and geography.
Instead of studying associations between single selected environmental and phenotypic variables, we employ a fully multivariate approach. We search for multivariate patters that jointly underlie the association between all ecological, geographic, and morphological variables. This exploratory approach requires careful interpretation but allows us to identify the actual complexity (number of independent factors) within this association, without prior specification of the number of variables considered. This approach also allows us to investigate the role of a wider range of potentially relevant variables not (often) explored elsewhere. While partial least squares analysis-our method of choice to relate ecological to morphometric variables-has been used in the literature already, the application of reduced rank regression to identify multivariate morphological clines is (to our knowledge) novel.
The effect of species' phylogenetic relatedness on observed relationships is often either not modelled or simply removed as part of the analysis. Here, we examine in detail how the patterns and magnitude of covariation between phenotype and ecogeography in macaques are influenced by phylogeny. In other words, we estimate the extent to which observed phenotype-environment associations in Macaca are consistent with shared ancestry in order to gauge adaptive interpretations such as evolutionary convergence. To this end, we carry out the between-species analyses both with and without phylogenetic correction.

| Morphometric and contextual data
Twelve macaque species were used in this study (Table 1), selected to capture the ecogeographic diversity across the genus in combination GRUNSTRA ET AL.

| 387
with practical considerations regarding their availability in museum collections. A total of 711 specimens (Table 1)  (Vienna). Preference was given to wild-caught and wild-shot specimens, although occasionally captive specimens or those without provenience data were included to obtain acceptable sample sizes.
A total of 46 linear measurements of tooth size were taken on the permanent maxillary and mandibular teeth. Tooth lengths and breadths were measured for all teeth, complemented by tooth height for the anterior dentition, following a standard approach (e.g., Swindler, 2002).
Teeth on the right side of the jaw were measured where possible; broken or missing teeth were substituted by the left antimere. The key to tooth variable names is described in Table 2, and tooth measurement protocols can be found in Tables S1-S3 in the Supporting Information.
In this article, the incisors and the canines are referred to as the anterior dentition, and the premolars and molars as the posterior (or postcanine) dentition. An additional 17 cranial, maxillary, and mandibular measurements were taken to represent the structure that houses the dentition and to provide a "morphological context" for the dental variables (Table 3). These measurements were only recorded on adult specimens showing full eruption of their third molars to minimize ontogenetic variation. Because the skeletal and dental measurements showed the same results separately as they did combined, we only present the results for the combined phenotypic dataset.
All measurements were taken by the same person (NDSG) with digital dental callipers (Mitutoyo, 573 series, Kawasaki, Japan) with an accuracy of 0.01 mm. Intraobserver mean measurement error (derived from a subsample of 50 specimens, with two replicates for each measurement) was 0.18 mm (an average of 2.2% error of the mean) and the intraclass correlation coefficient was 0.99. All morphometric data used in this study are available at doi:10.5281/ zenodo.182699.

| Between-species analysis
For the interspecific analysis, species means were computed from the phenotypic data by taking the arithmetic mean of the female and male means of each species. We pool male and female phenotypes in this analysis as the environmental and geographic variables are the same for both sexes, and strong sex-specific interaction effects between ecogeography and craniodental variation seem unlikely.
Contextual data of macaque ecogeography were collected from published sources. Henceforth, we use "ecogeography" to refer to geographic, environmental, and ecological parameters relating to macaque spatial distribution, climate, as well as habitat and dietary ecology. We use variables that represent average environmental conditions (e.g., mean temperature, minimum rainfall, median altitude, mean habitat productivity) and variables that reflect environmental heterogeneity itself (e.g., range size variables, climatic seasonality indices, and habitat and dietary breadth). All ecogeographic parameters are described below. were obtained from the PanTHERIA database (Jones et al., 2009). As several macaque species have successfully dispersed to insular  (Abegg & Thierry, 2002), as well as the number of specimens measured for this study (includes subadult individuals). Classification as per Groves (2005 Southeast and East Asia, species were assigned to one of the following categories: (1) island(s) only, (2) mixed, (3) continental mainland only.
Data can be found in Supporting Information Table S4.
Climate variables were selected from among the bioclimatic variables in the WorldClim database (Hijmans et al., 2005) and data were retrieved at a resolution of 2.5 arc-minutes. The variables used in this study are common climatic measures, namely mean, minimum and maximum temperature, annual, minimum and maximum precipitation, and the degree of seasonality in temperature and precipitation (defined in Supporting Information Table S5). Climate data were aggregated for each species by first retaining unique localities only to avoid pseudoreplication, and second, by averaging the climate data across these unique localities to obtain species means (Supporting Information Table   S6). Species' altitudinal median and range were derived from field observations reported in the literature (data and sources are in Supporting Information Table S7).
Macaques can broadly be divided into two ecological groups: species that predominantly occur in broadleaf evergreen (BE) forests, and those that exist in a wide range of forest and nonforest (non-BE) habitats (Fooden, 1982). In addition, habitat breadth, dietary breadth, the degree of frugivory and folivory, and the range in the proportion of fruits in the diet were collated from the literature. Finally, mean male and female adult body masses were retrieved to have a measure of overall body size. Data (including sources) on habitat, diet, and body mass are presented in Supporting Information Tables S8, S9, and S10, respectively.

| Within-species analysis
For the intraspecific analysis we used those species for which adequate spatial and climatic variation exist in our sample. These are Macaca nemestrina (N 5 43 from 28 unique localities), M. fascicularis (N 5 70 from 45 unique localities) and M. mulatta (N 5 44 from 33 unique localities). Only morphological data pertaining to adult, wild specimens with known provenience were used. Figure   1 depicts the geographical distribution of the sample for the three species. An interactive version of these maps showing the topography of the land surface and the sea bed can be accessed at https://nicolegrunstra.github.io/GeoMaps_3_species/. Because of damage to specimens there were missing data, which we substituted using Expectation-Maximization (EM) imputation (Dempster et al., 1977;Gunz et al., 2009). EM imputation uses an iterative algorithm that calculates maximum-likelihood estimates for missing values based on the covariance structure of the observed data.
Because of the paucity of data on ecological variation between populations within species, only elevation and climate were used in the  Tooth class I 5 incisor, C 5 canine, P 5 premolar, M 5 molar.

Dimension
MD 5 mesiodistal length of incisors and canines, LL 5 labiolingual width of incisors, H 5 height of incisors and canines, BL 5 buccolingual width of canines, L 5 mesiodistal length of premolars and molars (except P 3 ), OL and TL 5 occlusal and total length (includes the sectorial crest) of the lower third premolar (P 3 ), W 5 buccolingual width of premolars, AW and PW 5 anterior and posterior buccolingual width of molars.
| 389 intraspecific analysis. Both types of data were gathered from the WorldClim database (Hijmans et al., 2005) and analyzed on specimen level. The same eight climatic variables as in the between-species analysis were used. Latitude and longitude were derived from the specimens' locality data.

| Statistical analysis
Because the association between sexual dimorphism and ecogeography is outside of the scope of this article, males and females were pooled in the between-species analysis (the species samples were approximately In both the between-and the within-species analyses, we employed two-block partial least squares analysis (2B-PLS; Bookstein et al., 2003;Mitteroecker and Bookstein, 2007;Rohlf and Corti, 2000) to investigate the multivariate relationship between the environment (block 1) and the morphometric data (block 2). The 2B-PLS finds axes of successively maximum covariance between blocks. The first dimen- Even though 2B-PLS has also been used to study the multivariate association of morphological and geographic variables (Frost et al., 2003), we use a different method here. In ecology, such an association is typically construed as a spatial cline or gradient, represented by the slope of the surface that results from mapping a particular biological variable on a geographic map. Locally, this slope can be estimated by regressing the variable on both latitude and longitude. The two resulting partial regression coefficients (one for latitude, one for longitude) determine the spatial direction with maximum regression slope, i.e., with the steepest local gradient on the surface. In the current multivariate context, this translates into finding a linear combination of morphological variables that has maximum slope when regressed on a linear combination of geographic coordinates. This is achieved by a singular value decomposition of the p 3 2 matrix of partial regression coefficients of all p morphological variables on latitude and longitude (for a mathematical proof see Mitteroecker et al., 2016). The singular values equal the maximal slopes, and the singular vectors contain the morphological and geographic loadings that determine the corresponding LVs. This approach is similar to reduced rank regression (Izenman, 1975), hence we use this name to refer to our multivariate strategy here. Similarly to 2B-PLS, reduced rank regression yields two pairs of latent variables in this application, but they maximize the regression Because of their phylogenetic history, species' data are not statistically independent (Felsenstein, 1985). We therefore performed a phy- This second, primarily ecological factor represents the association of high rainfall and habitat productivity, low rainfall seasonality, a high percentage of fruits in the diet, low variation in the amount of fruits, and a low percentage of leaves in the diet with an antero-posterior dental contrast GRUNSTRA ET AL.

| R E SU LTS
3.1 | Between-species analysis

| Environment
Without phylogenetic correction, 2B-PLS yielded two latent variables that accounted for 63% (LV 1) and 34% (LV 2) of the squared covariance between the blocks of variables (see scree plot in Supporting Information Figure S2a). The correlation between blocks was strong along both dimensions, with LV 1: r 5 0.81, and LV 2: r 5 0.81 (we do not report p-values for the interspecific analyses as their meaning is limited in this small sample of selected species). After phylogenetic correction, however, 2B-PLS extracted only one latent variable (LV 1) with considerable covariance; in fact it accounted for 94% of the total squared covariance between blocks (r 5 0.67). The contribution of LV 2 decreased to 4% (see scree plot in Supporting Information Figure S2b), despite a strong correlation between blocks (r 5 0.85).
The PLS loadings of the environmental and morphological variables on LV 1 showed a very similar pattern before and after phylogenetic adjustment (Figures 2 and 3). Body size and temperature seasonality had high positive loadings, whereas temperature (mean, maximum, and minimum), geographic range size, habitat breadth, and ecological group had high negative loadings on LV 1. All morphological variables loaded The loading patterns of LV 2 also showed a highly similar pattern before and after accounting for phylogeny (Figures 4 and 5). Annual and minimum precipitation, minimum temperature, AET, and percentage of fruit in the diet had relatively high positive loadings for LV 2, whereas range variables, seasonality, measures of dietary variability, and percentage of leaves in the diet all loaded negatively on LV 2. The associated craniodental pattern showed a tooth size contrast (Figures   4b and 5b). Measurements pertaining to the anterior dentition loaded positively on LV 2, i.e., in the same direction as precipitation levels and percentage of fruits. Conversely, measurements of the posterior dentition loaded negatively on LV 2, in the same direction as precipitation seasonality and percentage of leaves. A larger anterior dentition is thus associated with fruit-eating and high degree of rainfall, whereas a larger posterior dentition is associated with more leaf-eating and drier, more seasonal environments. As mentioned, however, the association between blocks along this dimension was greatly diminished once phylogeny was taken into account.

| Geography
The results of the reduced rank regression after phylogenetic correction are presented in Figures 6 and 7. Without phylogenetic adjustment, the results were once again highly similar and they are therefore not presented here. North African M. sylvanus was omitted due to its outlying geographic location. The spatial vectors representing LV 1 and LV 2 strongly corresponded to latitude and longitude, respectively (Figure 6). Prior to phylogenetic correction, LV 1 accounted for 75% of the association (total squared regression slopes) between spatial geography and morphology (r 5 0.91) and LV 2 for the remaining 25% (r 5 0.41) (Supporting Information Figure S4). After phylogenetic correction, however, LV 1 accounted for nearly 100% of the association, whereas LV 2 has now become negligible (see Supporting Information Figure   S4). The effect of phylogeny is further shown by the drop in strength of the correlation coefficient along LV 1 (from 0.91 to 0.39).
Before as well as after phylogenetic adjustment, all craniodental measurements (with a few exceptions) were positively correlated with LV 1 (Figure 7a). Thus, macaque teeth and skulls tend to get larger along a roughly south-to-north gradient. LV 2 was positively associated with measurements pertaining to the anterior teeth and muzzle, and negatively with posterior tooth measurements and calvarium length ( Figure 7b). LV 2 thus discriminates between a relatively larger posterior dentition in the west and a relatively larger anterior dentition in the east, although the effect size of this association is very weak when phylogeny is taken into account.

| Environment
The 2B-PLS returned no significant linear combinations between climate, altitude, and morphology for any of the three species (LV 1: M.

| Geography
Reduced rank regressions also did not find any significant associations between latitude, longitude and morphology for any of the three spe- Results of the reduced rank regression after phylogenetic correction. The loading vectors prior to correction are highly similar and these results are therefore not displayed. M. sylvanus was omitted due to its outlying geographical position, and thus N 5 11. Latent variable (LV) 1, the direction with the steepest morphological cline, corresponds to a south(west)-to-north(east) gradient, and LV 2 to a (north)west-to-(south)east gradient GRUNSTRA ET AL.

| D I SCUSSION
We detected signals of geography, climate, and ecology in the interspecific variation of macaque craniodental morphology, although these patterns are variably mediated by phylogeny. Our between-species analyses demonstrated the presence of only two environmental and spatial gradients in the macaque craniodental phenotype, despite the diversity of variables in the ecogeographic and phenotypic datasets.
The first factor is dominated by overall craniodental size and varies (weakly) along a latitudinal cline, with a tendency for macaques to be smaller near the equator (e.g., M. sinica, M. fascicularis, and M. radiata) and larger at higher latitudes (e.g., M. assamensis and M. fuscata).
Concomitant with this latitudinal cline is the positive relationship of absolute craniodental size with male and female body mass, colder temperatures, and increased temperature seasonality. Taken together, these results are in agreement with a classic Bergmann effect (Millien et al., 2006), and match the positive relationship that has been found between macaque body mass and latitude (Ito et al., 2014;Harcourt and Schreier, 2009). The observed pattern along the first axis was minimally affected by phylogenetic correction: the PLS correlation was only slightly reduced from 0.81 to 0.67, and the effect size along this dimension-the percentage explained of the total squared covariance FIG URE 7 Results from the reduced rank regression with phylogenetic correction. This figure shows the loadings of the morphological variables onto latent variable (LV) 1 (a) and LV 2 (b). LV 2 is diminished in effect size after phylogenetic correction (see the scree plot in Supporting Information Figure S4), and therefore this pattern (b) warrants only limited interpretation. LV 1, on the other hand, describes the association between latitude and overall craniodental size (all morphological loadings are positive, with the exception of canine height) between blocks-remained large (it changed from 63% to 94%).
Furthermore, the loading on LV 1 of temperature, one of the most important variables to contribute to LV 1, even increased. The pattern in macaque craniodental size thus can barely be explained by phylogenetic effects, suggesting that convergent evolution and selection played an important role. Therefore, we interpret the variation in craniodental size-along with overall body size-as an adaptive response to variation in temperature along a latitudinal gradient, and suggest here that species differentiation in Macaca was associated with adaptive diversification in body size.
The second factor discriminates between species with a relatively larger anterior dentition and a more prominent muzzle, and species with a relatively larger posterior dentition and longer calvaria. (We point out that the dental contrast highlighted by LV 2 represents relative craniodental size, because differences due to overall size are captured primarily by LV 1.) The lower third premolar (P 3 ) is part of the CP 3 honing complex in Old World monkeys (Swindler, 2002), and indeed loaded in the same direction as the canines rather than the posterior dentition. The negative association between relative anterior and relative posterior tooth size may reflect the underlying architecture of genetic independence found between incisors and postcanine teeth in baboons and mice (Hlusko et al., 2011). Regardless of whether we account for phylogeny or not, a larger anterior dentition is associated with tropical climates and increased habitat productivity, less variable habitats, small elevational, longitudinal and geographic ranges, and a high subsistence on fruits (e.g., in all representatives of the silenus clade, namely M. nemestrina, M. silenus, and the two Sulawesi macaques). A larger posterior dentition, by contrast, is observed in taxa occupying more temperate regions. These include macaques that experience increased seasonality, occupy a larger variety of habitats and altitudes across larger geographic ranges, and which subsist on proportionally more leaves and highly variable amounts of fruit (e.g., M. mulatta, M. sylvanus, and M. fuscata). This association between relative tooth size and ecogeography may be mediated by the functional link between diet (as influenced by climate and habitat, and hence indirectly also by latitude) and teeth. The association of range size variables with LV 2 reflects the tendency for less frugivorous (and more folivorous) macaque species to also occupy a wider range of habitats, at varying altitudes, and across larger geographical areas.
The environment-craniodental contrast (LV 2) coincides with a longitudinal gradient as long as phylogeny is not accounted for. In fact, in contrast to the first factor, the second pattern can be explained almost entirely in phylogenetic terms. The environment and diet-related variance in craniodental morphology was greatly reduced following phylogenetic correction, resulting in a negligible effect size of LV 2. This is similar to the reduction in effect sizes of the relationships between climate, diet, and macaque craniofacial shape obtained by Ito et al. (2014) after controlling for phylogeny. Although we find no evidence of environmental adaptation in LV 2 once we correct for phylogeny, the pattern itself is similar to what has been found in capuchin monkeys; namely that relative tooth size (of primarily the postcanine dentition) is bigger in species living in relatively cooler, drier, and more seasonal climates (C aceres et al., 2014). C aceres et al. (2014) suggested that species with larger teeth for their body size are able to process a broader range of food items. However, the authors did not employ phylogenetic comparative methods and therefore it is unknown whether the relationship between tooth size and environment in capuchins mostly reflects processes of adaptation or shared ancestry.
Habitat productivity and rainfall patterns were not associated with variation in macaque body size and craniodental size, in contrast to what has been found for baboons (Jolly, 2012), vervets (Cardini et al., 2007), and sifakas (Lehman et al., 2005). This discrepancy between macaques and their cercopithecine relatives in Africa may result from different degrees of variation in environmental factors in different geographic regions. For example, in the tropics and subtropics, patterns of rainfall are more variable than temperature and, hence, may have stronger effects on morphology (DeMenocal and Bloemendal, 1995). Also, many African monkeys may vary morphologically more with longitude than with latitude, because they have wider longitudinal distributions and are therefore subject to environmental variation mainly in that direction.
In addition to the species-level analyses, we also investigated the presence of environmental and spatial gradients as well as isolation by distance (IBD) patterns within species to be able to infer what processes have been important in structuring intraspecific variation in macaques and whether these processes can also explain the variation between species. Intraspecific phenotypic variation that is correlated with environmental or geographic variation primarily reflects phenotypic plasticity, i.e., environmentally induced variation, or, if gene flow is low, genetic differences due to adaptation to local environments. An IBD pattern, by contrast, would result from strong genetic drift in the presence of reduced gene flow. IBD patterns can therefore reflect population history, an intraspecific equivalent to phylogenetic signal (Roseman and Auerbach, 2015). Intraspecific variation in body or skull length has previously been reported to correlate with latitude in M. nemestrina (Albrecht, 1980), M. mulatta (Fooden, 2000), and M. fascicularis (Fooden and Albrecht, 1993;Schillaci et al., 2009). However, we found no relationships between craniodental variation, climate, and spatial geography within the three species in our sample, indicating low phenotypic plasticity in both absolute and relative craniodental size. These results also suggest that local environmental adaptation in the craniodental phenotype is either weak, perhaps due to relatively homogenous environments, or that gene flow is strong, e.g., due to intensive migration.
The absence of detectable intraspecific plasticity supports our claim that the species differences along LV 1 are due to an evolved genetic basis.
Furthermore, we detected no evidence for IBD in craniodental size of M. nemestrina, M. fascicularis, or M. mulatta. The lack of an IBD pattern and spatial clines in these species is in contrast to the IBD found in recent modern humans (Betti et al., 2010)  | 397 unhindered by sea barriers and aided by their ability to move across a variety of habitats owing to their ecological flexibility, may exhibit strong male-mediated gene flow between populations in the sampled region (Figure 1c; Fooden, 2000;Melnick, 1988;Melnick and Hoelzer, 1992;Tosi et al., 2002Tosi et al., , 2003. Among the two ecogeographic gradients, the second (LV 2) is of particular interest. The depicted interspecific association between craniodental variation and diet is in agreement with early comparative work that has linked large incisors (relative to postcanine teeth) to frugivory, and smaller incisors to folivory (Hylander, 1975;Robinson, 1954). Likewise, there is a well-known and pervasive phenomenon among anthropoid primates that postcanine tooth size is often larger in folivores than in closely-related frugivores relative to body or facial size (e.g., Vinyard and Hanna, 2005;Scott, 2011). More recently, this dietmolar pattern has also been found in strepsirhines when adjusted for facial size (Scott, 2012). These two diet-related patterns of relative tooth size are often explained as an adaptive response to masticatory challenges posed by the external properties of food items; large incisors are useful for the ingestion of large, husky, and fleshy fruits, whereas a large postcanine occlusal surface benefits the consumption of small and hard food items (e.g., nuts) or tough, fibrous foods (e.g., mature leaves) (Lucas, 2004;Ungar, 2011). Even though we recovered this classic association between dental dimensions and diet in our analysis, the interspecific differences along this pattern were strongly aligned with phylogenetic relatedness: when phylogeny is statistically accounted for, the pattern of association remains intact, but its magnitude diminishes. Hence, whereas the convergence of LV1 (distantly related taxa are phenotypically similar) provides good evidence for adaptive evolution, the interspecific association for LV 2 can be explained equally well by adaptation and by common ancestry. Phylogenetic patterns are not inconsistent with adaptation per se because adaptive evolution can drive both phenotypic and phyletic divergence, but they can also arise from neutral evolution. In the latter case, phenotypic differentiation occurs mainly by genetic drift and closely related taxa inherit their phenotypic similarity from their common ancestor.
Based on the present data, an adaptive interpretation of the relationship between the relative size of the anterior and posterior dentition and diet is not sufficiently supported.
While significant relationships between diet and dental size have been recovered after phylogenetic correction on higher taxonomic levels (e.g., Scott, 2011Scott, , 2012, our results show that on lower taxonomic levels (like the genus level) shared ancestry may suffice to explain certain environment-phenotype associations. This point is particularly relevant in paleoanthropology where phylogenetic reconstruction is often complicated by the unresolved alpha taxonomy. The study of absolute and relative tooth size has been an important tool for the reconstruction of hominin diets and adaptive zones of closely related species (Kay, 1985;Organ et al., 2011;Robinson, 1954;Wood and Collard, 1999;and reviewed in Ungar, 2011). Our results offer a word of caution against adaptive interpretations of craniodental variation when the phylogeny of the studied taxa is either not known or not explicitly modeled in the analysis.