Introduction

Adaptive radiations are characterized by the rapid diversification of species and divergence of new forms, often as a result of an evolutionary innovation (Simpson 1953; Schluter 2000; Yoder et al. 2010). As lineages move into a variety of previously unoccupied niches, they can experience high rates of morphological change, increasing clade-wide morphological disparity and functional disparity. The term ‘evolvability’ refers to the potential for diversification and describes the capacity of a population to evolve in the direction of selection by generating adaptive genetic diversity (Kirschner and Gerhart 1998; Hansen and Houle 2008). Highly evolvable populations are predicted to have the potential to speciate more rapidly in response to ecological opportunity (Young et al. 2010; Goswami et al. 2014). Evolvability can be measured quantitatively from morphological data through tests of morphological integration, which examine whether components of an organism, such as the cranium, are composed of one or more parts (Olson and Miller 1958; Mitteroecker and Bookstein 2007; Klingenberg 2008, 2013). While a tightly integrated cranium is composed of a single unit (e.g., the cranium), a less integrated cranium may be composed of multiple components, termed modules (e.g., rostrum, basicranium, cranial vault), which are individually tightly integrated, but only loosely integrated with other modules (Klingenberg 2009).

Integration is thought to affect phenotypic divergence in two different ways, one that requires loosely organized modules and another that operates within highly integrated modules. Theoretically, when structures are not tightly integrated and are composed of two or more modules, constituent modules are able to evolve independently of one another, increasing a taxon’s ability to adapt rapidly to environmental changes, which in turn drives an increase in species diversity and morphological shape disparity. A number of studies demonstrate this association empirically (Drake and Klingenberg 2008; Marroig et al. 2009; Young et al. 2010; Collar et al. 2014; Arlegi et al. 2018). For example, Young and Hallgrímsson (2005) found lower magnitudes of integration of forelimb and hind limb shape in bipedal mammals (bats and primates) in comparison with quadrupedal mammals. In this case, the authors hypothesized that decreased integration allowed limb proportions to escape phylogenetic and developmental constraints, explore new regions of morphospace (e.g., bipedalism), develop novel morphologies, and thereby increase overall morphological disparity.

In contrast, other studies identify clades with high rates of speciation and morphological evolution that are characterized by high levels of integration (Hu et al. 2016; Sherratt et al. 2017), possibly as a result of selection occurring along a single trajectory that may represent a line of least evolutionary resistance (Schluter 1996; Marroig and Cheverud 2005; Felice et al. 2018). For example, Hu et al. (2016) found high rates of speciation and morphological evolution coupled with high integration and neutral levels of disparity in icefish crania. They hypothesized that high integration constrained icefish crania to a single evolutionary trajectory, which decreased disparity but allowed for increased speciation relative to less integrated sister clades because that trajectory allowed increased differentiation in craniofacial height and length.

New World leaf-nosed bats (Family Phyllostomidae) are a well-known example of an adaptive radiation that was spurred, or at least accompanied by, morphological innovations (e.g., Freeman 2000; Jones et al. 2005; Dumont et al. 2011, 2014; Shi and Rabosky 2015; Hedrick and Dumont 2018). Phyllostomids are among the largest bat families and have a broader range of dietary ecologies than any other mammalian family (insectivory, folivory, frugivory, sanguinivory, nectarivory, and omnivory). With the exception of Old World fruit bats (Family Pteropodidae), all other bats are primarily insectivorous. Recently, Rossoni et al. (2019) used a quantitative genetic approach to assess the evolution of morphological integration within the crania of phyllostomids. They found overarching commonalities in patterns of integration but striking differences among species with very different morphologies. They also found that the magnitudes of integration decreased substantially during three periods of increased species diversification, suggesting a link between low integration and speciation. While this quantitative genetic approach is very different from the morphometric analyses typically used to investigate integration, it recovered the same signal of strong selection for divergent morphologies associated with dietary specialization found in previous studies (e.g., Freeman 1981; Dumont et al. 2011, 2014; Rossoni et al. 2017).

