Cranial and mandibular shape variation in the genus Carollia (Mammalia: Chiroptera) from Colombia: biogeographic patterns and morphological modularity

Sample sites and specimen selection

A total of 286 specimens of Carollia (C. brevicauda = 108; C. castanea = 82; C. perspicillata = 96) from 143 different localities in Colombia were evaluated for this study (see Table S1). The criteria for specimen selection were: that only sites with at least one male and one female available were considered, and, to ensure adequate representation of all five biogeographic regions (Caribbean, Pacific, Andean, Amazonian, and Orinoquean) and independence between samples (sites separated by at least 30 km), that one locality only was selected per municipality for each species (Fig. 1).

All specimens were obtained from the Instituto Alexander von Humboldt (IAvH-M), Colección Teriológica de la Universidad de Antioquia (CTUA), Instituto de Ciencias Naturales de la Universidad Nacional de Colombia (ICN) and the Museo Javeriano de Historia Natural (MPUJ).

Morphological analysis

Photographs were taken with a Nikon D5100 mounted on a tripod; crania were photographed in ventral view and mandibles in lateral view. In order to optimize and standardize the photographs, focal distance was estimated using the method proposed by Blaker (1976) and different holders were used for crania and mandibles.

Following geometric morphometric principles, landmark configurations were established for crania and mandibles separately using type 1 and 2 landmarks (Bookstein et al., 1985). Modifying the methodology used by Jarrín & Menendez-Guerrero (2011), a total of 15 landmarks were used for the cranium (Fig. 2A), and following previous studies (Zelditch, Wood & Bonett, 2008; Jojić, Blagojević & Vujošević, 2012) 12 landmarks were used for the mandible (Fig. 2B) (see Table S1). Landmark digitalization was performed using TPSDIG version 2.16 (Rohlf, 2010).

Generalized Procrustes Analysis (GPA) was performed in order to superimpose landmark coordinates, obtaining the average coordinates of all landmarks in a tangent configuration; this was performed separately for the cranium and mandible datasets (Rohlf, 1990). GPA removes non-shape sources of variation resulting from scaling, rotation and translation (Rohlf, 1999). A tangent configuration is the configuration of landmarks projected from a nonlinear shape space into a tangent space in which parametrical statistical analysis can be performed. Using TPSRELW (Rohlf, 2010), a Relative Warp Analysis (RWA) was performed following the principle of the thin-plate spline technique, which allows the partition of the total variation among all specimens from the tangent configuration in two different components: affine components that describe differences in uniform shape variation (principal warps), and non-affine components that express local variation within the shape (partial warps) (Rohlf, Loy & Corti, 1996).

Relative Warps (RW) are the principal components of a distribution of shapes in a tangent space, comprising the majority of the variation in a few comprehensive components, which are easily visualized using a transformation grid (Rohlf & Bookstein, 2003). RW are non-biological variables used as a representation of affine and non-affine components that describe localized deformations in specific regions of the overall shape, and can be analyzed using conventional statistical methods (Klingenberg, 2013). RW were computed using the partial warps for further statistical analysis.

Patterns of interspecific variation

Interspecific differences in the intraspecific morphological variation were tested with a multivariate analysis of variance (MANOVA) and a paired Hotelling’s test using the RW pooled by species; these analyses were performed using PAST version 2.15 (Hammer, Harper & Ryan, 2001). Squared Mahalanobis distances were used as a measure of morphological distances between species to assess general patterns of variation for all species, P values were corrected with a Bonferroni correction for multiple comparisons α at = 0.05.

In order to detect specific regions where major morphological variation may be focused, RW were visualized using transformation grids for each species, comparing the morphological patterns of variation between each species for the cranium and mandible (Zelditch et al., 2004). Patterns of shape change were depicted using TPSRELW (Rohlf, 2010), and the grids were built with the Principal Components (PC) of the Procrustes coordinates using MORPHO J version 1.04a (Klingenberg, 2011).

Cranial-mandibular integration and modularity

