Sample, Measurements, and Repeatability

Skull measurements were obtained from 5222 crania of 16 genera and 110 species of Platyrrhini deposited at the following institutions: American Museum of Natural History, British Museum of Natural History, Field Museum of Natural History, Museu de Zoologia da Universidade de São Paulo, Museu Nacional do Rio de Janeiro, Museu Paranaense Emílio Goeldi, Museo de la Universidad Nacional Mayor de San Marcos, National Museum of Natural History, and the University of Tennessee, derived from the Marmoset Research Center, Oak Ridge Associated Universities' colonies. A complete list of examined specimens may be obtained from the authors upon request, and details of the taxonomic arrangement employed and sample sizes for each taxon are given in Marroig and Cheverud (2001). Only adult crania were used in the subsequent analyses. Specimens were considered adult when they had totally erupted and functional dentition as well as closed or fused spheno-occipital and/or spheno-ethmoid sutures.

Three-dimensional coordinates were recorded for 36 landmarks (Fig. 1) using a Polhemus 3Draw digitizer (Colchester, VT). The general procedure for measuring specimens and the landmark definitions follow Cheverud (1995). A set of 70 linear measurements describing cranial morphology was calculated from the coordinate values. This set was reduced to 39 measurements, averaging the measurements present on both sides of the skull (for the list of 39 traits see Marroig and Cheverud 2001, 2004). Repeatabilities for our measurements are generally above 0.90 when considering variation within a single population (Cheverud 1995, 1996; Ackermann and Cheverud 2000; Marroig and Cheverud 2001).

Overall Strategy, Summary, and Symbols

Tables 1 and 2 present a summary of the symbols used for the variables and analyses performed. Most analyses described below are based on the comparative method and make use of the NWM phylogeny to reconstruct ancestral estimates for each node in the tree. Thus, the overall strategy is to focus on phylogenetically independent changes and associations occurring within each internode (along branches) of the tree. In addition, we use a more familiar method to test for association among our variables using the independent contrasts (IC) to account for the nonindependence of phylogenetically structured data (Garland and Ives 2000). We use the module PDAP (Garland and Ives 2000) within the Mesquite package (Maddison and Maddison 2003) to obtain the correlation among variables under the IC approach.

Size Vector, Line of Least Resistance, and the Direction of Evolution

Pooled within-group phenotypic variance/covariance matrices (below referred to as covariance matrices or W for simplicity) were estimated for each genus using the general linear model (GLM) SYSTAT 11.0 (Systat Software, Inc., Richmond, CA) routine to control for sexual dimorphism and species differences whenever appropriate (for details see Marroig and Cheverud 2001). Averages for the 16 genera were obtained as least squares averages from the GLM module in SYSTAT with genera and species nested within genera as independent factors and the skull trait set (n = 39) as dependent variables. Throughout the paper we use vector correlation to measure the similarity of any two vectors. The vector correlation (Blackith and Reyment 1971) is a measure of vector orientation similarity in a p-dimensional space (p being the number of traits). The correlation between two vectors is equal to the cosine of the angle θ formed between them. The expected range of correlations commonly occurring among 39-element vectors by chance alone is −0.4 < r < 0.4 (Ackermann and Cheverud 2000), but only absolute values need to be taken into consideration because an angle of θ > 90° is equivalent to one of 180° − θ. We use a broken-stick model to obtain 1000 random 39-element vectors from a uniform distribution to then correlate each random vector to a fixed isometric vector (all elements equal 0.160). Our sample showed an average correlation of 0.127 and a SD of 0.095, both values being the same no matter which fixed vector was used for comparisons with the random vectors. These two statistics, based on a random sample of 39-element vectors, allow us to test whether the correlation of any two observed vectors is significantly different from the correlation between two vectors expected by chance.