Rossoni et al. (2019) provided an informative assessment of integration and morphological disparity within and among phyllostomids, but the study did not address whether and how phyllostomids differ from their close relatives. Here we extend the analysis of integration and disparity in the shape of the cranium to all families within the Superfamily Noctilionoidea: Noctilionidae, Mormoopidae, Mystacinidae, Furipteridae, and Thyropteridae. We use three-dimensional geometric morphometrics to investigate three questions: (1) Do insectivorous non-phyllostomid noctilionoids and primarily insectivorous generalist phyllostomids evolve along separate evolutionary trajectories? Given the similarities in their diets, fundamental differences in shape between generalist phyllostomids and non-phyllostomid noctilionoids would suggest that factors beyond diet-driven functional divergence influenced morphological change within Noctilionoidea. We then evaluate disparity and integration in the Noctilionoidea and identify patterns to suggest how these factors have shaped noctilionoid evolution. (2) Is the shape of the cranium in phyllostomids more disparate than the shape of the cranium in other noctilionoids? Phyllostomids are often assumed to be more disparate because of their tremendous ecological and apparent morphological diversity. Nevertheless, few studies have compared them to other families of bats, especially with respect to noctilionoids as a whole (Freeman 1981, 2000; Hedrick and Dumont 2018). (3) What is the relationship between disparity in cranial shape and integration within phyllostomids and within other noctilionoids? We predict that high disparity will be coupled with low integration, and that this relationship will be expressed most strongly by phyllostomids.

Materials and Methods

Sampling and Data Collection

We selected bat species for study from across the New World noctilionoid tree that represent all constituent families of the superfamily Noctilionoidea (Phyllostomidae, Noctilionidae, Mormoopidae, Mystacinidae, Furipteridae, and Thyropteridae), as well as all major diet categories within the Phyllostomidae. These data included 141 individual bats from 65 species (Fig. 1). To account for individual variation, we included at least two specimens for most species. Several more common species were better sampled (e.g., for Artibeus jamaicensis, n = 5). Prior to all analyses, the mean shape of within-taxon individuals was taken and analyses were run on species rather than individuals, as in previous studies (e.g., Sherratt et al. 2017; Hedrick and Dumont 2018) due to low within-species cranial shape variation in noctilionoids (Dumont et al. 2014; Hedrick and Dumont 2018).

Fig. 1
figure 1

Pruned phylogeny of taxa included in the analyses from Shi and Rabosky (2015). Node colors refer to clade (as in Fig. 3a). Colored nodes indicate families (dark orange = Mystacinidae, green = Furipteridae, purple = Noctilionidae, yellow = Thyropteridae, pink = Mormoopidae, dark green = Phyllostomidae). Example cranium silhouettes shown to demonstrate high variability in cranium shape within the Noctilionoidea. Scale = 50 mm. From top to bottom (Noctilio albiventris, Mormoops sp., Macrotus waterhousii, Glossophaga soricina, Phyllostomus sp., Carollia sp., Centurio senex, Artibeus lituratus)

Full size image

Bat crania were scanned using a Nikon Metrology (X-Tek) HMXST225 MicroCT system at the Center for Nanoscale Systems at Harvard University. Image stacks were generated using proprietary software associated with the X-Tek scanner (CTPro, Nikon Metrology Inc., Japan), and 3D reconstructions created in VGStudio Max 3.0 (Volume Graphics Inc., Germany), and then segmented in Mimics v. 16.0 (Materialise, Leuven, Belgium). Crania were exported as PLY files and imported into Landmark Editor (Wiley et al., 2005) to generate three-dimensional landmark configurations (PLY files to be uploaded to Morphosource). Nineteen individual landmarks and two semi-landmark curves composed of 12 and 16 semi-landmarks, respectively, were used to characterize cranium shape (Fig. 2, see ESM 1 for a list of landmarks definitions and their biological relevance).

Fig. 2
figure 2

Three-dimensional scan of Artibeus jamaicensis with landmarks overlain in (a) lateral, (b) dorsal, (c) ventral, (d) rostral, and (e) caudal views (see ESM 1 for landmark definitions). Semi-landmark curves in red. All landmarks were placed on left side, but are visualized on the right in some cases here to more effectively show the landmark position

Full size image

Exploratory Analyses and Evolutionary Modeling