Based on previous findings of functional modularity in mammals (Zelditch, Wood & Bonett, 2008; Monteiro & Nogueira, 2011; Jojić, Blagojević & Vujošević, 2012), two different a priori hypotheses were considered for evaluating morphological modularity, one for the skull and one for the mandible (Fig. 3). The first divided the skull into two functional modules, neurocranium (muscle insertion and brain development) and splachnocranium (feeding and biting behavior); the second divided the mandible also into two functional modules, the ascending ramus (muscle insertion) and the alveolar region (supporting the teeth).

These hypotheses were evaluated with the Escoufier’s RV coefficient using MORPHO J version 1.04a (Robert & Escoufier, 1976; Klingenberg, 2009; Klingenberg, 2011). This method takes the RV coefficients of the a priori hypothesis and compares it with coefficients of multiple alternate partitions, and hypotheses with coefficient values closer to zero are not rejected. Delaunay triangulations were considered during module construction among landmarks (Berg et al., 2000). For this study we set 10,000 alternate partitions to compare with each a priori hypothesis, and this procedure was applied for each species.

Studying cranial-mandibular integration allowed us to evaluate whether the cranium and mandible together behave as a functional unit, covarying morphologically in their shape (Klingenberg, 2008). To do this, partial least square analysis (PLS) was performed, which explores patterns of covariation between different blocks of variables. RW were pooled by structure (cranium and mandible) and species, performing a PLS for all species where cranium and mandible shape were assigned as different blocks; this analysis was performed using TPSPLS version 1.18 (Rohlf & Corti, 2000).

Geographic patterns vs. morphological variation

RW of each species were pooled, differentiating biogeographic regions (Caribbean, Pacific, Andean, Amazonian, and Orinoquean); this was done for the cranium and the mandible separately. MANOVA and paired Hotelling’s tests were used to assess morphological differences between populations from different geographic regions, and to test if the Andes represent a barrier that divides morphological differences among populations of the same species, separating populations of different biogeographic regions morphologically. P values were corrected with a Bonferroni correction for multiple comparisons α at = 0.05.

Patterns of interspecific variation

The MANOVA differences between species were significant for both cranium (λ = 0.5185; df1 = 18; df2 = 548; F = 11.84; p = 3.05E–26) and mandible (λ = 0.5966; df1 = 18; df2 = 546; F = 8.937; p = 1.68E–1) data sets. All pairwise comparisons were significant with P-values ≤ 1.05E–03.

Based on squared Mahalanobis distances we found C. castanea to be the most morphologically different species, being most distinct in its cranial morphology from the rest of the species (Table 1).

As a general pattern, for all species the majority of the variation was concentrated in the neurocranium, around the suture of the occipital and temporal bones, as well as the area comprising the vomer and the palatine (Fig. 4). Each species showed species-specific variation patterns within these regions (Figs. 4A–4C).

For C. brevicauda, the highest deformation in the neurocranium is displaced towards the mastoid due to a constriction of the occipitomastoid suture and the tympanic part of the temporal bone (Fig. 4A). On the other hand, C. castanea exhibited an expansion in the region of the suture towards the occipital and a reduction of the length of the vomer (Fig. 4B). Finally, morphological variation in C. perspicillata was evident in the basicranium, between the foramen magnum and the vomer, and at the occipital and temporal bones. Variation in both regions showed a general contraction of such bones, leading to a general reduction in the length of the neurocranium (Fig. 4C).

Regarding mandibular morphology, the three species varied in the same regions, but the way in which they varied was different between species. Most interspecific variation was concentrated in the middle region of the ascending ramus and the alveolar region (Fig. 5).

When comparing variation across species, C. brevicauda showed greater variation in the lower border of the ramus, between the condyloid and angular processes (Fig. 5A); for C. castanea the mandibular tooth row and the base of the ramus expanded, resulting in a constriction of the medium region between the ascending ramus and the alveolar region (Fig. 5B). Carollia perspicillata showed the same pattern in the lower border of the ramus, but in this case the mandibular tooth row was shortened, in contrast to C. castanea (Fig. 5C).

Cranial-mandibular integration and modularity

All a priori hypotheses for functional mandibular modularity were not rejected as they had the lowest RV coefficients, dividing the mandible into two different modules (ascending ramus and alveolar region) according to their functional specializations (Table 2). However, 913 different partitions, including the a priori hypothesis, are compatible with the data for C. castanea, which could mean that, although the evaluated hypothesis was not rejected, there are other factors that affect mandibular modularity in this species. Results indicate the a priori hypothesis for cranial modularity was rejected in all cases, finding alternate partitions with lowest RV coefficients (Table 2).

