The evolutionary relationship among beak shape, mechanical advantage, and feeding ecology in modern birds *

Extensive research on avian adaptive radiations has led to a presumption that beak morphology predicts feeding ecology in birds. However, this ecomorphological relationship has only been quantified in a handful of avian lineages, where associations are of variable strength, and never at a broad macroevolutionary scale. Here, we used shape analysis and phylogenetic comparative methods to quantify the relationships among beak shape, mechanical advantage, and two measures of feeding ecology (feeding behavior and semiquantitative dietary preferences) in a broad sample of modern birds, comprising most living orders. We found a complex relationship, with most variables showing a significant relationship with feeding ecology but little explanatory power. For example, diet accounts for less than 12% of beak shape variation. Similar beak shapes are associated with disparate dietary regimes, even when accounting for diet‐feeding behavior relationships and phylogeny. Very few lineages optimize for stronger bite forces, with most birds exhibiting relatively fast, weak bites, even in large predatory taxa. The extreme morphological and behavioral flexibility of the beak in birds suggests that, far from being an exemplary feeding adaptation, avian beak diversification may have been largely contingent on trade‐offs and constraints.

beak phenotypes might be retained if they are efficient for processing nonfavored resources, particularly if the favored resource is periodically limited (i.e., Liem's paradox;Liem 1980;Tebbich et al. 2004). Furthermore, in addition to feeding and foraging, birds use their beaks for a plethora of other tasks, such as preening (Moyer et al. 2002;Clayton et al. 2005), vocal modulation (Podos 2001;Herrel et al. 2009), thermoregulation (Tattersall et al. 2009;van De Ven et al. 2016) and water balance (Greenberg et al. 2012), tool use (Weir et al. 2002;Wimpenny et al. 2009;Laumer et al. 2017), nest construction (Hansell 2000), and as a display structure (Navarro et al. 2009). This functional and behavioral flexibility implies that multiple selective pressures likely played important roles in shaping beak evolution. Understanding the relative importance of trophic adaptation to beak morphological diversification in modern birds is therefore vital to understanding avian evolution, and to make accurate ecological inferences in extinct taxa (Lauder 1995;Rubega 2002).
Although the main patterns of beak shape evolution at a broad macroevolutionary scale in birds have been effectively characterized (Cooney et al. 2017), the extent to which such patterns are related to feeding ecology, or to biomechanically relevant traits such as the mechanical advantage (MA) of the jaws, remains largely unexplored. Besides Darwin's ground finches (e.g., Grant and Grant 2006), quantitative evidence evaluating the link between feeding ecology and beak shape in birds is limited to a handful of avian clades (Rubega 2002). These few studies have found strong associations in several families of passerines (Gosler 1987;Benkman 1988;Price 1991;Peterson 1993;Bardwell et al. 2001), anseriforms (Olsen 2017), and a few charadriforms (Barbosa and Moreno 1999), but weak associations among birds of prey (Bright et al. 2016). Biomechanical modeling is similarly limited taxonomically, but in Darwin's finches, it has been shown that skull and beak shapes are adapted to the mechanical demands of feeding (Soons et al. 2010;Soons et al. 2015).
Here, we use geometric morphometrics (GM) to quantify beak shape variation and its relationship with feeding ecology in a broad sample of birds. Shape analysis based on GM provides the analytical tools to partition the sources of beak shape evolutionary variance, as well as to test the strength and pattern of correlation with independent variables (Monteiro, 1999;Rohlf and Corti, 2000;Marugán-Lobón et al., 2013). Ecology is characterized by three components of feeding: we quantify the MA of the jaws as a functional trait related to the ability to transfer force or movement through the skull system (high MA describes efficient force transfer, low MA defines less efficient force transfer but faster jaw movement) (1); tabulate biological role by documenting use of the beak during feeding (2); and recompile detailed semiquantitative dietary data (3) for each of the studied species. We use multivariate statistics and phylogenetic comparative methods to test for correlations between these variables, while also accounting for the effect of size (i.e., evolutionary allometry) on beak shape, force transfer, and diet.

PHYLOGENETIC HYPOTHESIS
Our study includes 175 species from 94 families of extant birds, encompassing 38 of the 39 living orders, excluding only Mesitornithiformes, a Madagascan clade of three species (Hackett et al. 2008;Del Hoyo et al. 2017) (Table S1). A maximum clade credibility phylogeny of the 175 species was generated using TreeAnnotator (Rambaut and Drummond 2013) from a population of 10,000 "Hackett's backbone 'stage 2' trees" downloaded from www.birdtree.org (Jetz et al. 2012) (Fig. 1). Branch lengths were set equal to 'Common ancestor' node heights.
The feeding autecology (the presumed main biological role of the beak) of each species was characterized using two sources of ecological information, namely, semiquantitative dietary preferences and the use of the beak during feeding (UBF) (Fig. 1). The dietary data for each species were sourced from EltonTraits 1.0 (Wilman et al. 2014). This data was coded as a matrix of estimations of the relative importance of 10 main dietary categories translated from species-level dietary descriptions in the literature (Fig. 1, Table S1) to the overall diet of each species. These estimations were coded as bins of 10 units of percentage (i.e., 0, 10, 20, 30 . . . 100%) (Fig. 1, Table S1). A detailed description of the specific food items included in each category is included in the metadata archives in Wilman et al. (2014). To obtain a Euclidean representation of this noncontinuous data, we calculated a symmetric similarity/distance matrix (Euclidean distances) from the original 175 (species) × 10 (dietary items) matrix to conduct Principal coordinates analysis (PCoA) in PAST v.3.15 (Hammer et al. 2009) and used the scores from the PCoA for downstream analyses (following Legendre and Anderson 1999).
The use of the UBF was categorized by applying a simple dichotomous key (Fig. S1) to published observations of foraging and feeding behavior of each of the studied species (Del Hoyo et al. 2017). This allowed us an alternative means to subdivide feeding autecology given that dietary categories at such a wide phylogenetic scale often include very different foraging and feeding behaviors. For instance, the Atlantic puffin (Fratercula arctica) and the osprey (Pandion haliaetus) both feed almost entirely on fish (Wilman et al. 2014;Del Hoyo et al. 2017), but while the former feeds by underwater pursuit diving and grabs individual fish directly with the beak, the latter plucks fish from the water with the talons, and uses the beak instead to tear off chunks of meat before consumption (Del Hoyo et al. 2017). The UBF categories for these examples are therefore scored as "Grabbing/gleaning" and "Tearing" (Table S1). Every species in our dataset except the American flamingo (Phoenicopterus ruber, a specialized filter feeder) fits in to one of five categories (tearing, cracking/ biting, pecking/grazing, grabbing/gleaning, and probing; Fig. 1;  Fig. S1).

BEAK SHAPE AND SIZE
The skull of each species (without the rhamphotheca, the corneal sheath that covers the bony beak, which is commonly missing in museum specimens) was photographed in lateral view (Table S1), and the complete outline of the beak was digitized using a set of three fixed landmarks and two curves (Fig. 2), the latter comprising 50 evenly spaced semilandmarks (25 along the dorsal profile of the bill (culmen), and 25 along the left dorsoventral  Figure S1). Numbers correspond to clades as detailed in Table S2. edge of the beak (tomium)). The landmarks and semilandmarks were digitized in tpsDig2 (Rohlf 2006). The Minimum Bending Energy sliding method (Bookstein 1996, Bookstein 1997) was used to slide the semilandmarks in tpsRelw (Rohlf 2010), as this is more reliable when morphological variation is large (Perez et al. 2006;Fernández-Montraveta and Marugán-Lobón 2017). Shape data (i.e., Procrustes coordinates) were extracted using a full Procrustes fit and imported to MorphoJ (Klingenberg 2008), PAST v.3.15 (Hammer et al. 2009), and the R package geomorph v.3.0.6 , in which all the subsequent analyses were performed. Preliminary analyses revealed that slender, straight beaks are consistently associated with higher values of log-centroid size (CS; Fig. S3, Table S5). This is undesirable as it may erroneously exaggerate allometric effects particularly when variance is very skewed towards one direction, impeding our ability to reliably test for allometry using CS (Bookstein 1991). Beak allometry was therefore assessed using species average body mass (BM) data taken from Wilman et al. (2014). MA is a metric derived from lever mechanics (e.g., Uicker et al. 2011) and a well-established functional trait describing the tradeoff between bite force transmission and jaw closing speed during biting in vertebrates (e.g., Westneat, 1994;Anderson et al. 2008;Sakamoto 2010). Given the same force input, a high MA indicates a relatively more forceful bite; low MA indicates a relatively less forceful but faster bite. MA is calculated as the ratio of the length of the in-lever divided by the length of the out-lever (Uicker et al. 2011) and was determined for each species' skull at two different bite points (Fig. 2). The in-lever arm here is defined as the orthogonal distance from the mandibular articular facet of the quadrate (the fulcrum) to the intersection point with the midline of the fossa temporalis between the postorbital and zygomatic processes of the skull, where the midline of the adductor mandibulae group lies, which is the main adductor muscle group in modern birds (i.e., m. adductor mandibulae externus medialis/superficialis (m. AMEM/S), Sustaita 2008; Lautenschlager et al. 2014) (Fig. 2). The out-lever arms are defined as the linear distance from the articular facet of the quadrate to the tip of the bony beak (i.e., landmark 1; anterior out-lever) or to the midpoint on the tomial curve bisecting landmarks 1 and 3 (posterior out-lever; Fig. 2). This approximates the mechanics of avian jaw closure as a 2D, third-class lever system, although the three-dimensional lever system is often more complex than this (Olsen and Westneat 2016). Lever arm measurements were taken for each species using ImageJ (Rasband 1997). As anterior and posterior MA values (as defined here) show a strong correlation (Fig. S2), for simplicity we only used anterior MA for all the subsequent analyses.