Landmark data were imported into RStudio v. 0.99.902 (R Core Team 2017) and opened in the geomorph package (Adams and Otárola-Castillo 2013). Data were subjected to Generalized Procrustes Analysis (GPA) with semi-landmarks slid according to the bending energy criterion (Perez et al. 2006). GPA translates, rescales, and rotates all specimens into a common shape space, minimizing all non-shape information (Zelditch et al. 2012). Following GPA, we took species means of all specimens using custom code and the mshape function in geomorph. An exploratory principal component analysis (PCA) was used to visualize how species and dietary regimes occupy morphospace. We ran an additional PCA with ancestral states and the phylogenetic tree mapped over species observations. We compared the cranium shape of insectivorous phyllostomids, non-phyllostomid insectivores, the most recent common ancestor (MRCA) of the Phyllostomidae, and the MRCA of the clade including phyllostomids and their sister family, the mormoopids.

Prior to examining the disparity and integration structure within the Noctilionoidea, taxa were generally placed into the dietary guilds following Dumont et al. (2014) to evaluate the effect of diet on cranium shape and specifically, whether non-phyllostomid noctilionoid insectivores and generalist phyllostomid insectivores grouped together. Dumont et al. (2014) split the Phyllostomidae into four adaptive regimes: nectarivores, generalists (primarily insectivorous), frugivores (stenodermatines excluding short-faced bats), and short-faced bats. Those analyses focused on adaptive regimes for mechanical advantage, particularly on the basis for the exceptionally high bite force of short-faced bats (Dumont et al. 2009). However, as we are more interested in the relationship between phyllostomids and their outgroups, we combined frugivores and short-faced bats into a fig-eating specialist category (Stenodermatinae). We also added: sanguivorous bats and non-phyllostomid insectivores. Therefore, we define five groups in total: [1] non-phyllostomid insectivores (k species = 12), [2] phyllostomid generalists (k = 23), [3] nectarivores (k = 12), [4] fig-eating specialists (k = 17), and [5] sanguivores (k = 1).

We examined the relationship between phyllostomids and their outgroups using three methods prior to evaluating integration and disparity to determine whether or not phyllostomids (particularly generalist phyllostomids) and non-phyllostomid noctilionoids evolved along different trajectories. First, a phylogenetic generalized least squares (PGLS) regression tested if phyllostomids and their outgroups had significantly different cranium shapes. Second, an exploratory canonical variates analysis (CVA) was conducted using the hypothesized adaptive regimes as groups in the R package Morpho (Schlager 2017) to establish whether the proposed groups separated based on cranium shape. We calculated the Mahalanobis distances between groups and assessed significant differences in groupings through permutation. To visualize differences between groups, canonical variate 1 (CV1) was plotted against canonical variate 2 (CV2) and a dendrogram based on Mahalanobis distances was generated. Classification success was estimated using a discriminant function analysis. Finally, evolutionary modeling analyses were then used to assess the effects of diet on cranium shape evolution.

To evaluate how noctilionoid cranium shape evolved and the nature of the relationship between the hypothesized adaptive regimes and cranium shape, we examined seven different evolutionary models (ESM 2). Sanguivores were only represented by a single taxon in this dataset, Desmodus rotundus, and hence were not included in evolutionary modeling analyses. The first model was Brownian motion (BM), which assumes that evolution follows a random walk and is not drawn to any particular optimum. This represents the null hypothesis of no directed evolutionary change. Multiple Ornstein-Uhlenbeck (OU) models were also examined. OU models describe evolutionary change in a continuous trait, in this case, principal component scores of cranium shape, as a deterministic component describing the strength of selection to one or more optima (Hansen and Houle 2008; King and Butler 2009). This optimum is characterized by θ, with α representing the strength of the attraction toward that optimum, and σ representing the intensity of random walk fluctuations. A single-peak OU model represents stabilizing selection, with change drawn to a single optimum, whereas multi-peak OU models hypothesize evolution towards different adaptive optima.