Partitions recovered for mandibular modularity had the same structure for all species. Similarly, for cranial modularity, the same general partition pattern, dividing the cranium into two modules representing the neurocranium and the splachnocranium, was recovered. However, the structure of these partitions varied between species, each species having different modularity patterns, and such differences being present in the sphenoidal section of the basicranium (Figs. 6A–6C). Cranial modularity results for C. brevicauda showed that the neurocranium module comprises the zygomatic process of the temporal bone (landmarks 3–10), while the splachnocranium module comprises the palatine (landmarks 9–10) and vomer bones (landmarks 3–9) (Fig. 6A). For C. perspicillata the neurocranium module comprises the zygomatic process of the temporal bone and the vomer and the splachnocranium module comprises the palatine (Fig. 6C). Carollia castanea showed the most distinct modularity patterns where the neurocranium module extends anteriorly covering the zygomatic process of the temporal and the posterior section of the palatine, while the splachnocranium module extends posteriorly covering the vomer (Fig. 6B). For all species the first three dimensions of the PLS explained around 80% (C. brevicauda 78.32%, C. castanea 84.84% and C. perspicillata 76.91%) of the cranial-mandibular morphological integration, R values were always positive (ranging from 0.37 to 0.65), and the coefficient of determination (r²) values corroborated the significance of the results (Table 3).

Geographic patterns vs. morphological variation

MANOVA results were not significant for morphological differences in the mandible between specimens of the five biogeographic regions in any of the species; C. brevicauda (λ = 0.6674, df1 = 36, df2 = 354, F = 1.120, P = 0.2971), C. castanea (λ = 0.9243, df1 = 9, df2 = 72, F = 0.6552, P = 0.7461), C. perspicillata (λ = 0.6555, df1 = 36, df2 = 309, F = 1.024, P = 0.4358). For the skull, MANOVA found significant differences between biogeographic regions in C. brevicauda (λ = 0.5182, df1 = 36, df2 = 357.7, F = 1.906, P = 0.0018) and C. perspicillata (λ = 0.3862, df1 = 36, df2 = 309, F = 0.5179, P = 0.0469), and no significant difference was found for C. castanea (λ = 0.5468, df1 = 36, df2 = 260.3, F = 1.264, P = 0.1538).

Paired Hotelling’s test supports these results, finding significance only for C. brevicauda and C. perspicillata. Results revealed that for both species only specimens from the Andean region were different from the rest; Andean specimens of C. brevicauda were statistically different from Amazonian and Caribbean specimens, and for C. perspicillata Andean specimens were different from those of the Pacific region (Table 4).

Patterns of interspecific variation

Results confirmed that despite the presence of intraspecific variation in all species, the mode of this variation differs between species (Jarrín, Flores & Salcedo, 2010). Among these, C. brevicauda and C. perspicillata (larger species) are most similar, and C. castanea (smaller species) is the most divergent (Table 1). This is consistent with phylogenetic analysis in this genus that shows that C. brevicauda and C. perspicillata are sister species and the most recently diversified, while C. castanea is the oldest species (Hoffmann & Baker, 2003).

Previous studies have shown that major cranial morphological variation in these species is present in the neurocranium, specifically in the region that comprises the occipital bone and the squama portion of the temporal bone (Jarrín & Menendez-Guerrero, 2011), supporting our findings of major cranial morphological variation in the occipital and temporal bones (Fig. 4). Quantifying differences in dietary specialization and breadth between species (Dumont, 1999), as well as the specific characteristics of consumed items, such as object hardness and size, could shed some light on the mechanisms shaping the differences found in the patterns of intraspecific variation (Dumont & Piccirillo, 2005).

In phyllostomid bats, mandibular shape has evolved independently of mandibular size, the direction of shape variation being instead associated with diet and feeding behavior (Monteiro & Nogueira, 2011). Frugivorous bats have similar patterns in loading behavior and pressure point resistance in bones related to the masticatory apparatus that differentiate them between hard-heavy-item consuming species (short and flatten rostrum) and soft-light-item consuming species (elongated and narrow rostrum) (Santana, Grosse & Dumont, 2012).