Similarity in the p_max direction might be an indication of a shared LLER, but that does not mean that evolution necessarily follows that path. To test whether the LLER influenced the direction of morphological diversification, we need to compare the orientation of g_max to Δz. We first estimate ancestral states for our continuous trait set using a maximum likelihood (ML) approach (Schluter et al. 1997) implemented in the program ANCML available from D. Schluter's home page ( http://www.zoology.ubc.ca/~schluter/ancml.html). Although we use ML estimates, ancestral states reconstructed using maximum parsimony (MP) with the squared cost assumption (Maddison and Maddison 2003) result in virtually the same values. Tree topology and branch lengths were obtained from the same four-gene combined dataset used by Schneider et al. (2001). The Jukes-Cantor model was used to obtain the phylogenetic tree using the neighbor-joining method (Kumar et al. 2001). Our branch lengths (assumed to be roughly proportional to time) had a correlation of 0.993 with those reported by Schneider et al. (2001), with the largest difference being only 0.003 between the two datasets. After the estimation of ancestral data, the direction of evolution (Δz) within each branch was obtained simply as the difference vector in the 39 averages of each taxon pair. It is important to keep in mind that ancestral nodes are reconstructions and not true ancestors relying on assumptions of character evolution and on the accuracy of the phylogeny (Omland 1999). While the tree we are using is highly supported statistically (Schneider et al. 2001), uncertainty about the root of the NWM tree is a cause of concern (see Schneider et al. 2001). We deal with this uncertainty using the two alternative hypotheses for the rooting of the NWM tree and redoing all analyzes to check the consistency of the results.

We compare the direction of evolution with both versions of LLER (g_max and p_max). While the genetic line of least evolutionary resistance (g_max) is unavailable for most taxa used here, we do have an estimate for it in the genus Saguinus (Cheverud 1996). This estimate of g_max was obtained as the first principal component (PC1) of the pooled within-groups additive genetic variance/covariance matrix of the two tamarin species studied by Cheverud (1996). We correlate this g_max with the p_max of each genus. If p_max is a reasonable proxy to g_max we expect that those correlations should be high. Throughout the paper whenever we refer to g_max we are referring to the PC1 obtained from the Saguinus additive genetic W matrix.

The direction of p_max was calculated from the W matrices as the PC1 for each taxon. This procedure is straightforward in the terminal branches of the tree (genus) but a bit less certain for the deepest nodes. The W covariance matrices were obtained for each node in the Platyrrhini tree, pooling all descendent genera from that node to estimate W and structuring the pooling according to the phylogeny. For example, the W matrix of node 30 (Fig. 2) involves the pooling of the W matrices of the genera Callithrix and Cebuella. Then, the W matrix of node 28 was obtained pooling the node 30 W matrix with the Callimico W matrix and so on through the whole phylogeny. Note that we are not taking into account the branch lengths in the estimation of ancestral W values, just the topology of the tree is used to direct the pooling. This is justified by our previous finding of very similar W values among NWM and the independence of covariance matrix structural similarity relative to phylogenetic distance (Marroig and Cheverud 2001).

Finally, with both Δz and p_max (or g_max) estimated, we obtain the angle (θ) formed between them as a measure of how closely morphological diversification follows the LLER. When evolutionary change occurs primarily along the LLER, θ will approach zero and, conversely, as Δz deviates more from g_max, θ will increase. We used the cosine of the angle θ (vector correlation) to measure the relationship between Δz and g_max or p_max and therefore the smaller the angle the larger the correlation between vectors. Our strategy here was to focus on the evolutionary changes happening along each of the 30 branches on the tree and compare the direction of that change to the phenotypic and genetic lines of least resistance.

Lines of Least Evolutionary Resistance and the Amount and Pace of Morphological Change

Another question is whether the LLER influenced the amount and pace of morphological diversification. To answer this question, we first obtained a measure of the amount of morphological differentiation (m_amount) as the sum of squared differences of the averages for the 39 traits for each of the 30 internodes. A measure of the pace of morphological differentiation (m_pace) was obtained dividing the m_amount by the length of the corresponding branch on a genetic distance scale. We then compare m_amount and m_pace with the Δz-p_max and Δz-g_max using Pearson product moment correlation. If selection favors divergence in a direction markedly different from that of the LLER, the expectation is that the evolutionary amount and pace should be relatively small and slow respectively (Schluter 1996). The reason is straightforward: directions other than the g_max (or p_max) have less genetic variance, reducing the amount and rate of the evolutionary response expected along those other lines given the same intensity of directional selection. Therefore, a negative correlation of the amount and pace of morphological differentiation and θ (or a positive correlation if the cosine of the angle Δz-p_max is used) is expected.