OU peaks were chosen in accordance with the four hypothesized adaptive regimes: phyllostomid generalists θG, fig-eating specialists θF, nectarivores θN, and non-phyllostomid generalists θNP. Three two-peak OU models were examined. OU2 (Fig-eating bats) split fig-eating specialists from all other groups to test whether these bats – which have been shown to have an accelerated rate of speciation relative to other bats (Dumont et al. 2011; Shi and Rabosky 2015) – represent an optimum within the Noctilionoidea while other groups converge on a single optimum (θF ≠ θN = θG = θNP). We ran a similar model for nectarivores (OU2 (Nectarivores), θN ≠ θF = θG = θNP). Nectarivorous bats have developed some of the most derived bat skull shapes and demonstrate the only instance of convergent evolution to an adaptive regime within the dataset. Finally, we examined non-phyllostomid noctilionoid insectivores (OU2 (Outgroups), θNP ≠ θN = θG = θF). Non-phyllostomid noctilionoid insectivores use oral instead of nasal echolocation, and phyllostomids may represent an adaptive radiation, unlike non-phyllostomid noctilionoid outgroups (Freeman 2000). Therefore, it is of interest to determine if phyllostomids and other noctilionoids have separate optima within Noctilionoidea. Finally, we tested our hypothesis that diet is the primary driver of shape directly with a model in which fig-eating specialists and nectarivores represent separate optima while phyllostomid generalists and non-phyllostomid insectivores converge on the same optimum (OU3, θF ≠ θN ≠ θG = θNP). This would suggest that basal phyllostomids and generalists have cranium shapes similar to other noctilionoids. We then tested if all four groups diverge to separate optima (OU4, θF ≠ θN ≠ θG ≠ θNP).

Evolutionary modeling was carried out in the R package ouch (King and Butler 2009). These models were run for PC1 and PC2 separately. Models were fit using the brown (BM) and hansen (OU) functions, and Akaike Information Criteria modified for sample size (AICc) were calculated to assess model fit.

Disparity Analyses

Disparity was calculated for the phyllostomids and their outgroups separately both by estimating the Procrustes variance of each group without phylogenetic corrections and by calculating the Procrustes distance for each taxon from the mean shape of the data using phylogenetic corrections. Significance between group disparities for Procrustes variance was assessed with 999 permutations (Zelditch et al. 2012). More permutations did not alter the significance of the test. A phylogenetic ANOVA in the package phytools was used to compare Procrustes distances among groups (Revell 2012).

Modularity and Integration Analyses

The mammalian cranium has been hypothesized to be composed of two primary modules across the mammalian tree, a rostral and a caudal module (Porto et al. 2009). This is the case even in bats that have developed cranial elaboration for nasal echolocation (Santana and Lofgren 2013). These studies suggest that rather than differences in the number of cranial modules, differences in magnitudes of integration between the two modules drive mammalian cranial variation (Porto et al. 2009). Landmarks 1, 3, 4, 5, 8, 10, 12, 16, 17, and 20–37 were included in the rostral module and landmarks 2, 6, 7, 9, 13, 14, 15, 18, 19, and 38–47 were included in the caudal module (ESM 3). Additionally, a zero-module hypothesis and two-module hypothesis were tested using EMMLi (Goswami and Finarelli 2016). EMMLi uses a maximum likelihood approach to compare module patterns using Akaike Information Criteria (AICs) rather than testing a single hypothesis of modularity a priori. We note that several other studies have found that mammal crania are made up of as many as six modules (Cheverud 1995; Goswami 2006; Goswami and Finarelli 2016). Our landmark coverage was designed to test two modules and therefore higher-level modularity hypotheses were unfortunately not possible with our dataset.

After assessing modularity and prior to comparing levels of between-module integration among groups, we ran phylogenetic partial least squares (PLS) analysis (Adams and Felice 2014) for all taxa, and for the phyllostomids and their outgroups separately. Data for the rostral block were permuted 500 times across the tips of the phylogeny and the resulting distribution was compared to the observed covariation between the rostral and caudal cranium blocks to assess significance. These data were also compared with the results from the EMMLi analysis, which assesses within- and between-module integration through the rho statistic. Following this, we compared the levels of between-module integration between the phyllostomids and their outgroups by calculating the effect sizes of the two PLS analyses using the compare.pls function in geomorph (Adams and Collyer 2016). We additionally compared the levels of between-module integration within the Phyllostomidae by comparing stenodermatines to non-stenodermatine phyllostomids. Rather than using the PLS correlation coefficient, this analysis does two-sample z-tests, which are robust to differences in sample size and number of landmarks.

The datasets generated and/or analyzed during the current study (including PLY files) are available in the Morphosource repository.

Results

Exploratory Analyses and Evolutionary Modeling