Ecomorphological studies have demonstrated that morphological variation in bats is majorly associated with trophic specialization, and owing to the fact that bat skulls are under selective pressure to reduce their mass (i.e., reduction of skull mass to meet energetic demands of flight), their morphology might be optimized to meet functional demands (Dumont, 2007). Based on this, our findings might reflect interspecific ecomorphological differences in response to biological specializations for optimizing resource exploitation of soft and light items like Piperaceae fruits, one of the principal components of the diet in Carollia (Nogueira, Peracchi & Monteiro, 2009; York & Billings, 2009).

Evidence of niche differentiation and diet specialization for avoiding ecological competition and niche overlap has been reported in phyllostomid bats (Aguirre et al., 2002; Giannini & Kalko, 2004). Species-specific patterns of intraspecific morphological variation found in our study support the hypothesis of interspecific ecomorphological differentiation, which in Carollia is especially evident in sympatric species, where differences in diet breadth and composition have been studied (Lopez & Vaughan, 2007; York & Billings, 2009).

However, given that recent evidence suggests that more historical processes such as niche conservatism also influence the composition of assemblages in phyllostomid bats (Villalobos, Rangel & Diniz-Filho, 2013), to reach a greater understanding of the mechanisms underlying assemblage composition in this genus, it is advised to combine morphometric and phylogenetic approaches (i.e., community phylogenetics). The latter would give a more comprehensive understanding of the role that both historical, and ecological processes have in shaping the structure of modern geographic patterns of coexistence (Villalobos, Rangel & Diniz-Filho, 2013).

Cranial-mandibular integration and modularity

Cranial-mandibular integration was tested to determine whether the structures work together as a functional unit. Hypotheses tested in this study have been successfully studied in other mammals, revealing the importance of functionality in ecomorphological specialization and differentiation in mammals (Klingenberg & Mebus, 2003; Jojić et al., 2007; Zelditch, Wood & Bonett, 2008; Jojić, Blagojević & Vujošević, 2012). Our results indicate that cranial and mandibular modularity has different, independent patterns. Mandibular modularity was the same for all species, so that patterns in this trait were evident at the genus level, while cranial modularity patterns were species-specific. The lack of variation in mandibular modularity is consistent with findings that modularity patterns in the mandible are genetically patterned, which has been suggested to explain the highly conserved module identity (Klingenberg, Leamy & Cheverud, 2004). The variability found for the cranial patterns may align with evidence that cranial modularity can shift on relatively short time scales in relation to selective pressure (Beldade, Koops & Brakefield, 2002; Monteiro & Nogueira, 2010) and requires further, future investigation in the context of Carollia.

Mandibular modularity has so far not been tested in bats and in this first approach our results agree with those reported previously in other mammals, specifically identifying mandibular modularity as a two-module partition defined by functional traits (Klingenberg & Mebus, 2003; Monteiro & Bonato, 2005; Zelditch, Wood & Bonett, 2008; Jojić, Blagojević & Vujošević, 2012). Presence of these modules (ascending ramus and alveolar region) represents differences in functional specializations between different regions of the jaw for biting and food manipulation (Hiiemae, 2000; Badyaev & Foresman, 2004). The shape of the ascending ramus has evolved to support muscle insertion of masseter, pterygoid and temporal muscles which are related to jaw movement and mastication (Herring et al., 2001). The alveolar region specializes in supporting the dentition and loading capacity, which are important for the masticatory apparatus to resist tension-compression forces applied to the bone structure (Cox, 2008). Finding the same results for all species could indicate that ecomorphological plasticity of the jaw does not affect its modularity, also suggesting that this partition is evolutionarily stable and functionally appropriate for the ecology of these species (Koyabu et al., 2011).