Allometric Size and Absolute Size

It is important to keep in mind that size has two meanings here. Absolute size refers to the overall body size of an animal (or an estimator of it like the PC1 score), whereas p_max or allometric size refers to the allometry of one variable with respect to the size variation in a population of organisms of different absolute sizes, being therefore a measure of the relative contribution of each variable to absolute size. In other words, p_max is a measure of variation, whereas absolute size is a measure of the average skull (or body) size. Absolute skull size was obtained as the score on the PC1 extracted from the between-populations V/CV matrix L (L is the V/ CV constructed from the 16 living genera and 15 ancestral nodes population means). This between-taxon PC1 accounts for 92.8% of the total variation among groups. For the 16 living genera, the correlation between estimated skull size and their body weights (Fleagle 1999) is 0.987 (P < 0.0001), indicating that the PC1 scores are a good proxy for overall body size. Absolute skull size differences (s_dif) were obtained from the squared differences in the PC1 (extracted from L) scores of each pair of taxa. This measure (s_dif) is used here because is quite possible that large size differences are related to niche differentiation and time elapsed as well as with the LLER (Table 2). Allometric size (p_max) was compared in two ways. First, we estimated the pairwise vector correlations among the PC1 values of each living genus. This is a measure of similarity of the two vectors and an explicit measure among NWM genera of the existence and strength of the LLER and of their persistence through time when coupled with phylogenetic information (Maynard Smith et al. 1985; Arnold 1992; Schluter 1996). Second, we obtained PC1s for the W matrices of each node in the tree (p_max) and focus on the similarity of those PC1s structuring our comparisons at each of the 30 internodes on the tree in the same way we did for comparisons between Δz and p_max or g_max (see above).

Phylogeny and Diet

Diet similarity and phylogenetic distances were previously reported by us (Marroig and Cheverud 2001). In addition, we also estimate the diet of ancestral states using the same procedure (ML) described above for the quantitative traits. With the averages estimated for the five dietary variables (fruits, insects, gum, seeds, and leaves) we also obtained measures of the amount and pace of evolutionary change in diet for each branch in the tree. The amount of dietary change (d_amount) was obtained as the sum of the squared differences of the averages for the diet variables for each of the 30 internodes. Likewise, the pace of dietary change (d_pace) was obtained by dividing d_amount by the corresponding branch lengths. Finally, from the dietary similarity matrix we obtained a combined measure comparable to the absolute size measure (PC1). This was done by performing a multidimensional scaling analysis (MDS) on the dietary similarity matrix and using only one axis to represent this ecological variable. This allows us to test whether size and diet evolution are significantly associated.

Lines of Least Evolutionary Resistance Association with Size and Diet

Correlation between diet (d_pace and d_amount) morphological (m_pace and m_amount) and absolute skull size differences (s_dif) as well as their associations with Δz-g_max and Δz-p_max were obtained. Because the changes occurring within each branch are independent from the other branches, a standard Pearson correlation coefficient was used. Many of these variables were nonnormal and therefore some transformations were applied before analyses. Fisher's transformation was used to normalize the Δz-g_max and Δz-p_max correlation values, and a square-root transformation was used to normalize the diet pace and amount variables as well as branch lengths. Morphological pace and m_amount were transformed by obtaining their natural logarithms. Skull size difference (s_dif) was normalized by multiplying the original values by a constant (10,000) and then obtaining its natural logarithm. To correct for multiple testing, we use the approach outlined in Cheverud (2001) to estimate the number of effective independent comparisons and then applied the Bonferroni correction to adjust the alpha level.

Diet can be conceived as being organized into broad adaptive zones in NWM (Figs. 2, 3). The invasion of new diet-based adaptive zones occurred four times during the Platyrrhine diversification. All NWM eat fruits, allowing for variation in the proportion consumed and of the other items used to supplement their diets. Each of the four other diet categories (seeds, insects, gum, and leaves) define four distinct clades within Platyrrhini. If the invasion of new dietary adaptive zones was one key factor unleashing adaptive radiation along LLER in Platyrrhini, we expect that the largest transitions in absolute size as well as the largest Δz-p_max and Δz-g_max correlations to occur in those branches associated with the invasion of new diet zones. Invasion of new dietary zones occurred at the following internodes (see Fig. 2): 10– 3 (the seed zone), 10–11(leaf zone), 22–24 (gum zone), 2– 18 (insect zone), and 2–10 (leaf-seed). We used one-tailed t-tests with unequal variances to detect whether the values at these internodes were significantly different from those of the others.