In principal component (PC) morphospace, PC 1 accounts for 43.2% of total variance and PC 2 accounts for 33.8% of total variance. PC 3 and subsequent PCs account for 5% of total variance or less and are not considered further in this study (ESMs 4, 5). PC 1 primarily divides non-phyllostomid taxa from phyllostomids (Fig. 3a) while PC 2 distinguishes different dietary modes within phyllostomids (Fig. 3b). Taxa with strongly positive PC 1 scores (e.g., Mormoops spp.) have an upturned rostrum, while taxa at the negative end of PC 1 have flattened rostra similar to that of many phyllostomid generalists (e.g., Micronycteris spp.). Thyroptera, Mystacina, Furipterus, and Noctilio all plot towards the center of morphospace between the more exaggerated rostra of mormoopids and the flattened rostra of phyllostomids (Fig. 3c). Taxa with strongly negative PC 2 scores have extremely shortened rostra as exemplified by Centurio and Sphaeronycteris. Nectarivorous taxa with elongate rostra characterize the positive end of PC 2. Non-phyllostomids all have intermediate PC 2 values with non-phyllostomids showing neither the extreme rostral elongation nor the rostral shortening found in phyllostomids. The morphospace thus represents two qualitatively different trajectories, one leading from the center of morphospace to derived mormoopids (along PC 1) and another along PC 2 dividing ecological groups within the Phyllostomidae. Finally, a PGLS of cranium shape showed that phyllostomids and their outgroups had significantly different cranium shapes (p = 0.002), but group membership explained a relatively small amount of total shape variation (R2 = 0.054; Table 1).

Fig. 3
figure 3

a Principal component analysis (PC1 x PC2) of all species (n = 65). Each point represents a single species. Species for which multiple specimens were scanned are represented by the mean of the specimens. Colors reflect node colors in Fig. 1. Note that the basal insectivorous phyllostomids (dark green) plot away from all insectivorous non-phyllostomids. Wireframes show major changes along the first principal component axis. b As Fig. 3a, but species are colored based on hypothesized dietary regimes (see branch colors in Fig. 1). c Cartoon morphospace showing overall trends in cranium shape using sample taxa. d Canonical variates analysis of major dietary regimes. All groups are fully separated. Inset shows dendrogram of Mahalanobis distances

Full size image
Table 1 Phylogenetically informed Procrustes ANOVA of shape and group (Phyllostomidae, non-phyllostomid)
Full size table

The CVA fully distinguishes all five dietary groups in morphospace (Fig. 3d). Canonical variate (CV) 1 summarizes 41.3% of total variance and distinguishes non-phyllostomid insectivores from phyllostomid groups (fig-eating specialists, nectarivores, phyllostomid generalists, Desmodus). CV2 summarizes an additional 31.6% of total variance and distinguishes nectarivores and fig-eating specialists from Desmodus. Based on Mahalanobis distances, Desmodus is farthest away from other groups. Phyllostomid generalists, fig-eating specialists, and nectarivores are more similar to one another than to non-phyllostomid insectivores (Fig. 3d). Furthermore, a discriminant function analysis shows all five groups significantly differ from one another (ESM 6).

Seven adaptive regimes were modeled and AICcs were used to compare between models (Table 2; ESM 2). For PC1, the best-supported model is the OU4 model. The OU2 (Outgroups) model has a ΔAICc of 3.67 suggesting it was relatively well supported, but less so than the OU4 model. For PC2, the BM, OU1, OU2 (fig-eating specialists and others), and the OU2 (phyllostomids and outgroups) models are all equally well supported with OU4 being poorly supported (Table 2). The PCA, CVA, and evolutionary modeling analyses, combined, strongly suggest that while the outgroup taxa and generalist phyllostomid taxa predominantly eat insects, they have diverged in cranium shape.

Table 2 Evolutionary models (see ESM 2) showing AICcs for PC1 and PC2. For PC1, the OU4 model fit the data best. For PC2, the BM, OU1, OU2 (fig-eating and others), and OU2 (phyllostomids and outgroups) models were all equally successful
Full size table

Disparity Analyses