Regarding cranial modularity, differences found in module partitions among species could reflect ecological differences in foraging behavior and niche partitioning and their relation with morphological specializations reported for these three species (Giannini & Kalko, 2004; York & Billings, 2009). These modules (neurocranium and splanchnocranium) represent functional specializations in different areas of the skull; the neurocranium exemplifies morphological specializations for muscle insertion and brain development, and the splachnocranium for biting biomechanics and masticatory activity (Hallgrimsson et al., 2004; Goswami & Polly, 2010; Wellens & Kuijpers-Jagtman, 2013). The latter is reported to be in turn related to morphological diversification in the dentition (Santana, Strait & Dumont, 2011) and rostrum (Nogueira, Peracchi & Monteiro, 2009; Santana, Dumont & Davis, 2010; Santana & Dumont, 2011). Other tested hypotheses that evaluated alternative sources of variation that could explain the presence of these modules in bats (e.g., developmental, genetic or ecological) have been rejected, suggesting a strong correlation between evolutionary conservatism in these modules and its functionality (Goswami, 2007; Santana & Lofgren, 2013). Modifications in the neurocranium are associated with differences between trophic guilds in such a way that cranial structure influences functional importance and recruitment of masseter, pterygoid and temporal muscles during biting (Herring et al., 2001). Additionally, the neurocranium is related to brain development that, in bats, co-varies with foraging behavior and mating systems (Pedersen, 2000; Reep & Bhatnagar, 2000; Pitnick, Jones & Wilkinson, 2006).

Functionality of the masticatory apparatus will depend on the correlation between cranium and mandible structure (Hiiemae, 2000), this correlation was evident from the PLS results, which showed that cranial-mandibular integration explained approximately 80% of the shape variation in all species (Table 3). This integration is due to multiple factors that divide the morphological correlation into regions specialized for muscle insertion (neurocranium and ascending ramus) and regions specialized for biting biomechanics (splachnocranium and alveolar region); these regions together comprise the functional and morphological aspects of trophic diversification and fitness in mammals (Freeman, 1998; Cornette et al., 2013). Morphological integration between the neurocranium and the ascending ramus relates to muscle recruitment, and, depending on the feeding behavior and characteristics of the diet, the functional importance of specific muscles will change, altering the morphology of the skull and jaw in order to work as a functional unit and produce the optimal bite force for each species (Santana, Dumont & Davis, 2010). Consequently, it can be deduced that the morphology of the neurocranium and the ascending ramus will vary jointly, forming a component of a functional unit that will correlate with variation in the rostrum, and that is more important in loading capacity and pressure resistance during biting (Cornette et al., 2013).

Rostrum shape variation in rhinolopid bats has been attributed to evolutionary processes of ecological specialization resulting in niche partitioning among ecomorphologically similar species (Santana, Grosse & Dumont, 2012). These processes respond mainly to functional requirements based on an organism’s alimentary and nutritional needs, which relate to shape diversity for exploiting particular resources (Nogueira, Peracchi & Monteiro, 2009; Labonne et al., 2014). The splachnocranium and alveolar region form the rostrum. These modules correlate functionally with biting biomechanics (Dumont & Herrel, 2003), generating functional convergences in load capacity of pressure points in both the cranial and mandibular structures (Herring et al., 2001; Badyaev & Foresman, 2004).

It is established that in these points of pressure the relationship between the proportional importance of tension-compression forces is the same in the cranium and mandible, integrating the two structures (Herring et al., 2001). Accordingly, cranial-mandibular morphological integration found in this study reveals the presence of a functional unit of the skull and jaw, subdivided into two different modules reflecting the functional requirements for both muscle insertion and biting biomechanics (Santana, Dumont & Davis, 2010; Cornette et al., 2013).

Geographic patterns vs. morphological variation

Our findings may be explained on the basis of two hypotheses that reflect different aspects of the evolutionary history of the genus Carollia. Our results reveal: (1) morphological differences at a phylogroup level for these species, which could be an indicator of ongoing processes of speciation, and (2) geographic patterns of morphological variation in these species are influenced by geographic isolation of populations occurring in the Andes.