STATISTICS
A Principal Component Analysis (PCA) on the Procrustes shape data was performed in MorphoJ to explore the main patterns of beak shape variation. We mapped the phylogeny onto the PCA scores in MorphoJ using the weighted squared change parsimony method (Maddison 1991) to visualize changes in beak shape along the phylogeny (i.e., in the terminals and internal nodes). The phylogeny was also mapped over the anterior MA values to visually explore the changes in MA in MorphoJ using the weighted squared change parsimony method. Anterior MA values were also mapped as isoclines over the PC1-3 phylomorphospace plots using the software MATLAB (Grant et al. 2008).
We used phylogenetically informed (phylogenetic generalized least squares [PGLS]) regressions to test for potential correlations between our trophic data, MA, size, and beak shape variation using the R package geomorph v.3.0.6 ). Specifically, we tested six pairwise relationships ( Fig. 2): (1) beak shape variation and log-BM, to test if beak shape variation is allometric; (2) MA and log-BM, to test if MA variation is allometric; (3) the relationship between beak shape and MA; (4) the relationship between beak shape and dietary preferences; (5) the relationship between MA and dietary preferences; and 6) the relationship between BM and dietary preferences. PGLS regressions with dietary preferences as the independent variables also included UBF categories as a factor to account for the complex relationship between the dependent variables (i.e., beak shape, MA, and log-BM), dietary preferences (i.e., matrix of diet), and feeding behavior (i.e., UBF categories).
Phylogenetic MANOVAs were conducted in the R package geomorph v. 3.0.6 to test for pairwise differences in: (1) beak shape; (2) MA; and (3) BM between UBF group means. Because our variables are unevenly dispersed across our phylogeny (e.g., specialized piscivorous taxa belong mostly within particular clades, Fig. 1), which can severely reduce statistical power of linear models , we used randomizing residuals in a permutation procedure (10,000 iterations implemented in geomorph v.3.0.6  to assess statistical significance for all PGLS regressions and Phylogenetic MANOVAs, as this has been shown to be more robust to group-clade aggregations . Furthermore, because dietary preferences and UBF categories covary with each other (R 2 = 0.05547, F = 1.9848, Z = 2.2061, P = 0.023; e.g., taxa that use the beak for tearing tend to consume a higher percentage of vertebrates [e.g., raptors], Fig. 1), we used type II (conditional) sums of squares to assess the statistical significance of those PGLS linear models including both dietary preferences and UBF groups .
Current implementations of PGLS regressions assume a Brownian Motion mode of evolution. To test if our data meet this requirement, we compared the relative fit of the estimated residuals of shape, MA, and BM to three different models of evolution: Brownian Motion, Ornstein-Uhlenbeck, and Early-Burst. We used the residuals of the PGLS linear models conducted in this study and the AICc criterion to ascertain which model best fits the data in each case (the one yielding the lowest AICc value). For shape data, fitting these models requires reducing its dimensionality, therefore we used the first nine PCs (accounting for ß99% of the variance in all the PGLS models in which shape is the independent variable). Brownian Motion is only preferred over the other models in the PGLS model of MA as a function of diet. For the remaining PGLS models, the Ornstein-Uhlenbeck model is preferred, with only a small difference in AICc value in all the cases (except for the two PGLS allometric models, which are either nonsignificant or significant but explain little shape variance in our sample; Table 1; Table S6). We therefore interpret that our data do not greatly deviate from a Brownian Motion model of evolution, and thus meet the expectations of the PGLS linear models. Nevertheless, these results must be taken cautiously, as recent research suggests current modelfitting methods based on maximum likelihood are prone to exhibit ill-conditioned covariance matrices that could lead to errors of interpretation (Adams and Collyer 2017). The implementation of more complex evolutionary models for analyses of high dimensional data is not fully developed (Monteiro 2013), therefore it is a methodological endeavor that goes beyond the scope of this article.
Variation along shape vectors is displayed as thin-plate spline deformations of an outline diagram based on the lateral beak outline of the plush-crested jay (Cyanocorax chrysops, Corvidae, Passeriformes), the species which is most similar to the Procrustes mean. The coefficients from the PGLS regressions with shape as the dependent variable were used to calculate the beak shape differences along the regression vectors. The R code used for all the analyses is provided in the Supplementary Materials.

ALLOMETRY
The first three principal components (PCs) explain 92.54% of the total shape variance in our sample, implying that few dimensions underlie beak shape variation. The main axes of beak shape recovered in this study (Fig. 3 & 4, and Supplementary Materials) are roughly equivalent to those recovered by a crowdsourced study encompassing the 3D beak shapes of more than 2000 species of modern birds (Cooney et al. 2017), suggesting that discarding the third dimension and rhamphotheca produces comparable patterns of avian beak disparity at this macroevolutionary level. Namely, our PC1 describes the same lateral shape change (thin and straight, to deep and down-curved). Similarly, our PC2 (thin and curved, to deep and straight) and PC3 (downcurved to slightly upturned) explain similar shape changes to Cooney et al.'s PCs 2 and 4. Although some groups of birds cluster within restricted areas associated with deeper and curved beak shapes (e.g., Accipitriformes, Strigiformes, Falconiformes, and Psittaciformes), several species or clades widely diverge from their sister groups to different areas of the PC space (e.g., Semnornis, Piciformes; Podargus, Caprimulgiformes; Phoenicopterus, Phoenicopteriformes; the family Anatidae) or to cluster within the deep and curved scatter (e.g., Carduelis, Passeriformes; Musophaga, Cuculiformes; Figs. 3 & 4). PGLS regression of beak shape on log-BM is not significant (P = 0.362) (Table 1, Fig.  S4) revealing that beak shape allometry across birds as a whole is negligible.
MA  Table S1). However, MA values are generally low, and 80% of the taxa possess anterior MA values <0.14 (Figs. 3B & 5, and Table S1). Plotting MA over the PC1-3 space (Fig. 3A) reveals a broad trend between shape and MA: low MA values in positive PC1 (thinner, straighter beaks) and higher MA values in negative PC1 (deeper, more curved beaks). However, the trend is not linear, and there are islands of high MA, meaning that two taxa separated by small Procrustes distances may have quite different MA values. This biomechanical decoupling is particularly noticeable between tearing (i.e., mostly raptors) and cracking birds (i.e., mostly parrots). For instance, the boreal owl (Aegolius funereus, Strigiformes) and the hyacinth macaw (Anodorhynchus hyacinthus, Psittaciformes) show a Procrustes distance of only 0.073 between their beak shapes but they show extremely different anterior MA values (Fig. 3). Anterior MA values show a significant but weak (R 2 = 0.03479, P = 0.014) correlation with BM (Table 1; Fig. S4).
Although MA data show a statistically significant phylogenetic structure (P < 0.0001), most internal nodes are constrained to a narrow range of relatively low MA (Fig. 3B). Only two lineages clearly diverge from this: parrots (Psittaciformes), which explore more than half of the upper range of MA values; and sandpipers, snipes, and phalaropes (Scolopacidae), with extremely low MA values (Fig. 3B & Table S1). Some pheasants (e.g., Perdix) also exhibit high values of MA within the range of Psittaciformes, along with some specialized cracking/biting passerines such as the Northern cardinal (Cardinalis cardinalis) (Figs. 3 & 5). Clustering near the Psittaciformes with lower values of MA are mainly herbivorous taxa such as the snow goose (Chen caerulescens), the common linnet (Carduelis cannabina), the Western capercaillie (Tetrao urogallus), and the least seedsnipe (Thinocorus rumicivorus), as well as the Andean condor (Vultur gryphus). The latter represents a clear deviation from the general low MA values of Accipitriformes (Figs. 3 & 5), due to a ventral deflection of the beak tip that shortens the out-lever of New World vultures (Cathartidae) relative to the Old World vultures (Accipitridae).
PGLS regression of beak shape on anterior MA values exhibits a significant (R 2 = 0.133, P < 0.0001) correlation (Fig. 5). The shape differences described by this regression vector are remarkably similar to those described by PC1: thin, straight, long beaks (positive PC1) show the lowest values of MA, whereas deep, curved beaks (negative PC1) show the highest. Deviating from this general trend with much lower values of MA than predicted by the regression is the majority of the tearing group, composed of the Accipitriformes; the northern crested caracara (Caracara cheriway, Falconiformes); and Strigiformes  1). The remaining Falconiformes cluster closer to parrots than to other raptors exhibiting higher values than the rest of raptors (Fig. 5).