Analyses compared levels of disparity in phyllostomids and in their outgroups. Despite having fewer taxa represented (53 phyllostomids versus 12 outgroup taxa), the outgroups have significantly higher shape disparity, as measured by Procrustes variance, than phyllostomids. Indeed, shape disparity is more than three times higher among outgroup taxa (Procrustes variance = 0.0381) than phyllostomids (Procrustes variance = 0.0112, p = 0.004). In contrast, phylogenetically corrected Procrustes distance suggested that the two groups did not have significantly different levels of disparity (p = 0.139). However, Procrustes distance between the two groups is significantly different without phylogenetic corrections (p < 0.001) suggesting phylogeny plays a strong role in structuring disparity in noctilionoids. Further, phyllostomids have a mean Procrustes distance of 0.094 while outgroups have a mean Procrustes distance of 0.178 mirroring the same pattern suggested by Procrustes variance.

Modularity and Integration Analyses

The EMMLi analyses (Goswami and Finarelli 2016) suggest the crania are modular with the two-module hypothesis more strongly supported than no modularity (Table 3). EMMLi further shows between-modules integration is high (rho = 0.4), but lower than within-module integration (module 1: rho = 0.55; module 2: rho = 0.63). The rostral and caudal aspects of the cranium are significantly correlated (PLS correlation coefficient = 0.862, p < 0.001) for all species (Fig. 4a). Block one (rostral component) changes from an elongated rostrum (e.g., nectarivorous phyllostomid) to a foreshortened rostrum (e.g., short-faced bat). Block two (caudal component) displays minimal changes, but those present relate to the shortening of the cranium roof relative to the ventral aspect of the basicranium. The rostral and caudal aspects of the cranium are also significantly correlated when the outgroups are evaluated on their own (PLS corr. = 0.964, p < 0.001). In the outgroups-only analyses, changes in both the rostral and caudal aspect of the cranium relate to how upturned the rostrum is with taxa such as Mormoops spp. having more negative PLS scores for both blocks, and taxa with less upturned crania having more positive PLS scores for both blocks (Fig. 4b). Results of analysis of the phyllostomid-only dataset mirror those of the all-species between-module integration analysis in the changes that are displayed along blocks one and two (Fig. 4c). As in the other analyses, block one and two are also significantly correlated for the phyllostomids-only dataset (PLS corr. = 0.873, p < 0.001).

Table 3 Modularity and integration results for non-phyllostomid noctilionoids and phyllostomids separately
Full size table
Fig. 4
figure 4

Phylogenetic partial least squares analysis comparing degree of integration of the rostral and caudal portions of the cranium for (a) all specimens, (b) only non-phyllostomid noctilionoids, and (c) only phyllostomids. Warp grids show gradation in shape from positive to negative values for each block. PLS correlation coefficients and p values are inset

Full size image

Although the PLS correlation coefficient is higher in outgroup taxa than phyllostomids in the within-group PLS analyses, between module integration effect size is nearly twice as high in phyllostomids than in the outgroup taxa (Table 3). This difference is significant (p = 0.024). Within the Phyllostomidae, non-stenodermatine phyllostomids have a higher integration effect size (7.038) than stenodermatines (3.563). However, this difference is not significant (p = 0.27).

Discussion

Although a number of studies have examined the relationship between morphological disparity, species richness, and the magnitude of integration with theoretical modeling (Klingenberg 2005; Marroig et al. 2009; Goswami et al. 2014; Randau and Goswami 2017; Felice et al. 2018), the relationships among these factors have rarely been investigated with empirical data (but see Hu et al. 2016; Sherratt et al. 2017; Felice et al. 2018). To empirically test hypotheses regarding these relationships, we generated a 3D cranium dataset for the Phyllostomidae and outgroups within the superfamily Noctilionoidea. Phyllostomidae is one of the most species-rich families within Chiroptera (Shi and Rabosky 2015), and, with an unparalleled number of dietary types of any mammalian family, is often cited as a textbook example of an adaptive radiation (Jacobs 2016). This group therefore provides an excellent opportunity for investigating how integration, or the lack thereof, can drive such radiations. Theoretical studies have posited that increased integration leads to evolution along a single trajectory into a more limited morphospace, whereas decreased integration leads to greater morphospace occupation and thus higher disparity (Goswami et al. 2014; Felice et al. 2018). Likewise, several empirical studies have suggested high integration leads to evolutionary constraints (Young and Hallgrímsson 2005; Marroig et al. 2009; Young et al. 2010; Drake and Klingenberg 2008; Collar et al. 2014; Arlegi et al. 2018), while others have found that high levels of integration can lead to niche expansion and increases in morphological disparity (e.g., Hu et al. 2016) via evolution along a single trajectory, potentially representing a line of least evolutionary resistance. In this study, we found decreased morphological disparity and increased integration in the crania of phyllostomids relative to those of their outgroups in spite of the relatively high rates of speciation, high species richness, and greater functional breadth that characterize the Phyllostomidae (Dumont et al. 2014; Shi and Rabosky 2015).