Retrospective Selection, g_max, and the Direction of Evolution

The appraisal of the similarity between Δz and g_max as evidence for genetic constraints rests on the assumption that the direction of divergent natural selection is random with respect to the direction of g_max (Schluter 2002). An alternative hypothesis is that selection, rather than genetic constraints, is responsible for the bias in the direction of divergence (Schluter 2002). In such a scenario, the direction of evolution (Δz) is closely aligned with g_max because selection corresponds to g_max. The evolutionary response of a set of quantitative traits is described by the equation Δz = Gβ, where Δz is the vector of differences in means between generations, β is the selection gradient vector, and G is the additive genetic V/CV matrix (Lande 1979). Rearranging the evolutionary response equation, the pattern of selection responsible for species differences can be reconstructed from observed mean differences using the following relationship:

The following relationship was used to obtain the selection gradients:

Allometry, Line of Least Resistance, and Evolutionary Divergence

NWM genera share the same basic allometric patterns (p_max) of variation in skull morphology. This allometric vector represents the line of least evolutionary resistance. Given the considerable differences in absolute size among NWM and their deep history (at least 30 million years of evolutionary diversification), there is surprising conservation of shape variation associated with size. Allometric size accounts for 20–40% of all variation within each taxon. Differences in absolute skull size, as measured by the between genera PC1, accounts for 92.8% of the total variation among genera and is very similar (total p_max in Table A1, available online) to the within genus PC1s (vector correlation X̄ = 0.86, SD = 0.066). This simple result suggests that morphological divergence between genera follows the within-population allometric line. The conservation of allometric patterns is even more striking when we consider that it is not only a phenotypic pattern but also reflects the genetic architecture underlying it. This is reflected in the correlation between g_max and p_max (X̄ = 0.852, SD = 0.051) uncorrected in this case for sampling error, but still fairly high. Given the striking similarity in g_max and p_max found here, we consider both interchangeable for the purpose of the discussion below (Table 3). A genetic LLER therefore appears to exist in NWM and to have persisted for at least 30 million years. But did that LLER influence the paths of evolutionary diversification for NWM crania?

An inspection of the Δz-p_max in Figure 2 (see also Table 3) reveals that most of the diversification in NWM follows the LLER defined by size, with the Δz-p_max correlation usually in the range of 0.85 to 0.97. Moreover, the high correlation between absolute size differences and m_amount also shows that most of the morphological diversification in NWM was related to size diversification (a result also consistent with the high correlation between total p_max and the p_max values, Table A1, available online). While in NWM these two measures of morphological differentiation are very highly correlated, there is no a priori reason why they should be. Keep in mind that m_amount is a measure of how much two groups differ in overall morphology (not only size), whereas s_dif is a measure of the absolute size differences between two groups. This result is robust to alternative topologies of the phylogenetic tree (results not shown), where the root is placed in the Atelinae or in the Pithecinae clades before reconstruction of the ancestor states, supporting that evolution followed p_max or g_max (Table 3) during the NWM radiation.

