david.jacobs@uct.ac.za (DSJ)

Abstract

The relative contributions of adaptation and drift to morphological diversification of the crania of echolocating mammals was investigated using two horseshoe bat species, Rhinolophus simulator andR. cf. simulator as test cases. We used 3D geometric morphometrics to compare the shapes of skulls of the two lineages collected at various localities in southern Africa. Shape variation was predominantly attributed to selective forces; the between population variance (B ) was not proportional to the within population variance (W ). Modularity was evident in the crania of R. simulator but absent in the crania of R. cf. simulator and the mandibles of both species. The skulls of the two lineages thus appeared to be under different selection pressures, despite the overlap in their distributions. Selection acted mainly on the nasal dome region of R. cf. simulator whereas selection acted more on the cranium and mandibles than on the nasal domes of R. simulator . Probably the relatively higher echolocation frequencies used by R. cf. simulator , the shape of the nasal dome, which acts as a frequency dependent acoustic horn, is more crucial than in R. simulator , allowing maximization of the intensity of the emitted calls and resulting in comparable detection distances. In contrast, selection pressure is probably more pronounced on the mandibles and cranium ofR. simulator to compensate for the loss in bite force because of its elongated rostrum. The predominance of selection probably reflects the stringent association between environment and the optimal functioning of phenotypic characters associated with echolocation and feeding in bats.

Introduction

Understanding the relative contributions of drift and adaptation to organismal diversification is fundamental to studies of evolutionary ecology. To avoid overestimation of selection, drift should always be explicitly accounted for (Bettiet al. , 2010). However, quantifying the relative contributions of these processes to phenotypic diversification is challenging because distinguishing the two processes and identifying their impacts on diversity is difficult (Brandon & Carson, 1996; Millstein, 2002, 2008; Brandon, 2005). Fortunately, there has been some progress in this regard (Millstein, 2008). Adaptation is deterministic and results in phenotypic patterns correlated to environmental/climatic clines (Millstein, 2008). In contrast drift is neutral and results from random processes affecting the genetic composition of populations (Millstein, 2008). In many cases, drift is assumed when evidence for selection is not found (Millstein 2008). However, mathematical approaches e.g. Lande’s model (Lande, 1976, 1979), that allow the quantification of the effects of drift on patterns of phenotypic variation, has made it possible to directly determine the relative importance of drift and selection to phenotypic variation.
Although the application of Lande’s model to phenotypic traits that vary seasonally (e.g. body weight) or are flexible (e.g. behaviour) is theoretically possible e.g. Mutumi et al . (Mutumi et al. , 2017), application of the model to such data might lead to different results depending on when the traits are sampled. In contrast, hard tissue e.g., bony skeletons including skulls provide a more permanent record of the evolutionary processes that a species has endured over its history. Several studies have therefore suggested the use of skulls and geometric morphometrics for enquiries into the relative roles of drift and selection e.g., Evin et al . (Evin et al. , 2008).
Skulls serve functions crucial to the fitness of organisms and their diversification is likely primarily through adaptation (Santana et al. , 2012). The neurosensory system (brain), diet acquisition structures, olfactory system, visual system, speech and sound systems are integrated and housed in the skull. Skulls are therefore subject to diverse selection pressures imposed by the environment on these systems (Cheverud, 1982; Pedersen, 1998; Klingenberg, 2008). For example, the evolution of increased head height, prominent temporal ridge and huge jaw adductor muscles in Chamaeleonid lizards were associated with strong bite force (Herrel & Holanova, 2008). The association between skull morphology and bite force has also been demonstrated in many other vertebrates (Cleuren et al. , 1995; Freeman & Lemen, 2008; Curtiset al. , 2010; Davis et al. , 2010). For example, elongated snouts in some fish appear to be an adaptation which facilitates feeding through suction (Westneat, 2005). Besides dietary adaptations, other behaviours relevant to fitness have shaped the evolution of skull shape. These are grooming (Rosenberger & Strasser, 1985), fighting with conspecifics (Huyghe et al. , 2005), building shelters (Zuriet al. , 1999; Hansell, 2000; Santana & Dumont, 2011) and sensing the environment (Oelschläger, 1990; Ross & Kirk, 2007).
The role of drift was demonstrated in the evolution of human skull form and shape (Rogers Ackermann & Cheverud, 2002; Betti et al. , 2010) using quantitative models. Smith (Smith, 2011) showed that some parts (basicranium, temporal bone, and face) of the skull evolved neutrally whereas the mandible evolved through selection. Quantitative and population genetic methods have shown that isolation between Neanderthal and modern human populations led to cranial diversification through genetic drift rather than the commonly proposed adaptive explanations (Weaver et al. , 2007). Similarly, Ackermann and Cheverud (Rogers Ackermann & Cheverud, 2002; Ackermann & Cheverud, 2004) applied Lande’s model (Lande, 1976, 1979) to variation in the shape and size of human and monkey skulls and found that drift played a significant role. The role of selection may thus be exaggerated if drift is not accounted for quantitatively (Smith, 2011). This is especially important because studies within the same genus have yielded conflicting results. For example, drift was proposed as the cause of phenotypic convergence and divergence in two horseshoe bats, Rhinolophus darlingi (Jacobs et al. , 2013) and Rhinolophus monoceros(Chen et al. , 2009), respectively. In contrast, selection was implicated in the divergence within two other horseshoe bat species,Rhinolophus capensis (Odendaal et al. , 2014), and R. ferrumequinum (Sun et al. , 2013). Thus, two of the four studies on horseshoe bats (genus Rhinolophus) suggest that selection is the predominant driver of diversification but the other two suggest that drift is the main factor. A rigorous test of the processes behind phenotypic diversification should therefore employ models that weigh the relative contributions of adaptation and drift to determine which is the more dominant process shaping phenotypic variation.
The evolution of skull morphology in animals that rely on acoustic signals for communication or navigation (e.g. bats, dolphins, whales, rodents and birds) is particularly interesting because it adds a whole suite of selection pressures on the skull besides those associated with diet and the other five senses (Santana & Lofgren, 2013). For example, there are prominent resonant chambers (forming the nasal dome) in the nasal region of the skulls of horseshoe bats (Rhinolophidae) which acts as an acoustic horn (Hartley & Suthers, 1988; Pedersen, 1998) allowing echolocation call frequencies to be filtered and emitted at high intensity.
Using 3-D geometric morphometrics and Lande’s model we investigated the relative roles of adaptation and drift in two African horseshoe bat lineages, Rhinolophus simulator and R. cf. simulator (Doolet al. , 2016) that are of similar size but differ markedly in the frequency of their echolocation calls. R. cf. simulator was previously classified as R. swinnyi but genetic analyses, using six nuclear markers and an mtDNA fragment (Dool et al. , 2016), indicated that R. swinnyi to the northeast of South Africa was indistinguishable from R. simulator , despite marked differences in the frequency of their echolocation pulses. The frequency of echolocation pulses have a direct impact on the operational range of echolocation and is generally inversely correlated with body size in bats (Jones, 1996, 1999; Jacobs et al. , 2007; Jacobs & Bastian, 2018) and with the volume of the nasal dome in the Rhinolophidae (Jacobset al. , 2014). R. cf. simulator uses higher frequency echolocation calls which are more affected by atmospheric attenuation and probably must emit its calls at greater intensity to achieve the same operational range as R. simulator . We therefore hypothesised that selection rather than drift should be the predominant process in the evolution of skull shape because of the vital sensory and foraging functions of the skull. We predicted: 1) significant deviation from proportionality between the within and between population trait variance in both species (Rogers Ackermann & Cheverud, 2002); 2) modularity should be more prevalent in the crania of both species than in the mandible because of the central role of echolocation to the survival and reproduction of bats. Independence between the cranium and muzzle allows for relatively more flexible response to sensory driven selection. Additionally, the existence of modularity would indicate that the skull is under directional selection because drift and stabilizing selection are inefficient at creating modularity (Melo & Marroig, 2015).