Ecological shifts have almost certainly influenced phyllostomid evolution, as the magnitude of selection in phyllostomids has been shown to be strongest on branches leading to novel feeding habits, followed by a decrease in the magnitude of selection after novel dietary niches have been occupied (Rossoni et al. 2017). Therefore, before evaluating relative levels of integration and disparity among phyllostomids and their outgroups, we examined whether primarily insectivorous, generalist phyllostomids differed from their obligatory insectivorous sister groups in cranium shape. If this were the case, then it would suggest that factors beyond functional divergence (such as differing magnitudes of integration) played a role in the generation of the morphological differences within Noctilionoidea. We identified two qualitatively different trajectories in morphospace separating phyllostomids from their outgroups. While the crania of insectivorous phyllostomids are closest in shape to those of outgroups, they remain distinct (Fig. 3b, d), a result supported by PGLS. Furthermore, the insectivorous MRCA of the Phyllostomidae plots away from the outgroup morphospace and within the phyllostomid morphospace (ESM 7) suggesting a shift in cranium shape likely occurred as the ancestor of extant phyllostomids evolved. We also found evidence that phyllostomid generalists and their insectivorous outgroups evolved to occupy separate adaptive peaks (Table 2). Although the morphospace and adaptive peak analyses together suggest that morphological change in the Noctilionoidea has responded to other factors besides selection for dietary function, they do not elucidate what those factors may have been. To gain further insights into these factors, we analyzed the levels of integration and disparity in their crania.

Despite higher rates of speciation and greater functional breadth among phyllostomids (Dumont et al. 2014; Shi and Rabosky 2015), their crania have lower morphological disparity and greater integration than those of their outgroups (Fig. 4, Table 2). These results seem, on the surface, contrary to the observation that the Phyllostomidae represents one of the most functionally diverse chiropteran families (Freeman 2000), and that stenodermatine phyllostomids are the only clade within Chiroptera to show increased rates of speciation. We had predicted that phyllostomids would exhibit increased disparity relative to their less species-rich outgroups (although previous studies on the phyllostomid palate found high levels of integration––Sorenson et al. 2014). Other studies have found a reduction of between-limb integration and integration of cervical vertebrae in humans relative to great apes, which may have contributed to the emergence of bipedalism (Young et al. 2010; Arlegi et al. 2018), decreased limb integration in some non-quadrupedal mammals compared to some quadrupedal mammals (Young and Hallgrímsson 2005), and a decrease in magnitudes of integration in some biting eels allowing for more varied morphologies as a result of a relaxation of constraints (Collar et al. 2014).

While some authors have found evidence that reduced integration may not correlate with functional specializations (Botton-Divet et al. 2018), others have also recently found that high levels of integration can be coupled with high levels of speciation and morphological evolution (Hu et al. 2016; Sherratt et al. 2017), suggesting that increased evolutionary potential may be achieved via a single axis of evolutionary change (Schluter 1996; Marroig and Cheverud 2005; Goswami et al. 2014; Felice et al. 2018). Still others have found that high levels of integration do not necessarily limit or increase morphological evolution over long time scales (Goswami and Polly 2010). These empirical findings have been supported by some theoretical modeling studies as well. For example, Goswami et al. (2014) noted that high magnitudes of integration evoke large responses to selection, albeit along paths determined by trait covariances such that perfectly integrated structures evolve along a single axis. These authors further showed that strong integration may limit morphological disparity and can lead to the evolution of both extreme morphologies and convergence in morphology. These findings raise the question of whether a similar case might apply to phyllostomids whereby phyllostomids have reduced overall disparity in comparison with their outgroups, but have evolved extreme morphologies and multiple convergent shapes.