However, there are also many instances where the direction of evolutionary change is remarkably distinct from the size dimension. Examples are Leontopithecus, Lagothrix, Chiropotes, and Pithecia, all showing values of Δz-p_max close to or below that of random expectation (Figs. 2, 5). These four cases are exactly those where the size differences between ancestors and descendents are the smallest (Fig. 6) observed in the tree (t = 3.217, df = 3.8, P = 0.034). Therefore, whenever evolution proceeded in a direction remarkably different from the LLER, the absolute size differences were also comparatively small. It is also interesting to note that those four cases with lower Δz-p_max values also had a slower pace of morphological change (m_pace, t = 4.019, df = 7.6, P = 0.002) compared to the remaining cases, but while the amount of morphological change on average is smaller in those four genera (m_amount, t = 1.672, df = 18.7, P = 0.055) this last result is not significant. The difference in these two results is basically due to node 22, which has a very small amount of morphological change relative to its last common ancestor (node 18) but has a very high pace of morphological change due to its very short branch length. Excluding internode 22– 18 from the comparison, the t-test results are: m_amount, t = 2.238, df = 15.9, P = 0.020 and m_pace, t = 4.156, df = 7.7, P = 0.002. These results suggest that when evolution does not occur along the LLER, the amount and pace of morphological change are small and slow respectively, which is also apparent in the association analysis results presented in Tables 4 and 5. This is in agreement with the expectation advanced by Schluter (1996) that evolution should be slower when selection favors divergence in a direction markedly different from that of g_max. Therefore, taking these results together, there is support for the view that the LLER influenced the path, rate, and amount of morphological evolution in NWM at a macroevolutionary scale.

We now address two final and related questions about the LLER: How much diversification is allowed along such lines and for how long do those lines persist? Conceivably, natural selection can change the patterns of morphological integration either by actively building new adaptive trait ensembles (Berg 1960; Baker and Wilkinson 2003; Bégin and Roff 2003) or indirectly while changing the averages of evolving populations and consequently affecting correlation among traits (Marroig and Cheverud 2001; see also Jones et al. 2004). By altering morphological integration patterns, selection should also change the lines of least resistance. Therefore, in an adaptive radiation it is reasonable to expect that at some point the lines of least resistance will diverge in orientation. Schluter (1996) predicted that the bias in evolutionary direction should be temporary and diminish with time. Our results here are mixed. First, there is indeed a trend of longer branches being associated with lower Δz-p_max correlations (Fig. 5, Table 4), but those correlations are not significant. Also, there is a negative correlation of the p_max-p_max vector correlations (Table 3) and their associated branch lengths (Fig. 7), indicating that p_max itself is diverging with time. One possible explanation for such a pattern is that any correlation between these vectors (Δz-p_max or p_max-p_max) at the tips of the tree is compounded at deeper nodes given our strategy of reconstructing ancestral states. However, a comparison of deeper nodes against tips in the tree after excluding the outliers discussed above (see also Figs. 5, 7) show no significant difference between the two groups with respect to the similarity in p_max orientation, suggesting that compounded error is not responsible for this pattern. Moreover, using the IC approach to compute the correlation between Δz-p_max (Table 5) or p_max-p_max (R = −0.139, P = 0.605) with branch length reveals no association among them. Conversely, longer branches are concentrated at the tips of the NWM tree (t = −5.28, df = 27, P < 0.0001), which might well explain why deeper nodes present higher correlations (less time to diverge). Second, while this trend occurs within branches, on a macroevolutionary scale the g_max and p_max conservatism discussed above suggests that allometric size acts as a LLER during the whole period of the NWM radiation, with no trend to diminish with time's arrow. This conclusion is also reinforced by the low and nonsignificant matrix correlation followed by a Mantel test (R = 0.013, P = 0.552) observed between phylogenetic distances and similarity in p_max (Table A2, available online). The amount of cranial morphological diversification occurring in NWM along the LLER might be quantified as the relative difference in skull size between the extremes of variation observed in Figure 2. Brachyteles' skull is 3.57 times bigger than the skull of Cebuella, which corresponds to an 81-fold difference in body mass (Fig. 2). Another measure of the amount of morphological diversification occurring in NWM along the LLER is the between groups PC1 eigenvalue, which inform us that 92.8% of the total variation is due to size changes (given the similarity between total p_max and the p_max values). Therefore, in NWM there is weak support for the prediction that the bias in the direction of evolution imposed by the LLER should be temporary and diminish with time. Our results indicate that not only the LLER endured on a macroevolutionary scale but also allowed for a large amount of diversification without changing its orientation.

Adaptive Radiation, Size Evolution, and Diet