Materials and Methods

Study sites and animals

Skulls were extracted from voucher specimens of both lineages collected in support of two other studies, Mutumi et al . (Mutumi et al. , 2016) and Dool et al . (Dool et al. , 2016). These skulls were supplemented with museum specimens of both lineages (S1 Table). The distributional ranges of the two focal speciesRhinolophus simulator (four localities) and R. cf. simulator (four localities) follow a latitudinal gradient ranging from 16°S to 32°S in south eastern Africa (fig. 1 in Mutumi et al . (Mutumi et al. , 2016)). Both R. simulator and R. cf. simulator lineages have pulses dominated by a constant frequency but at different frequencies with means of 80 and 107 kHz, respectively, when at rest (see fig. S1 in Mutumi et al. (Mutumi et al. , 2016)). The two lineages occur in seven woodland types; eastern half of southern Africa, ranging from DRC in the north, through Zimbabwe and Botswana into South Africa in the south. Woodland types include the Central Zambezian Miombo woodland in DRC and Zambia, the Zambezian and Mopane woodlands, Southern Miombo woodlands, and the Eastern Zimbabwe Montane Forest-grassland Mosaic (Olson et al. , 2001). The southern-most populations occur within Highveld grasslands. In Botswana the sampling site occurred in an ecotone of three woodlands: Kalahari Acacia-Baekiaea, Kalahari Xeric Savannah, and Southern Africa bushveld. Botswana sites experience the driest climate and the Eastern Zimbabwe Montane Forest-grassland Mosaic, the wettest (Olson et al. , 2001).
The specimens were grouped according to the geographic location where they were captured (S1 Fig a and b; Table A1). These locations included north-eastern South Africa (NE), northern Zimbabwe and combined southern Zambia (NZ), Democratic Republic of Congo (DR), south-eastern South Africa (SE), southern Zimbabwe and northern South Africa combined (SZ; S1 Fig a and b; Table A1).
3-D images of each skull were captured through micro-focus X-ray tomography at the South African Nuclear Energy Corporation (NECSA, Pretoria, South Africa; (Hoffman & De Beer, 2012)) following the same procedures as in Jacobs et al . (Jacobs et al. , 2014). All images were imported into the 3-D imaging software, Avizo (version 8.0; Visualization Sciences Working Group, Merignac, France) as volume files. After creating iso-surfaces from the volume files in Avizo, files were saved in ‘Stanford ply’ format and opened in Meshlab (version 1.3.3, Visual Computing Lab of ISTI - CNR, Italy) for placing landmarks. Landmarks were chosen depending on their homology (common and repeatable points on all skulls for each lineage). One skull and one mandible were first tested by repeating the land-marking process 10 times to determine the precision at which landmarks could be placed. Among the tested landmarks the most precise were selected, i.e., 24 landmarks for the skull and 15 for the mandible (S2 Table). Landmarks were placed on only the right half of the skull and the right mandible to control for possible asymmetry (Jacobs et al. , 2014). Each landmark in the 3D space had three coordinates (x, y and z). These sets of three coordinates were used in MorphoJ (version 1.7.0_45; (Klingenberg, 2011)) to analyse shape variation in skulls and mandibles of the two lineages across different localities.
Landmark co-ordinates were analysed as follows. Firstly, a Procrustes superimposition was done on the coordinates to remove variation because of differences in orientation and scale and to standardise the landmarks in a common coordinate system (Adams et al. , 2004). Outliers were checked and extreme cases were double-checked against the original volume files. Where necessary the landmarks were re-inserted on the skull images. A covariance matrix was generated from the Procrustes coordinates, on which a principal components analysis was performed to explore variation in skull shape amongst the different localities for each species. A Procrustes ANOVA (provided in MorphoJ software) was used to test the significance of the differences in skull shapes across localities and between sexes. To visualise the shape differences, a Canonical Variate analysis (CVA) was used. Shape changes in the skulls and mandibles were visualised using the wireframe outlines in MorphoJ which compares shape variations against the average skull shape along each Canonical Variate (CV) with the outlines at the extremes of each CV.
Modularity was also investigated using a priori hypotheses according to Klingenberg (Klingenberg, 2009). Modularity is the differential evolution of different complexes, each complex consisting of groups of traits that evolve together but relatively autonomously from other such complexes (Cheverud, 1996; Wagner, 1996; Klingenberg, 2005). Processes contributing to modularity can be genetic, developmental, functional, or environmental (Klingenberg, 2005). The mandible was divided into subsets of seven (ascending ramus) and eight (alveolar region) landmarks and the cranium was divided into subsets of 10 (basicranium) and eight (rostrum) landmarks as in Jojic et al . (Jojić et al. , 2015). The strength of association between hypothesized modules and all alternative partitions were tested by the CR – covariance ratio in R statistics according to Adams (Adams & Otárola-Castillo, 2013). The CR measures the strength of association between two blocks, i.e., the two modules identified by the covariance matrices of their landmark coordinates compared with the two hypothesised modules (Adams & Otárola-Castillo, 2013). The CR varies from 0 (completely uncorrelated data) to 1.0 (correlated). The strength of the modularity was also measured by the Zcrcoefficient which measures the strength of modularity in each structure – the more negative the coefficient the higher the strength of modularity.