Cranial morphology primarily differs among phyllostomids along a single axis based on the relative length of the rostrum. Taxa with highly positive PC2 values have elongate rostra (e.g., Choeroniscus minor) and taxa with highly negative PC2 values have foreshortened rostra (e.g., Centurio senex) (Fig. 3). Morphological variance among phyllostomids along PC1 is extremely low. Instead, phyllostomids have evolved extreme elongation and extreme shortening of the rostrum relative to other bats, allowing them to access novel diets within Chiroptera. Finally, although there are not many examples of convergence within the Phyllostomidae, nectarivorous bats do converge in morphology (Datzmann et al. 2010; Rojas et al. 2016). Hence the Phyllostomidae appear to follow quite closely the predictions made by Goswami et al. (2014).

Although we now have a substantial amount of data concerning phyllostomid cranial morphology (Nogueira et al. 2009; Rossoni et al. 2017; Hedrick and Dumont 2018), the genetic and developmental mechanisms at the origin of their diversity remain poorly understood. Many candidate genes are linked to variation in face length based on mouse mutants (reviewed in Usui and Tokitam 2018) including BMP2, BMP4, and Msx1 (Bonilla-Claudio et al. 2012). In addition, some genes such as Wnt5a have been shown to directly regulate crucial stages in face formation such as the migration of the ectomesenchyme (Kaucka et al. 2016). These genes are suspected to be potential candidates for the origin of the short face for shorted-faced bats (Bonilla-Claudio et al. 2012). Recently, variation in amino-acid composition of the transcription factor Runx2 (Runt-related transcription factor 2) implicated in bone differentiation, has been correlated to variation in face-length in carnivorans, primates, and phyllostomids (Sears et al. 2007; Ritzman et al. 2017; Ferraz et al. 2018). However, this correlation remains indirect and Runx2 likely acts in concert with other factors. It is possible that changes in some of these genetic and developmental mechanisms have led to the dietary diversity and species richness in phyllostomids without decoupling the rostrum from the back of the cranium. By identifying distinct groups in terms of face-length and morphology, our study could serve as a starting point to investigate which changes in the developmental program drive face-length variation.

In addition to diet, another major difference between phyllostomids and other noctilionoids is their mode of echolocation. Phyllostomids are primarily nasal echolocators, albeit with some plasticity (Griffin and Novick 1955; Schmidt 1988), whereas non-phyllostomid noctilionoids are oral echolocators. This difference, therefore, provides a potential confounding variable in our study (Pedersen 1998). However, it appears that facial measurements do not strongly correlate with echolocation parameters when comparing nasal and oral echolocators (Goudy-Trainor and Freeman 2002). Expanding empirical data for integration magnitudes and disparity levels beyond the Noctilionoidea will be a fruitful area of future research.

This study builds upon the assessment of modularity and integration in the mammalian cranium providing support for a two-module system in comparison with a fully integrated cranium. Although the two-module hypothesis is often supported (e.g., Porto et al. 2009; Santana and Lofgren 2013), other studies have found stronger support for a six-module mammalian cranium (e.g., Cheverud 1995; Goswami 2006; Goswami and Finarelli 2016). Although we did not test hypotheses of more than two modules because there were few landmark locations, it remains possible that the internal structure of the cranium contains additional internal modules not captured in our analyses. The highly integrated external form may allow for high bite force, while decreased levels of integration of internal structures such as sensory structures (e.g., cochlea, eye globe) lead to functional divergence (Goswami and Finarelli 2016). Future work should analyze individual sensory structures in the head using a combination of both internal and external landmarks to assess whether finer-scale modules exist within the phyllostomid cranium.

In conclusion, we propose that the tightly integrated crania of phyllostomids were capable of evolving into unique dietary ecologies along a single trajectory in morphospace, which possibly represents a line of least evolutionary resistance. Although phyllostomid cranial morphology only expanded along a single axis and did not spread into as large a region of morphospace as their immediate outgroups, evolving along a single axis allowed phyllostomids to rapidly evolve exaggerated morphologies and convergent shapes capable of making use of unique diets. This study builds on many previous theoretical and empirical studies examining modularity in the mammalian cranium (Cheverud 1995; Goswami 2006; Drake and Klingenberg 2008; Porto et al. 2009; Marroig et al. 2009) and provides one of the best empirical examples of Goswami et al. (2014) and Felice et al. (2018)‘s predictions of the relationship between morphological disparity, the magnitude of integration, and evolvability, expanding our understanding of how integration, and not just function, can lead to unique morphologies as well as potentially spur adaptive radiations.