For this genus, a phylogroup is defined as a group of individuals that share evolutionary history and a geographic location (Hoffmann & Baker, 2003). In C. brevicauda two different phylogroups have been identified. Both are distributed in Colombia: one covers the Andean, Pacific regions and a portion of the Amazonian region; and the second covers the Caribbean region and a portion of the Orinoquean region. Carollia perspicillata includes three different phylogroups, two of which are present in Colombia, one covering the Pacific and Caribbean regions, whereas the other covers the Andean, Amazonian and Orinoquean regions (Hoffmann & Baker, 2003). In Colombia, only one of the four phylogroups described for C. castanea is present; that phylogroup is present in the Pacific region and it is suggested that another phylogroup could be present in a small portion of the Amazonian region (Pine, 1972; Hoffmann & Baker, 2003).

The patterns found in this study fit with the distribution of these phylogroups, rising the hypothesis that morphological differences between phylogroups can be detected based on the geographical distribution of their morphological variation, further suggesting that our results might shed light on ongoing processes of speciation within C. brevicauda and C. perspicillata (Marchiori & Bartholomei-Santos, 2014). This is supported by the idea that for phyllostomid bats the processes of speciation and diversification in the neotropics are related to the orogeny of the Andes (Hoffmann & Baker, 2003; Velazco & Patterson, 2013). However, our results only give preliminary evidence to this conclusion due to the uncertainty of the exact genetic compatibility between our specimens and proposed phylogroups, so it is suggested for further studies to combine both morphometric and molecular techniques to evaluate this particular hypothesis.

It has been postulated that several species in this genus are different species complexes that remain unsolved (Jarrín, Flores & Salcedo, 2010), so our results could provide insight into this topic. Nevertheless, it will be necessary to perform more detailed studies testing the link between intraspecific morphological differences and the distribution of the phylogroups in the neotropics, in order to detect the presence of undescribed species.

Our second hypothesis focuses on intraspecific ecological differences. Limiting similarity has been described as the main factor that determines the composition of species in the genus Carollia; this contends that species that are more similar ecomorphologically will tend not to coexist thereby avoiding competitive exclusion, and hence more morphologically dissimilar species will coexist (York & Papes, 2007). More recent studies have invalidated this hypothesis, showing that morphologically similar species share environmental space, and that dissimilar species coexist less often (Jarrín & Menendez-Guerrero, 2011). Our results agree with those reported by Jarrín & Menendez-Guerrero (2011) in Ecuador, revealing that C. castanea is the species with the most differentiated ecomorphology and distribution of its morphological variation; conflicting with the limiting similarity hypothesis for this genus in the northern Andes. Jarrín & Menendez-Guerrero (2011) propose that the Andes represent a geographic barrier for C. castanea, isolating populations and generating morphological differences between them. As a conclusion, they stipulate that large body size is a buffer that allows large-bodied species to cross the Andes, maintaining the gene flow and morphological similarities among populations. Our study does not support this. Our results are contrary to those from Ecuador in two ways: (1) we found that for larger species (C. brevicauda and C. perspicillata) not only Andean populations are the only ones morphologically differentiated from other populations across the country, but also populations on opposite versants of the Andes are similar; (2) all C. castanea populations across the country showed the same patterns, such that the Andes do not represent a geographic barrier isolating populations from different regions. Based on our results, we hypothesize that only populations present in the Andes are different in their cranial shape from populations in the rest of Colombia. In this way, the northernmost region of the Andes acted more like an independent and isolated environmental region rather than a barrier splitting lowland areas. Inconsistencies between our results and those reported for bats in Ecuador (Jarrín & Menendez-Guerrero, 2011) may be due to environmental differences between the central and northern Andes. The Andes of Ecuador form one single mountain range, while in Colombia the Andes form three mountain ranges, leading to major ecosystem heterogeneity in the interandean valleys of Colombia (Josse et al., 2009). This could represent a wider range of environments to which species may adapt, occupying greater niche diversity without competition (Bloch, Stevens & Willig, 2011; Pereira & Palmeirim, 2013).

Finally, by comparing results from Ecuador with ours we do not support the hypothesis that large body size favors larger species to cross altitudinal barriers, stabilizing genetic pools and morphologies among populations. Our results elucidate that Andean populations of large-bodied species are morphologically different from populations at lower altitudes, which could be a consequence of gene flow interruption between them. Recently, an inverse relation between body size and altitude was discovered in C. perspicillata, where body size decreases along an altitudinal gradient (De Barros & da Fortes, 2014), supporting our conclusion that large species in this genus do not have a competitive advantage in this regard.