Evolutionary patterns and processes in the radiation of phyllostomid bats
© Monteiro and Nogueira; licensee BioMed Central Ltd. 2011
Received: 1 November 2010
Accepted: 23 May 2011
Published: 23 May 2011
Skip to main content
© Monteiro and Nogueira; licensee BioMed Central Ltd. 2011
Received: 1 November 2010
Accepted: 23 May 2011
Published: 23 May 2011
The phyllostomid bats present the most extensive ecological and phenotypic radiation known among mammal families. This group is an important model system for studies of cranial ecomorphology and functional optimisation because of the constraints imposed by the requirements of flight. A number of studies supporting phyllostomid adaptation have focused on qualitative descriptions or correlating functional variables and diet, but explicit tests of possible evolutionary mechanisms and scenarios for phenotypic diversification have not been performed. We used a combination of morphometric and comparative methods to test hypotheses regarding the evolutionary processes behind the diversification of phenotype (mandible shape and size) and diet during the phyllostomid radiation.
The different phyllostomid lineages radiate in mandible shape space, with each feeding specialisation evolving towards different axes. Size and shape evolve quite independently, as the main directions of shape variation are associated with mandible elongation (nectarivores) or the relative size of tooth rows and mandibular processes (sanguivores and frugivores), which are not associated with size changes in the mandible. The early period of phyllostomid diversification is marked by a burst of shape, size, and diet disparity (before 20 Mya), larger than expected by neutral evolution models, settling later to a period of relative phenotypic and ecological stasis. The best fitting evolutionary model for both mandible shape and size divergence was an Ornstein-Uhlenbeck process with five adaptive peaks (insectivory, carnivory, sanguivory, nectarivory and frugivory).
The radiation of phyllostomid bats presented adaptive and non-adaptive components nested together through the time frame of the family's evolution. The first 10 My of the radiation were marked by strong phenotypic and ecological divergence among ancestors of modern lineages, whereas the remaining 20 My were marked by stasis around a number of probable adaptive peaks. A considerable amount of cladogenesis and speciation in this period is likely to be the result of non-adaptive allopatric divergence or adaptations to peaks within major dietary categories.
The Phyllostomidae (leaf-nosed bats) is the dominant family of bats in Central and South America. This family has undergone an adaptive radiation unparalleled among other mammals in terms of ecological and morphological diversity . Starting from an insectivore ancestor in the late Eocene [2–4], the 53 extant genera in this family have diversified into specialized forms for insectivory, carnivory, frugivory, granivory, nectarivory, and sanguivory (although many species have mixed diets) [5–7]. This ecological diversity and specialisation originates from an intricate partitioning of resources [8, 9], and is probably responsible for the high local species richness (ranging between 31-49 syntopic species [10, 11]) observed for leaf-nosed bats. The evolution of specialised diets created functional demands, apparently determining changes in cranial and mandibular shape [1, 12]. The magnitude of skull form (shape + size) variation among phyllostomid lineages is large and has been assessed by correlational studies both qualitatively and quantitatively, using measurement ratios of functional relevance or traditional distance measurements [1, 12–16]. One important aspect of phenotypic variation in the family is the snout elongation associated with nectarivory. This elongation is thought to be responsible for a trade-off between two functional demands: tongue support and bite force [8, 9, 17]. Bats with longer snouts might have longer operational tongue lengths , but are known to have weaker bites than bats with short snouts, what might restrict the dietary range accessible to them  and lead to seasonal migrations . Nogueira et al. (2009)  have shown that apart from the general elongation, other skull and mandible shape changes are associated with (size-independent) bite force, such as the relative size of mandibular processes (coronoid, angular), zygomatic arch position and robustness and the relative position and sizes of tooth rows. The feeding behaviour is also relevant to the understanding of the biomechanics of feeding and resource partitioning in bats [14, 20], but is less studied than morphology.
The ecomorphological diversification of phyllostomids has long been considered the result of an adaptive radiation, but no specific tests of the responsible mechanisms have been performed, apart from the correlational studies mentioned before (most of them not using comparative methods). An analysis of diversification rates indicated a significant shift at the base of the phyllostomid tree , that could be associated with an adaptive diversification, but increased speciation is not an unequivocal evidence of adaptive radiation [22, 23]. Monteiro and Nogueira (2010)  provide indirect evidence of adaptive evolution, assessing the integration patterns in the phyllostomid mandible during evolutionary shape changes. The interspecific integration patterns (correlated shape changes among mandibular components) were independent from pooled within-species integration patterns (which mirrored mammalian developmental genetics patterns), as expected during adaptive evolution on an adaptive landscape [25, 26]. Long term selection on species means (macroevolutionary changes) is expected to be independent from the structure of genetic correlations among morphological variables, depending only on selection gradients [25, 26] associated with specific adaptive peaks. In this context, we can move beyond the usual phenotype-ecology correlation approach, assessing the likelihood of different evolutionary scenarios through recent model-based approaches for comparative analyses [27, 28]. It is possible to select between alternative adaptive models with different postulated selective agents or adaptive peaks . It is also possible to estimate optimal phenotypic values for each adaptive peak to use as a basis for testing hypotheses regarding the adaptive evolution of lineages , and to cast light on the evolutionary processes responsible for species diversity .
We focus our study on the variation of mandible form (shape and size), which is a model for the evolution of complex morphological structures . The mandible can be used as a proxy for the facial skull due to developmental integration of jaw (mandible-maxilla) parts . Unlike other parts of the skull that harbour different functions (protecting the brain and sensory organs), the mandible's main functional demands are related to feeding, and dietary changes are expected to be the main selective agent for this structure. We combine ecomorphological correlations, patterns of disparity through time and a model-based comparative analysis to test evolutionary mechanisms and scenarios during the radiation of phyllostomid bats.
Multivariate PGLS regression results for the shape PCs on size and diet.
R2 = 0.570
R2 = 0.807
R2 = 0.444
R2 = 0.440
R2 = 0.159
The disparity through time (DTT) plots provide further evidence of a two-mode evolutionary history for shape and diet (Figures 5C and 5E). These plots show average sums of morphometric and dietary distances (relative to the total sum or total disparity) for lineages of a given (relative) age range through the phylogenetic tree. The early evolution of phyllostomid lineages is marked by larger disparity (blue solid line) than expected by 1000 neutral evolution simulations (median and 95% confidence intervals of simulations depicted as solid and dashed red lines, respectively). During this initial period, there is a high average disparity within existing lineages. Both shape and diet disparities peak above the expected neutral disparity around 25 Mya. There is also a turning point after this early diversification (around 20 Mya), when observed disparity within lineages is smaller than expected under neutral evolution, suggesting that stabilising selection was the main evolutionary force constraining diversification within lineages after the main ecological divergence had taken place. Size disparity through time follows a similar pattern initially, but there is a second peak above expected median disparity between 16-7 Mya. This second peak is unique for size disparity and takes place after the main dietary categories were established. The relevance of this pattern is not so clear because the confidence limits for univariate characters are considerably large and much of the size disparity fluctuation within lineages falls inside the confidence intervals. The five panels (Figure 5A, B, C, D, E) are complementary in the sense that they are indicating a duality of evolutionary processes at different times during the phyllostomid radiation, and that shape and size did not have a correlated pattern of disparity through time. The steep slopes observed in the right hand extremity of Figures 5A and 5B are consistent with the larger than expected shape disparity in the deeper parts of the tree (between 30 and 20 Mya - Figure 5C), and both are predicted under directional selection. The near zero slopes in the left hand extremity of Figures 5A and 5B are consistent with the lower than expected disparity within younger clades (Figures 5C and 5D), and both patterns are likely under stabilising selection.
Performance of alternative models for the evolution of mandible shape using principal components.
Models with 3 PCs
Models with 5 PCs
Performance of alternative models for the evolution of skull length.
The combination of ecomorphological correlations, disparity analysis and model-based comparative techniques provided strong evidence for the evolutionary mechanisms responsible for the phyllostomid radiation. Although there are a few instances of convergence in phyllostomids [3, 4], the dietary divergence presents a strong phylogenetic structure [24, 36], where most specialisations occurred only once and many of the main lineages (recognised as subfamilies in most taxonomic studies ) are relatively homogeneous morphologically and ecologically. Convergence is an important part of the study of adaptation , and its absence makes it harder to infer the role of natural selection in biological diversification. One has to be particularly careful with the comparative methods and evolutionary model assumptions, for using the wrong models or methods are likely to lead to erroneous conclusions [28, 38, 39]. For example, in the present study, forcing a Brownian model of evolution into the phylogenetic regression model would lead to non-significant associations among shape and diet principal components (results not shown). Flexible approaches allowing for different evolutionary models [27, 28, 40] provide more sophisticated and interpretable results , as well as more interesting questions to be explored. The model-based approach used in this study shed light on the evolutionary mechanisms responsible for the evolutionary radiation of phyllostomid bats and makes it possible to generate predictions of optimal shapes (or ratios of functional relevance) for dietary specialisations that can be used in the future to test hypotheses of biomechanical optimisation .
The major axis of shape variation among phyllostomid species (first shape principal component) ordinates species with contrasting degrees of rostral elongation. This shape change is regarded by Freeman (2000)  as one of the "cheap tricks" of mammal diversification, and studies on other groups of mammals such as domestic dogs  suggest that this pattern of shape variation can arise in a microevolutionary scale (~10000 years) as a response to selective processes. Rostral elongation in phyllostomids is associated with a trade-off between support for the elongated tongue and bite force (species with shorter rostra usually have stronger bites) [17, 43, 44]. This reduced bite force is a result not only of the increased out-lever , but also of the relative decrease of muscle insertion areas and robustness of elongated skulls and mandibles . The bite force constraint caused by morphological specialisation is known to limit the dietary scope [8, 9] and might have important ecological consequences, particularly regarding patterns of resource use and foraging strategies within guilds .
Mandible and palate length are generally good predictors of operational tongue length in most nectarivores  (but see [45, 46] for an alternative adaptation in Anoura fistulata). Longer tongues are considered an adaptation for nectar extraction, not only because these will reach flowers with longer corollas , but because they allow the nectar specialists to explore a larger variety of plant species more efficiently (avoiding seasonal changes to insects and fruit when local nectar abundance decreases) [19, 47], and possibly to maintain a longer distance from the flower during feeding (reducing predation risk) . It is currently established that nectarivory has evolved twice independently within phyllostomids [3, 4], and despite a superficial phenotypic convergence (both lineages do present the mandibular elongation), there are many anatomical differences, particularly in tongue morphology . All this morphological variation suggests that, in this system, there are many ways to achieve the same function (many-to-one relationships ) and that there is more variation in diet and foraging strategies among nectarivores than implied by the common use of such category. Brachyphylla is a singular nectarivore in the sense that it presents a short mandible (the blue sphere with negative score on shape PC1 - Figure 3A, see also Additional file 2) and is phenotypically similar to Phylloderma and Sturnira, being ordinated between frugivore and insectivore species. The similarity between Brachyphylla and the frugivore Stenodermatine lineage has been recognised for other morphological characters in the dentition and cranial shape [1, 49] and could be the result of dietary changes to exploit open niches during island colonisation, as suggested by Griffiths (1985) .
The second principal component of interspecific shape variation depicts a common pattern of shape change for diets with low mastication (shared features between nectarivores and sanguivores) versus high mastication (shared features between animalivores and frugivores). Both sanguivore and nectarivore specialists do present lower bite forces than expected for their sizes . Although nectarivores and sanguivores are in extreme opposites regarding the mandibular elongation pattern of shape PC1, they do present similarities in the relative reduction of mandibular processes (particularly the coronoid) and molar series (more pronounced in sanguivores). These shared features are consistent with developmental consequences of selection for morphological specialisations . This is clear for nectarivores, where selection for a longer, narrower skull has limited the relative size of muscles and their areas of attachment . The explanation for the pattern of shape changes in sanguivores is not as direct as that for nectarivores, as the evolution of different mandibular components towards this adaptive peak has been shown to disagree with expected developmental patterns , where tooth bearing components and mandibular processes form almost independent modules. This is because the two tooth bearing components seem to change their shapes independently: a relative increase in the anterior alveolar (incisives and canines), combined with a relative decrease in the posterior alveolar (molar series). This emphasis in anterior dentition is also observed in tooth development as well, and was hypothesised to be the result of selection to sanguivory . Therefore, even though a reduction of bite force and muscle masses is observed in both sanguivores and nectarivores [8, 51], the developmental consequences of these changes are more localised (in the posterior mandible) in sanguivores than in nectarivores. This is probably due to stronger stabilising selection in the anterior region of the mandible of sanguivores . The evolution of sanguivory has been discussed in the literature, but the current limitations of available data make it difficult to get past the delineation of general scenarios and hypotheses  to actually testing functional predictions about agents of selection and phenotypic variation. Sanguivore skulls did not evolve according to expected developmental integration patterns , and do not meet traditional biomechanical predictions based on models of mastication . The functional demands on the skull associated with sanguivory are not known in detail, but should include not only biting off flesh and lapping blood, but also sensorial tasks in finding and approaching potential bloodmeal sources. Sanguivores are unique among phyllostomids in that there are no known intermediate forms (both in diet and phenotype), whereas for the species that specialised in plant items (frugivores and nectarivores) there are plenty of intermediate species with mixed diets. This pattern suggests that the divergence of frugivores and nectarivores might have occurred along ridges and local maxima on the adaptive landscape, but the evolution of sanguivores has probably occurred via a leap over a fitness valley in a short period of time (between 31 and 21 Mya ). Because the transition period was very short, useful fossils documenting it would need to be obtained from this specific window on the late Oligocene. However, all known fossils of vampire bats are Plio-Pleistocenic Desmodus and Diphylla , which are the same or very similar to extant vampire species.
The first and second major axes of variation in the mandibular shape space of phyllostomids mostly contrast groups that present high and low masticatory demands. The phenotypic differences between dietary specialists with high masticatory demands (frugivores and animalivores) are depicted only on the third and fourth principal components of shape. These differences are associated with relative positions and sizes of mandibular processes (lower condyle, smaller angular in frugivores), tooth row lengths (longer molar rows in animalivores), and the anterior region of the mandible (mandible shorter and deeper, forming a "chin" in frugivores). Aside from large beetles (which seem to be harder than any other food item processed by bats ), insects and fruit do not present significantly different hardnesses and there is a lot of variation within these general groups, mostly due to size differences of the dietary items (larger items are harder [9, 15, 53]). Animalivores and frugivores do not present noticeable differences in size-independent bite forces [9, 17], and there is a considerable amount of behavioural plasticity that allows them to modulate the bite force according to food hardness . The extensive shape changes in the skull, mandible and dentition [1, 13] observed in the evolutionary transition towards frugivory did not seem to be associated with changes in bite force . Different from sanguivores and nectarivores, frugivores needed to maintain the bite strength for mastication, while adjusting for their specific functional needs (which are varied, according to the main type of fruit consumed [6, 13, 15, 54]). The morphological changes observed in frugivores seem to sacrifice precision of dental occlusion (a feature apparently more important in animalivores) to favour skull and mandible robustness , while maintaining a small body size , particularly in the shorter, deeper mandibles observed in the small clade of specialist frugivores (the Stenodermatina).
Allometric variation is a controversial theme in phyllostomid evolution and ecology. There is no doubt that body size is the main factor causing differences of absolute bite force  and will have a consequence on the dietary scope of particular species [9, 12]. On the other hand, small frugivores, such as Centurio are able to produce higher bite forces than expected for their sizes  due to cranial shape changes. Carnivory is traditionally associated with an increase in body size and is considered an allometric extrapolation of insectivory . One small group of strict frugivores (Stenodermatina, see above) presents small average size, but it is not clear whether their size evolution was associated with diet or with the postulated insular origin of the clade  (see  for a review of insular size patterns in vertebrates). No obvious size trends can be associated with nectarivores or sanguivores. The ecological evidence for size as a structuring influence in phyllostomid communities is equivocal , even when considering more detailed dietary compositions. Our results do show that the evolutionary trajectories followed by diverging phyllostomid lineages in mandible shape space were too complex to be explained by simple allometric changes, and the shape changes associated with size differences, as measured by skull length (towards extreme frugivores and carnivores) were small, when compared to the drastic shape changes associated with nectarivory and sanguivory, which did not involve significant size changes. Shape and size divergence patterns seem to be decoupled in phyllostomid evolution (see also discussion below), but size is nevertheless a relevant dimension in the phyllostomid adaptive landscape. A measure of body size would probably be more suited for such ecological associations than cranial length.
The patterns of morphological and ecological divergence over time were informative regarding possible evolutionary mechanisms in the radiation of phyllostomid bats. When considering the parallel evolution of mandibular shape and diet, one sees clearly a two step radiation, marked by strong directional selection and divergence in the first phase (30-20 Mya) followed by stabilising selection and stasis in the second phase, after the main phyllostomid lineages (and diet specialities) appeared. The scatterplots comparing morphometric distances with genetic distances and time since divergence can be qualitatively compared to patterns from Polly's simulations  of evolutionary mechanisms. Different mechanisms are expected to leave distinctive imprints in the patterns of morphological divergence over time. In real data, noisy patterns are expected from the juxtaposition or superimposition of different evolutionary processes and it is harder to extricate the processes from the patterns. The scatterplots of genetic distances (percentage of sequence divergence), time since divergence and morphometric distances in Figures 5A, B show a two-phase pattern, which is more obvious in Figure 5A (genetic distances) than in Figure 5B (time since divergence) because of species that are phenotypically and genetically similar (insectivores) but have diverged in beginning of the radiation. The disparity through time plots for mandible shape and diet suggest also a two-phased radiation, marked initially by larger disparity than expected by neutral evolution simulations and changing lo lower disparity than expected by neutral evolution after around 20 Mya (the transition between Oligocene and Miocene - see discussion of the time frame for the phyllostomid diversification in [3, 4]). At this transition point, the main lineages leading to dietary specialisations were established and most shape and diet disparity is observed among these lineages. Size disparity presents a different pattern from shape and diet, with at least two distinctive peaks at different times in phyllostomid evolution. The first peak occurring between 30 and 20 Mya, and it was probably caused by the large carnivores sharing the same branch with other groups until nearly 20 Mya (these are very different in skull length but not so different in mandible shape). The second peak happens between 16-7 Mya and is likely to be caused by disparity within Stenodermatinae  (The clade in Figure 1 joining species from Sturnira to Ametrida). The Stenodermatina (the clade joining species from Ariteus to Ametrida) have the shorter skulls among the phyllostomids (all of them smaller than 14 mm average - see Additional file 1), the within lineage disparity only decreases after they separate from the other Stenodermatinae (which have skull lengths comparable to the remaining lineages) around 10 Mya. Because the confidence limits for the neutral simulations are very wide, the fluctuation of size disparity through time never exceeds the expected disparity. On the other hand, it is not clear whether these confidence limits are too conservative or realistic. Therefore we have interpreted and discussed the size disparity results. This pattern of size disparity is clearly not associated with the mandibular shape disparity and dietary disparity. From the other analyses performed in this study, it is clear that evolutionary allometry is a significant but minor influence in mandibular shape diversification. On the other hand, the lack of correspondence between diet disparity and size disparity might also be explained by the lack of detail of the diet data within major categories.
The radiation of phyllostomid bats has commonly been referred to as an adaptive radiation [1, 16]. However, a complete test of this hypothesis has not been thoroughly performed (although different criteria have been examined separately). According to Schluter (2000) , there are four criteria to determine if a radiation was adaptive: common ancestry, phenotype-environment correlation, trait utility and rapid speciation. The monophyly of phyllostomids (common ancestry criterion) was long established by a number of studies using molecular and morphological data sets [3–5, 21, 33]. The rapid speciation criterion has been assessed by Jones et al  whose findings indicate a statistically significant shift in diversification rates in phyllostomids (a hypothesis that would not be testable with the incomplete tree used in the present study). These authors found two significant shifts in phyllostomids: one just after the vampire bats (Desmodontinae) branch out and one within the genus Artibeus. The first diversification rate shift observed by Jones et al  within phyllostomids coincides with the morphological and ecological disparity peak observed by us at the base of the phylogenetic tree. The phenotype-environment correlation and trait utility criteria have been discussed before in the literature and in the present study. The radiation of phyllostomids in mandibular shape space occurred along axes leading to postulated adaptive peaks determined by the relative importance of dietary items. The shape changes along these main axes of variation have a clear functional relevance for the acquisition and processing of food, as discussed above and in the literature [1, 16, 17, 24, 58].
The model-based analysis presented in this study provided quantitative evidence of evolutionary mechanisms and scenarios. The models describing phenotypic evolution on an adaptive landscape with five different peaks fit the data better than the neutral evolution model (Brownian motion) and simpler models with less peaks. The favoured scenario combining all analyses would be one early burst of phenotypic and ecological diversification caused by directional selection towards five different adaptive peaks (insectivory, carnivory, sanguivory, nectarivory and frugivory). In our interpretations, we gave more weight to the Akaike information criteria corrected for sample size because in simulations, it performs better than the Schwarz information criteria when reality is assumed to be infinitely dimensional and the true model is not in the candidate set . The results for models with different numbers of variables (shape PCs) show how sensitive the methods are to the number of parameters being estimated. In the models with the first three shape PCs, the number of parameters being estimated (DOF in Table 2) is never larger than the number of species (49), whereas with the first five shape PCs, OU.4 and OU.5 (the best fitting models otherwise) do exceed the number of species and are severely penalised by the more conservative criterion (SIC, which assumes that the true model is in the candidate set and is low dimensional ). As a consequence, Brownian motion is suggested as the best fitting model because it requires estimation of the smallest number of parameters. These results indicate that we were working at the limit of model complexity given our sample size, and more complex models would require considerably larger data sets (in terms of number of species) to avoid problems in estimation and model comparisons.
When using the first three shape PCs, model OU.4 was favoured, whereas model OU.5 was favoured by the models with five PCs. The separation of the carnivore from the insectivore peak is only accomplished when the fifth PC is added to the models. One might expect that five groups would span a space with four dimensions, however, because principal components do not maximise among group differences, there will be processes other than diet-related adaptive changes contributing to the observed interspecific variation patterns. These could be associated with neutral evolution or changes related to factors not measured in this study. In any case, it is a noticeable point that Horn's parallel analysis used to choose the number of principal components to be considered in the comparative analyses, independently indicated the number of components needed to discern all possible adaptive peaks in the most complex model. Parallel analysis is often shown by statistical simulations to be the best performing method to choose principal components [60, 61].
The early burst of phyllostomid divergence was followed by consistent stabilising selection keeping mandible shape relatively constant around postulated dietary optima. The apparent phenotypic and ecological stasis that persists from the early Miocene (20 Mya) to the present was not followed by a lack of speciation, and the lineages of specialist nectarivores and frugivores seem to be particularly speciose . This pattern suggests either a non-adaptive radiation  or agents of selection not specifically examined in the present study. Allopatry and biogeographic distribution patterns account for a considerable proportion of speciation within diverse phyllostomid genera where stabilising selection seems to constrain phenotypic and ecological variation [62–64]. In fact, niche conservatism and stabilising selection are expected to play a significant role in allopatric speciation processes , what seems to fit well with the environmental changes during the Tertiary [4, 66], when the global temperature decrease, tectonic processes (Andes uplift) and sea level fluctuations created large expanses of dry open areas in South America and are linked with the diversification of many other mammal clades [67–69]. On the other hand, further selective episodes cannot be discarded, as there might be smaller adaptive peaks within the main dietary groups and the temporal and spatial variation in resource availability (particularly fruits and flowers) might generate randomly fluctuating selection that is hardly discernible from a Brownian motion [25, 34]. Frugivore phyllostomid species are known to largely occur in sympatry [10, 15], possibly due to further specialisation, such as the dichotomy between ground-story frugivores and fig-feeders , or the recently discovered granivory  associated with functional and morphological specialisations . Nectarivore guilds can also be diverse  due to nomadic behaviour, seasonal and spatial changes in resource use [18, 47] and different degrees of specialisation and dietary item mixtures . The radiation of phyllostomid bats is actually a number of radiation episodes nested within each other, caused by a mixture of adaptive and non-adaptive evolutionary mechanisms. This is expected as real data is more likely to fall along a continuum caused by a mixture of evolutionary processes, rather than fit yes/no definitions for the adaptive radiation metaphor .
The radiation of phyllostomid bats is a unique example among mammals in terms of ecological and phenotypic diversity, and should be considered an important model system for studies of the evolution of functional optimisation due to the ecological diversity under the constraints imposed by powered flight . Further research should provide a deeper understanding of the myriad of evolutionary mechanisms at work in this lineage. The foundations of phyllostomid ecomorphology, based on rather limited morphological, ecological, phylogenetic and functional data painted an interesting picture and provided important references and new questions to be addressed. Recent contributions comprise functional studies [41, 73] that go beyond traditional biomechanics to examine the functional consequences of shape change in terms of energy efficiency and structural resistance, comparative analyses of detailed morphological, functional, and ecological data [17, 24], comparative analyses combining field data on function, behaviour and ecology [7–9, 20, 43, 44, 74] and phylogenetic diversification patterns using complete species-level trees . Whereas a lot of attention has been paid to the biomechanics of mastication and bite force, the functional demands associated with nectar and blood-feeding are still underrepresented in the literature. Of particular importance would be the identification of functionally relevant measurements (with fitness consequences) with regard to different dietary compositions to be associated with phenotypic data. When comparing the abundance of phyllostomid species in scientific collections, Freeman (2000)  has pointed out that species with specialised phenotypes and diets are rarer than the species with intermediate phenotypes and mixed diets. Are dietary specialisations equivalent to local optima on the adaptive landscape with lower fitness than optima for mixed diets? This is an important question for future research with possible implications for the dynamics of assemblages  and biodiversity conservation. Other promising areas are the ecological influences on cranial phenotypes outside dietary differences, such as echolocation and roost building [1, 75]. The ever increasing literature is producing a massive database of morphological, ecological, functional and phylogenetic data that will be instrumental to elucidate the questions in future studies of phyllostomid ecomorphology.
The patterns of phenotypic and ecological divergence are consistent with an evolutionary scenario with at least five main adaptive peaks during the early radiation of phyllostomid bats. Starting from an insectivore ancestor around 30 Mya, different lineages specialised into carnivores, sanguivores, frugivores and nectarivores. Some species are strict diet specialists with extreme morphologies, but a considerable number present intermediate diets and phenotypes, possibly as a result of adaptive ridges or multiple local optima (many phenotypes with equivalent or similar fitness) and geographical/temporal variation in resource availability. Nevertheless, some trade-offs are clear, such as the mandibular elongation observed in specialised nectarivores, which supports a longer tongue but decreases bite force changing the scope of usable resources (increased variety of flowers versus harder fruit and insects). Size and shape have evolved almost independently in this family, where the major axes of shape change (towards nectarivory and sanguivory) in phyllostomid mandibular morphospace were not correlated with size differences. On the other hand, the most impressive size differences (towards carnivory) are associated with minor shape changes. After the early burst of ecomorphological divergence, when directional selection was possibly a dominating evolutionary mechanism, some of the lineages went into a long period of stasis or non-adaptive radiations, but further selective episodes cannot be ruled out, particularly for size variation among frugivores, although the role of diet as selection agent is not entirely clear in this case. The radiation of phyllostomid bats was marked by a complex mixture of adaptive and non-adaptive mechanisms through a period of extreme environmental changes when new ecological niches were probably emerging and disappearing quickly. This led the leaf-nosed bats to present unparalleled morphological and ecological variation among mammals, which, together with the vast amount of information available in the literature, makes phyllostomids one of the most important model systems for the study of morphological, functional and ecological diversification.
A total of 443 specimens representative of 49 genera of phyllostomid bats were analysed (Figure 1). Sample sizes and species names are detailed in Additional file 1. All specimens were identified as adults based on the ossification of the basisphenoid region. Mandible images were obtained with a digital camera Nikon Coolpix 8700. All specimens were photographed under the same plane in relation to the camera lens (reference points in the specimens were used to align the mandible with the focal plane). To capture as much shape variation as possible, 11 landmarks, and 25 semilandmarks were used (Figure 2). The reference points were digitalised in the TpsDig software . The landmark configurations were superimposed according to a least squares criterion (Procrustes superimposition) [77–79], first calculating a mean for each species and then superimposing all configurations to calculate a grand average. The semilandmarks were slid to minimize the least squares criterion  in TPSRelw . Even though slid semilandmarks do not hold biological correspondence (just like type II or III landmarks), there is a gain in the information on shape variation (as ascertained from a comparison of superimposed warped images in the program TpsSuper  - results not shown), and the correlation of Procrustes distance matrices obtained from data sets with and without semilandmarks was 0.956. After Procrustes superimposition, the aligned coordinates present more dimensions than the actual shape space (aside from the degrees of freedom lost in superimposition, a maximum of 48 dimensions is needed to describe 49 species means). To reduce dimensionality, the aligned coordinates were transformed to shape principal components . Principal component analyses of aligned coordinates and the visualisation of shape changes associated with each principal component (shape PC) axis were performed in MorphoJ (importing the aligned landmarks and semilandmarks from TPSRelw) . We have used Horn's parallel analysis [60, 61] to choose the number of principal components to be retained. Parallel analysis is a method of consensus in the literature , and consists of generating a large number of random samples with the same number of variables (but uncorrelated) and specimens as the original data, and comparing the observed eigenvalues to the distribution of random eigenvalues. The observed PCs that presented eigenvalues larger than a percentile (we used 95%) of the random eigenvalue distribution are retained. We used the paran package  in the R environment  to perform the parallel analysis. Size was measured both as the centroid size (square root of summed squared distances from each landmark to the centroid of points) from mandible landmark configurations , and as Condylobasal Length (CL - from the anteriormost point on the premaxilla to the condyle). These two measures of size were highly correlated (r = 0.983), and all results were the same for both variables. Because cranial length is a more intuitive measure, and to facilitate further comparisons, we chose to report results using CL instead of mandible centroid size.
Dietary information was obtained from the literature (modified and updated from ). We organized the quantitative and qualitative information for the 49 bat species into five categories comprising the main food items identified in phyllostomid diets (insectivory [arthropods], carnivory [small vertebrates], frugivory, nectarivory [pollen and nectar], and sanguivory). With the exception of sanguivore bats, all phyllostomid species present mixed diets, with varying levels of specialization and predominance of items from a given category. To compare the relative importance of structurally different items, such as nectar and arthropods in the diet of a single species is a challenging task. We defined rank categories based on relative usage of food items (0 - absent, 1 - complementary, 2 - predominant, 3 - strict), and a rank was attributed to each dietary item for each species (see Additional file 1). The distribution of diet variables along the phylogenetic tree used  is depicted in Figure 1. The diet variables were used later as independent variables in a multivariate regression model. Preliminary analyses indicated strong multicollinearity among the diet variables, mostly caused by a positive correlation between insectivory and carnivory (sometimes lumped as animalivory), and a negative correlation between frugivory and animalivory. To avoid multicollinearity, we calculated the principal components (PCs) of a correlation matrix among the diet variables . Details of this analysis have been published elsewhere . The former approach allows for testing gradual shape changes in species groups with mixed diet. However, diet was also analysed as a discrete variable, using the predominant (or strict) diet item as category (see model-based analyses below).
Shape variation among species was summarised by principal components (the major axes of shape variation). The shape PCs were fitted to a phylogeny by estimating ancestral character states (nodes) for each PC and plotting the scores for OTUs (operational taxonomic units - the species at tree tips) and HTUs (hypothetical taxonomic units - the nodes) along with the branches linking the tree nodes. Each OTU and tree branch was plotted with a colour according to the predominant diet (estimated ancestral diet for the branches). A maximum likelihood algorithm  was used in the estimation of ancestral states, using the function ace in the APE package for the R environment for statistical computing [87, 84]. The phylogeny of Baker et al. (2003; 2010) [3, 33] (Figure 1), based on mtDNA and RAG2 sequences, was used as a framework for all the comparative analyses performed. The branch lengths used were the divergence time estimates obtained from Baker et al. (2010) . This phylogeny has been independently corroborated in a recent paper , although the time estimates are slightly different between the publications.
Principal components of species average shapes depict the major axes of variation among species in shape space. These linear combinations of the original superimposed coordinates ordinate species means according to gradients formed by unknown latent factors. Associations of PC scores with measured factors can be performed a posteriori to aid in the interpretation of ordinations. The linear associations between mandible shape and diet principal components were tested by multivariate regression (where each data point is a species average), using a phylogenetic generalized linear model (PGLS - [40, 88]). Generalized least squares (GLS) models allow for the phylogenetic non-independence to be incorporated into the error structure as an among species covariance matrix. Unlike spatially organised data (which share a similar problem of non-independence), where geographic distances among samples can be measured with small error, the phylogenetic covariances depend on the phylogenetic hypothesis they derive from and the evolutionary model assumed for character change. When phylogenetic covariances are based on the sum of branch lengths from the root of the tree to the most recent ancestral for each pair of species, a Brownian model of evolution is assumed, where the expected phenotypic differences among OTUs are directly related to their distance to the last common ancestor . The Brownian model of evolution has been increasingly considered an unlikely assumption for comparative methods (particularly when dealing with adaptive radiations [28, 38, 39]). We used the estimate of covariance structure proposed by Martins and Hansen , that results from an evolutionary model that incorporates stabilising selection as a constraint. According to this model, covariances among OTUs are determined as V ij = γ exp(-αt ij ), where α is a parameter determining the magnitude of a restraining force or pull towards an optimum phenotype, and γ is a parameter loosely interpreted as the interspecific variance in equilibrium (due to the balance between selective and neutral evolutionary forces). This model causes the influence of a common ancestor to decrease exponentially with the sum of branch lengths leading to that ancestor, as a consequence of stabilising selection "erasing" the evolutionary history. The value of α was interactively estimated via maximum likelihood during model fitting, providing flexibility of evolutionary models assumed. For the multivariate regression, the phylogenetic covariance matrix was estimated via the APE package  and the multivariate model fit with NTSYS-PC . The shape PCs were included as dependent variables and the first three diet PCs plus average cranial size (Condylobasal Length - CL) were included as independent variables. The partial correlation coefficients from the multivariate regression were used to produce biplots showing the species ordination within the shape PC space and the associated directions of variation of diet PC vectors.
Adaptive radiations commonly generate a distinctive pattern of temporal change in disparity, where morphological divergence increases quickly in the first phase, and slows down towards present time, as the niches are filled . The change in disparity patterns through time (DTT) is usually studied by the comparison of fossil taxa . However, Harmon and colleagues  have proposed and implemented (in the geiger package for the R environment ) a method to evaluate the changes in disparity patterns along a phylogenetic tree, using exclusively phenotypic measurements of terminal taxa. Using the full space of shape principal components (including all components with non zero eigenvalues) for the mandible shape, disparity was calculated from average pairwise Euclidean distances between species . For this data set, this is equivalent to using average pairwise Procrustes distances (summed squared distances among corresponding landmarks after Procrustes superimposition) . Using the full shape space is possible in this case because the data are reduced to pairwise distances regardless of dimensionality, and no parameters are estimated. However, the number of dimensions might have an effect in the simulations described below and the confidence intervals calculated, so it is more consistent to work with the full shape space. Size disparity was calculated using the same procedure, but using average Condylobasal Length as the single variable. Diet disparity was calculated as well, using the relative importance of dietary items (see Additional file 1) as data and average pairwise Manhattan distances as measure of disparity. Disparity through time (DTT) plots can be generated by first calculating the disparities for the entire tree and for each subclade defined by a node in the phylogeny. Relative disparities are obtained for each node by dividing its disparity by the tree disparity. Moving up from the root of the phylogeny, a mean disparity is calculated at each node (divergence event) as the average of the relative disparities of all lineages whose ancestral lineages were present at that point . The DTT plots allow for insight on the relative timing of phenotypic divergence while avoiding the need of ancestral character estimation. The inference of disparity deviations from a neutral evolution pattern was performed by a comparison of observed patterns with the average of 1000 simulated DTT patterns after evolving the shape variables according to Brownian motion along the phylogenetic tree [91, 93]. These simulations generate random phenotypes as phylogenetic independent contrasts along the tree (for each node) using a multivariate normal distribution with zero means and a covariance matrix based on the observed phylogenetic independent contrasts. A disparity distribution based on a neutral evolutionary process is then generated for each relative time (node) of the tree and the observed disparity can be compared to neutral confidence intervals.
A second approach to the investigation of evolutionary processes via a comparison of shape variation and relative time was based on Polly (2004)  simulation study. This author showed that different evolutionary models leave distinctive imprints on the distribution of shape distances versus time of divergence. We have compared the pairwise shape distances (Procrustes distances) to pairwise genetic distances (proportional to percent sequence divergence from Baker et al. (2003) ), and estimated time since divergence (from Baker et al. (2010) ) in a matrix scatterplot. This method allowed for a qualitative comparison of the observed pattern to the expected under different evolutionary modes by visual inspection .
The comparative methods based on GLS described above provide a direct test of linear relationships among major axes of diet and shape variation. Whereas this class of comparative methods might provide a statistically valid test of association, it provides no direct insight into the evolutionary processes responsible for the observed variation. A more recent approach, based on the simultaneous fit and comparison of evolutionary hypotheses combining evolutionary models, phenotypic and ecological information has been proposed by Butler and King (2004)  after Hansen's (1997)  method for the estimation of selective factors on adaptive optima (see also ). The model-based approach fits distributions predicted under a number of evolutionary hypothesis to observed data, allowing for simultaneous inference of model likelihood using information criteria.
The simplest model of phenotypic evolution is the neutral Brownian motion , where total evolutionary change in a continuous character is distributed normally, with a mean equal to the ancestral value for the entire tree and variance proportional to the amount of time elapsed from tree root to tips. An alternative evolutionary model (or class of models) for phenotypic change under stabilising selection that has gained popularity recently is based on the Ornstein-Uhlenbeck (O-U) process [27, 28, 40, 94]. These models incorporate a constraint as a constant pulling force (measured by a parameter α) towards a central point (θ) which can mimic the action of stabilising selection on the phenotype. The central point can be interpreted as an adaptive optimum. In the approach proposed by Butler and King , the O-U models can have any number of adaptive optima (θ i ), which are proposed according to hypotheses of ecological diversification. The models then incorporate a specific evolutionary process, a phylogeny with branch lengths (assumed to be at least proportional to time), and the selective regimes operative in each branch (derived from ancestral estimation of ecological changes).
where A is a p × p square symmetric matrix of α parameters that measure the strength of selection. the vector q is the optimum corresponding to the θ i for the particular branch, S is a p × p square symmetric matrix with sigma parameters that can be interpreted as the strength of random drift, and B(t) is a standard Wiener (Brownian motion) process . The alpha and sigma matrices have each p × (p + 1)/2 independent parameters (one for each variable and one for each pair of variables). These are estimated via an iterative optimisation procedure (we used identity matrices as starting forms) minimising the error of a multivariate GLS model that uses the alpha and sigma matrices to estimate the vector of optima q and the expected multivariate phenotype for each species (terminal branches) in the tree. Each model has to estimate 2 × p × (p + 1)/2 + p × n θ parameters, where n θ is the number of optima in the model. Because matrices of alpha and sigma parameters are calculated, the models do not constrain selection strength and the random drift disturbances to be isotropic in shape space. On the other hand, a simplification is accomplished because the same alpha matrix is used for all postulated adaptive optima, requiring the selective constraints to be the same for all peaks. This is true also for univariate models because they use a single α parameter for all adaptive peaks .
The ecological (dietary) information available allowed for the elaboration of five alternative models of evolutionary change (BM, OU.2, OU.3, OU.4, OU.5). The first model only incorporates neutral evolution expectations along the phylogenetic tree according to a Brownian motion (BM) process. The other models use the O-U process with a number of different optima (Figure 6) according to postulated selective regimes. We start with a simple discrimination of two adaptive optima related to high and low mastication regimes (Model OU.2). Model OU.3 breaks the low mastication selective regimes into nectarivory and sanguivory, recognising the fundamental ecological and functional differences between them. Model OU.4 separates all dietary optima, but lumps insectivory and carnivory together into the animalivory category. A justification for this is the strong correlation between insectivory and carnivory observed in the diet PCA . It has been suggested that carnivory is just an allometric extension of insectivory  and that these two categories share the same adaptive peak . Model OU.5 discriminates adaptive optima for five selective regimes associated with each dietary category, representing the most complex hypothesis about the adaptive radiation of phyllostomid bats.
for which ∑w i (AICc) = 1. The Akaike weight for a particular model M i can be interpreted as the probability of M i being the "best" model given the data and the set of candidate models . The same procedure was applied for the SIC to calculate SIC weights. The Schwarz information criterion is more conservative than the AICc, because it imposes larger penalties for more complex models. The SIC also assumes that the true model is within the candidate set of models, and is low dimensional, whereas the AICc assumes that the true model is not within the candidate set and is infinitely dimensional . Confidence limits (ellipsoids) for adaptive optima were estimated by a parametric bootstrap (resampled 5000 times). All model fitting procedures were performed using package OUCH (version 2.7-1)  for the R environment.
We are thankful to A. L. Peracchi (Universidade Federal Rural do Rio de Janeiro), M. de Vivo (Museu de Zoologia, Universidade de São Paulo, Brazil), and N. B. Simmons (American Museum of Natural History, EUA) for allowing access to specimens under their care, to J. Barros and E. Westwig for assistance in finding specimens in their respective institutions. Financial support was provided by FAPERJ and CNPq (Brazil) and The University of Hull (UK). A. Cardini, C. P. Klingenberg and an anonymous reviewer provided important comments and suggestions that greatly improved the manuscript. We also would like to thank L. Harmon for help and discussions about the package geiger.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.