BEAK SHAPE AND FEEDING ECOLOGY
PGLS regression of beak shape as a function of dietary preferences and UBF revealed a significant but weak correlation between beak shape and overall dietary habits (R 2 = 0.1156, P = 0.001; Table 2). The effect of UBF groups in beak shape variation is also statistically significant but the correlation is not strong (R 2 = 0.0923, P = 0.001) ( Table 2). Such results are largely congruent with visual inspection of the PC1-3 plot, where the main dietary groups overlap without any clear separation, and UBF groups exhibit only slightly clearer regionalization (Fig. 4). For instance, tearing and cracking/biting birds tend to occupy the same areas of the morphospace, being restricted to deep and curved shapes in the negative extreme of PC1 (Fig. 4). Probing birds are restricted to the positive side of PC1, exhibiting relatively thin and straight shapes. Pecking/grazing taxa are restricted to approximately 0.0-0.1 on PC3, exhibiting relatively straight and flat beaks (Fig. 4). However, Phylogenetic MANOVA shows vertebrates 100% of the time and is scored therein as "VertFishScav"; here, it was rescored as "VertEnd" (Table S1).   Thin straight beaks tend to be associated with a higher percentage of invertebrate consumption in birds, and deeper curved beaks are associated with consumption of more mechanically demanding food items such as vertebrates and seeds (Fig. 6). Thin and slightly curved beaks are also associated with highly piscivorous taxa (Figs. S5 & S7), which together with visual inspection of shape vectors associated with other axes of dietary variations underlines that similar beak shapes are associated with disparate dietary regimes (Figs. S5 & S7). Furthermore, regressions show that the relationship between beak shape and dietary preferences differs between UBF groups (  while there are diet-dependent allometric relationships in our data, these are not affected by UBF behavioral groups (Table S4).

ECOLOGY
PGLS regression of anterior MA values as a function of dietary preferences and UBF groups reveals a statistically significant correlation (R 2 = 0.1692, P = 0.001; Table 2) that is stronger than the relationship between beak shape and those measures of dietary ecology. Higher values of MA are consistently associated with cracking/biting taxa, and those whose diets rely heavily on plant matter, with large proportions of items such as fruits and drupes, seeds, bulbs, shoots, grass, or leaves (Fig. 6). Phylogenetic MANOVA revealed no pairwise differences between any of the groups based on MA values (Table 3). We found a strong significant interaction between dietary preferences and UBF groups (R 2 = 0.26376, P = 0.001) revealing that the relationship between diet and MA varies depending on the feeding behavior (Table 2; Fig. S6).

BODY MASS AND FEEDING ECOLOGY
PGLS regression of log-BM as a function of dietary preferences and UBF groups reveals a stronger correlation of body size with feeding ecology than that of both beak shape and MA with feeding ecology, with dietary variations explaining as much as 25% of log-BM variation (Table 2). Visual inspection of the regression scores of log-BM associated with the first axis of diet variation (PCo1) reveals that taxa with large amounts of invertebrates in their diet tend to be smaller, whereas some dietary groups such as scavengers tend to be associated with bigger sizes (Fig. 6).
UBF groups are only weakly associated with log-BM and none of the UBF groups are statistically different to any other in log-BM (Table S3), although significant diet/UBF interactions reveal that different behavioral groups exhibit different body size to diet relationships (Table 2; Fig. S6).

Discussion
Our analyses aimed to quantitatively test the common wisdom that feeding adaptation is one of the main drivers of beak morphological diversification in modern birds. Our results suggest that adaptation to dietary composition is not as fine-tuned as generally perceived, and there is not a close to one-to-one mapping of beak shape on feeding ecology. At a broad macroevolutionary scale, we found a more complex but weak overall covariation between beak shape and diet, with other factors such as biting MA and body size being stronger covariates for feeding autecology. Similar beak shapes are associated with the increased consumption of different food items (i.e., a one-to-many relationship between shape and ecology) and the relationship between beak shape and dietary preferences is different within different UBF groups, likely owing to the ecological heterogeneity of feeding behavior groups (i.e., many-to-one ecology to behavior relationships). For instance, probing birds in our sample are composed primarily of two very ecologically different groups: longirostrine waders (e.g., Numenius, Gallinago, Limosa) and the kiwi (Apteryx), and anseriforms (e.g., Aythya, Anas, Cygnus), which both use the UBF as a probing tool in (mostly) soft substrates (Figs. S1 & S6).
Our results suggest that the beak is generally used as a versatile, tweezer-like clamp. Mechanical preprocessing of food (i.e., tearing and cracking/biting feeding behaviors) is generally associated with deep and curved beaks, which are able to accommodate comparatively higher stresses than thinner, straighter beaks (Soons et al. 2010;Soons et al. 2015). Similarly, beaks well suited for sensing and probing in fluid or soft soils tend to be long and thin (Barbosa and Moreno 1999). Although such shapes represent the ends of a clear ecomorphological spectrum, it is difficult to predict where a given species should fall upon it, as species well suited for performing a certain feeding behavior may not actually use their beaks in the way we would expect given their morphology (e.g., the kakapo, Strigops, has a typically parrot-like beak well suited for cracking/biting, yet chooses to feed on soft leafy vegetation rather than fruits or seeds). Most of the species studied fell between these extremes in ecomorphology, using the beak for grabbing/gleaning or pecking/grazing and exhibiting a broad range of beak morphologies therein (i.e., many-to-one mapping of shape and behavior). Furthermore, the majority of bird taxa show values of anterior MA congruent with fast gapes and low bite force transmission, and many of these belong to the grabbing/gleaning behavioral group, which occupies virtually all of beak shape and functional space We found a significant relationship between beak shape and MA: increased values of anterior MA are strongly correlated with increased beak depth/length ratio, driven, in part, by shortening of the beak, and suggesting that enhanced biting force transmission requires a deeper beak to accommodate higher stresses and avoid fracture (Soons et al. 2010, Soons et al. 2015. However, this relationship differs between taxa, and thus indicates a many-toone relationship between shape and this functional trait. Raptorial birds are interesting, as they have much lower anterior MA values than predicted by the general regression. Initially this may be surprising, given the predatory nature of raptors, yet this result is congruent with previous research showing that Strigiformes and Accipitriformes rely heavily on talon adaptations to kill their prey (Sustaita 2008;Sustaita and Hertel 2010;Del Hoyo et al. 2017;Madan et al. 2017). Deep beak morphologies are, however, associated with enhanced biting MA in the two taxa representing falconin falconiformes (Falconinae, Falconidae; Falco and Herpetotheres). Falcons dispatch prey with their beaks rather than their talons (Sustaita 2008;Sustaita and Hertel 2010;Del Hoyo et al. 2017), which may explain why both falconid taxa differ from the other raptors and instead follow the general regression trend for all avians.
The evolution of faster gapes and comparatively weaker bite force advantage happen primarily within the Charadriiformes (i.e., Scolopacidae). Unique modes of cranial kinesis, such as distal and double rhynchokinesis (i.e., avian cranial kinesis characterized by additional bending areas in the tip of the beak, and in both the tip and the base of the beak, respectively; Zusi 1993; Estrella et al. 2007), appear in this clade of mainly probing taxa, and could further enhance gape speed. In contrast, comparatively slower gapes and enhanced biting force transmission evolve less frequently. Parrots (Psittaciformes) are the most notable and extreme example, especially when we consider that their MA values here may be underestimated, thanks to novel adductor muscles and skeletal adaptations which may enhance lever efficiency in some parrots (Zusi 1993;Tokita et al. 2007). Our results suggest that dietary transitions toward increased herbivory are correlated with evolutionary changes toward higher anterior MA, implying that herbivory imposes higher performance demands on the beak. This observation is congruent with previous ecomorphological studies on waterfowl (Olsen 2017) The transfer of grasping and manipulation behaviors from the forelimbs to the beak in bird evolution has necessitated that bird beaks be highly versatile, used in virtually every aspect of their biology, not just feeding and foraging (Bhullar et al. 2016). The complex evolutionary scenario demonstrated by our results suggests that diverse and multidirectional selective pressures were involved in beak morphological diversification, reflective of functional and behavioral multitasking. In this evolutionary context, a fast, generic grabbing tool could most easily fit the required compromise of functional versatility (i.e., trade-off between varied beak functions), explaining the prevalence of thin and straight beak shapes and optimization for low-force transmission highspeed gapes in our sample. More nuanced relationships between feeding adaptation and beak shape may be operating, with variable strength, within lower taxonomic levels, to accommodate different macroevolutionary regimes and trade-offs. For example, while a strong association between feeding ecology and beak shape characterizes the diversification patterns within waterfowl (Olsen 2017), skull CS, not diet, is a major driver of beak shape in diurnal raptors (Bright et al. 2016). Nevertheless, our data support the idea that beak shape and MA reflect the mechanical demands of specific feeding and foraging strategies (Bowman 1961;Schwenk 2000). This relationship may be best envisioned as a threshold rather than a one-to-one connection, with certain shapes and mechanical properties critically needed to perform certain functions and feeding behaviors (e.g., to avoid fracture). In agreement with these views, some species of Darwin's finches show dietary habits and feeding strategies that are more flexible than previously thought; their specialized beak phenotypes (e.g., cracking/biting) are still efficient in processing many other dietary resources, which might lead to the evolutionary retention of these phenotypes (i.e., Liem's paradox; Tebbich et al. 2004).
In conclusion, our results imply that the relationship between beak shape and feeding ecology at a broad macroevolutionary scale may be more complex than usually assumed. This is particularly important in fossil taxa, where trophic hypotheses are rarely testable (e.g., fossilized gut contents). In light of these results, it is important to evaluate the strength of the relationships between form, functional traits, and feeding behavior within a taxonomic context, before drawing trophic assumptions based solely on beak morphology. In doing so, we will open pathways for a more detailed understanding of the role of trophic adaptation in shaping avian diversity.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article.     Figure 5. The relationship between beak shape and dietary preferences. Figure 6. The relationship between beak shape, anterior MA, and body mass with dietary preferences. Figure 7. The relationship between beak shape and dietary preferences. Table 2. Number equivalences for clades in the phylogeny in Figure 1. Table 3. Phylogenetic MANOVAs. P values for the pairwise differences in effect sizes (Z scores) between groups. Table 4. Diet-dependent allometries. Table 5. Summary of the full PGLS linear model for Procrustes coordinates (beak shape) as a function of dietary preferences, log-CS and UBF categories (including main effects of both independent variables and their interaction). Cells in bold indicate statistical significance (P < 0.05). Effect sizes (Z) are computed as standard deviates of the F values' randomized sampling distributions. P values are calculated for the F values' randomized sampling distributions. Table 6. Evaluation of evolutionary models for the residuals from the PGLS linear regressions used in this study.