Adaptive radiation is the evolution of ecological and phenotypic diversity within a rapidly multiplying lineage (Schluter 2002). Four features might be used to detect an adaptive radiation according to Schluter (2002): common ancestry, phenotype-environment correlation, trait utility (fitness advantage), and rapid speciation. All NWM share a common ancestor and we have indirect evidence for cranial trait utility (Marroig and Cheverud 2004). The indirect evidence came from the analyses of morphological evolution in NWM showing that natural selection dominated their diversification, particularly at the levels of the origins of families, subfamilies, and genera, and especially upon size (Marroig and Cheverud 2004). Rapid speciation occurred when the three NWM families emerged nearly simultaneously (Schneider et al. 2001) at the Oligocene–Miocene boundary (26 million years ago) giving rise to three of the four major diet adaptive types in a burst of diversification (Atelidae, leaves; Pitheciidae, seeds; Cebidae, insects). The correlation between form and function is not totally established, at least not quantitatively, because association of skull morphology and environmental factors (like diet) are still sparse for NWM (Rosenberger et al. 1996). However, the association of diet with skull morphology is clear in our results (Fig. 3, see also Marroig and Cheverud 2001). Considering these four features, the evolution of NWM might be considered an adaptive radiation resulting today in more than 110 living species (Rylands et al. 2000) occupying many adaptive zones related to diet, locomotion, systems of mating, and habitats.

We suggested before that the adaptive radiation in NWM was triggered by an association of size evolution with dietary diversification (Marroig and Cheverud 2004). This association is apparent in our results here in the significant correlation between size and diet (Fig. 3). The largest NWM are leaf-eaters and the smallest are gum-eaters (Fig. 3). Tables 4 and 5 also show this relationship with size differences (s_dif) and morphological amount and pace of evolutionary change (m_amount and m_pace) being associated with both dietary change measures (d_amount and d_pace), suggesting that larger dietary differences are associated with larger size differences. The association of size and diet evolution is also apparent in the t-tests comparing the new invasion and stable groups. Significant differences were observed for size difference (s_dif), Δz-g_max and Δz-p_max, indicating that indeed the largest absolute size transitions are associated with the invasion of new dietary adaptive zones and that the direction of evolutionary change was close to the LLER in those invasions.

An empirical pattern found by Kay (1984) looking at the distribution frequency of primate species body sizes and diet is that dietary habits are correlated with body size (Fleagle 1999). Species that eat insects tend to be relatively small, whereas those that eat leaves tend to be relatively large (Kay 1984). But why should diet and morphological evolution, especially size, be associated? The answer involves the relationships between absolute size changes, ecological opportunity related to empty diet niches in the past, and metabolic, physiologic, and foraging issues related to scaling (for a thorough account of the ideas discussed below, see Fleagle 1999, ch. 9, and references therein). Perhaps the easiest way to visualize these relationships is to consider a simple question: Why don't we have large primates feeding predominantly on insects, or tiny vegetarian primates eating mainly tree leaves? All NWM are fruit eaters that supplement their diets with other items to meet the needs of a balanced diet in terms of energy and other nutritional requirements, such as proteins and trace elements (Fleagle 1999). Fruits are a good source of calories but are very low in protein content. To supplement their diets two abundant resources are usually exploited by primates: other animals (such as insects) and folivorous material (leaves, shoots, and buds). Insects (and animal material in general) are an excellent source of nutrients and calories per unit weight, but they are difficult to catch and the numbers of prey items ingested during an active period should be dependent on their local abundance. Small primates can live by preying on insects quite effectively but a large animal will not be able to capture enough insects to meet its needs for protein and energy. Unlike insects, leaves are neither cryptic nor difficult to catch, but instead pose other problems for their consumers. While a good source of protein, leaves are generally low in energy yield for their weight compared to fruits and insects. Consequently, large amounts must be ingested and processed by the gut. Leaves are also composed of large amounts of less digestible components, like cellulose or even toxins. A larger body helps a primate overcome these problems inherent to a leafy diet providing longer gut tubes, which allows proper digestion and detoxification. Thus, while the lower size limit of folivorous primates seems to be determined by metabolic and digestive parameters, the upper size limit of insect eaters seems to be imposed by the time to locate and catch prey (Fleagle 1999). If we consider that at some point in the past, close to the beginning of the adaptive radiation of NWM, there were empty dietary niches available for invasion, size evolution is a consequence of adaptation to those niches. We can imagine that the empty dietary niches correspond to peaks in the adaptive landscape and ascending those peaks involves changes in absolute size. Thus, the ecological force driving NWM adaptive radiation is dietary diversification or, under a Simpsonian view of nature (Simpson 1953), the invasion and further subdivision of diet-based adaptive zones. Natural selection operated to diversify diet in NWM, allowing species to exploit previously unexplored adaptive zones (Van Valen 1971). The cranial morphological correlates of that adaptation were such that genera diverged along a cranial size dimension that is also the dimension of greatest variance within species (LLER).