Lande’s Model

The relative contributions of drift and adaptation to the variation in skull and mandible shape was tested by applying the principles of Lande’s model (Lande, 1976, 1979) in the form of the β -test (Rogers Ackermann & Cheverud, 2002), which is described in detail in Mutumi et al . (Mutumi et al. , 2017). The β -test is based on the hypothesis of a log-linear relationship between the variation of phenotypic characteristics between (B) and within (W) populations. If the slope of this relationship is not significantly different from one the null hypothesis is accepted and the observed variations in phenotypic traits can be attributed to neutral evolutionary processes (mutation and drift). Otherwise, the null hypothesis is rejected, which implies that non-neutral evolutionary processes, such as natural selection, can be inferred as the dominant driver of diversification.
Successive landmark coordinates were used to generate Euclidean distances (D ) for successive pairs of landmarks using the following formula:
Where x, y and z are the 3-D landmark coordinates, the subscripts 1 and 2 denote successive positions, and Di is the Euclidean distance for landmark i.
The resulting multivariate response matrix comprisingDi was used to derive the within locality (W ) and between locality (B) variances following the procedure outlined in Mutumi et al . (Mutumi et al. , 2017). Briefly, the Di response matrix was fitted using MANOVA with localities and sex as the categorical predictors to generate a variance/covariance (V/CV) matrix for each species. A measure of the within- population variance W was then obtained in the form of eigenvalues derived from principal component analysis (PCA) on the V/CV matrix. The between population variation B was estimated through multiplication of the matrix of PCA-derived eigenvectors with the matrix of Di means of each locality. We then regressed the log-transformed within variance against the log-transformed between variance and carried out regression t-tests to test the hypothesis that there was no significant difference between the regression slope and one as a function of:
Where β0 is the intercept term and ε is the error (see Mutumi et al . (Mutumi et al. , 2017)).