Selection, Constraints, or Both?

One final question regards the role of selection versus constraints during the adaptive radiation. If size evolution resulted from selection operating to diversify diet in NWM, was the evolutionary change a result of the morphological integration patterns biasing the evolutionary divergence (Δz) in the direction of maximum additive genetic variance (g_max) or did selection operate directly upon both g_max and Δz? In other words, was the close link between g_max and the direction of evolution evidence of long-term genetic constraints or were both g_max and size under selection? Our results are surprising in revealing that the selection gradients (β) are not similar to either p_max (or g_max) or Δz values (Fig. 4). This suggests that the similarity observed between evolutionary change in NWM and the axis of maximum genetic variance (g_max) arises from the operation of long-term constraints. These constraints might be the result of stabilizing selection operating via a common developmental program that structure or modulate the expression of genetic and environmental variability (Marroig and Cheverud 2001). Therefore, our results so far suggest that the genetic correlation structure shapes macroevolutionary patterns by directing the evolution of mean phenotypes toward certain adaptive peaks and away from others, regardless the selection patterns operating (the constraints models in Björklund and Merilä 1994; Baker and Wilkinson 2003).

However, we also took another approach to test the adaptive model (Björklund and Merilä 1994; Baker and Wilkinson 2003). To test whether selection was operating primarily on size, we built an isometric hypothetical selection gradient (all 39 elements equal to 0.160 or 1/√39). This size selection vector (β_size) was then multiplied by the W matrices to give the hypothetical evolutionary response to selection on size alone (Δz_h) for each internode in the NWM tree. We then compare those Δz_h values with the observed Δz values using vector correlations. The hypothetical and observed evolutionary responses tend to be quite similar (X̄ = 0.77, SD = 0.22; Table 3). Moreover, there is a striking similarity between the Δz-Δz_h and the Δz-p_max (R = 0.988, P = 1 × 10⁻¹⁵), indicating that selection on size alone would produce evolutionary responses along p_max. The exceptions are the same as those observed in Figure 2 and represent those cases where evolution was not size based. A full exploration of these approaches would be to randomly generate a subset of allometric size selection vectors (not only isometric), apply them to the W matrices and then compare the Δz values produced with the observed ones and to a full set of responses produced from random selection vectors (not only allometric or isometric size). Because dietary diversification in primates is linked with body size changes, we expect that empty dietary adaptive zones would impact covariance patterns as a size selection gradient. Therefore, considering all the results presented here, it is perhaps not useful or possible to separate the influence of genetic patterns of covariation from the influence of natural selection on the evolution of NWM cranial diversity.

While it seems clear that NWM cranial diversification occurred primarily along the LLER, it is not clear whether that path was followed because of selection along it or because evolution was constrained to follow this path. Even though the general pattern is for evolution to follow the path of least resistance, several instances of diversification along other, distinct, morphological dimensions were also found (Leontopithecus, Lagothrix, Chiropotes, and Pithecia) indicating that constraint is not absolute in NWM. It is also clear that the path of least resistance is very conservative in the NWM, remaining consistent over 30 million years of evolutionary change. Genera that diverged in directions distinct from the LLER do not display especially distinctive patterns of within species allometry (Tables A1, A2 available online) or variance/covariance matrices (Marroig and Cheverud 2001), and thus their morphological diversification cannot be attributed to changes in constraints. Also, only 20–40% of the within-species variation is due to size (PC1). Clearly, there is substantial genetic variation in other dimensions for response to selection. Thus, it seems most likely that size was often the target of selection in NWM evolution and that cranial morphology evolved allometrically in response (Lande 1979; Zeng 1988).