High Duty Cycle Echolocation May Constrain the Evolution of Diversity within Horseshoe Bats (Family: Rhinolophidae)

Abstract

The phenotype of organisms is the net result of various evolutionary forces acting upon their lineages over time. When an innovative trait arises that confers a substantial advantage in terms of survival and reproduction, the evolution of adaptive complexes between such an innovation and other traits may constrain diversification of that lineage. The specialized echolocation system of the Rhinolophidae may represent such an innovation which affects other parts of the phenotype. We investigated a potential constraint on the diversity of phenotypes of several species of horseshoe bats within a phylogenetic framework. If phenotypic convergence stems from stasis as a result of the specialized echolocation system, phenotypes should converge not only among members of the same species and between sexes but also among species. We analyzed the phenotypic diversity of >800 individuals of 13–16 species. The phenotypes in the horseshoe bats did indeed converge. There was no sexual size dimorphism in mass, forearm length and wingspan within species and there was marked interspecific similarity in both wing and echolocation variables but marked variability in body mass. Furthermore, correlations of wing and echolocation variables with mass suggest that variability within horseshoe bats was largely the result of selection on body size with allometric responses in wing and echolocation parameters, a potential consequence of constraints imposed by their specialized echolocation.

1. Introduction

One of the major questions in biology is why do some clades consist of a single species and others of hundreds of species [1]. Certain general traits (e.g., having a skeleton or living in a terrestrial habitat) have been identified as predisposing certain phyla to diversify into many different species [1]. However, less is known about the causes of diversity within families or genera where all lineages are likely to share traits that favour or retard diversification. However, the number of species can vary markedly from lineage to lineage within orders. For example, the madrone butterfly (Eucheira socialis) is the only species in the monotypic genus Eucheira in a family of butterflies, Pieridae [2], that also includes the genus Delias with about 250 species [3]. Similarly, the bat family Rhinolophidae is a monogeneric (Rhinolophus) family of about 78 species compared to the family Vespertilionidae which includes several genera and more than 300 species [4]. One of these genera, the genus Myotis, has about 100 species [4].

The relatively limited divergence from a common ancestor in some genera may be due to insufficient genetic variation to respond to changing selection pressures, stable environments which exert little selection pressure to change, or the maintenance of existing advantageous traits and trait complexes which may constrain divergence or the evolution of novel ones [5,6]. An extreme example of such a preserved trait is oviparity in birds, the only vertebrate not to have evolved viviparity [7]. Oviparity which was inherited from their reptilian ancestors, may have precluded the evolution of viviparity in birds. The evolution of the intermediate egg retention stage, a possible pre-condition for the evolution of viviparity in birds [7] may have been precluded in birds as a result of some aspect of avian biology, e.g., the potential drop in fecundity that may have accompanied viviparity.

If such constraints are shared amongst lineages, especially as a result of sharing a fitness-enhancing innovation, their phenotypes may exhibit conformity as the conserved innovation prevents diversifying responses to other selection pressures. Traits essential for survival and reproduction such as for example body size, song in birds and frogs or echolocation and wing shape in bats may be particularly constrained by innovations in one or several of these traits preventing novel adaptation in the others [8,9]. For example, innovation in bat echolocation systems might constrain body size and wing shape preventing variation in the latter two traits because these traits are known to form adaptive complexes [9,10,11,12].

Bat echolocation is a unique sonar system used for orientation, foraging and potentially communication [13,14]. There are several kinds of echolocation systems but a particularly specialized kind is that of high-duty-cycle echolocation (HDC) [15]. HDC echolocation is an evolutionary innovation in three bat families, the Hipposideridae, the Rhinonycteridae, and the Rhinolophidae, and in one species of the Mormoopidae (Pteronotus parnellii). HDC allows these insectivorous bats to detect echoes of their emitted pulses reflected off flying insects in situations of high clutter, i.e., echoes from non-target objects such as background vegetation [15]. The specific characteristic of HDC echolocation is that it overcomes the masking effects of clutter through the generation of acoustic glints. Glints are amplitude and frequency peaks generated from the flapping wings of insects [16,17] which are superimposed on the echo from the background. The bats are able to analyse these glints through a combination of an acoustic fovea, a very narrow range of frequencies overrepresented in the auditory pathway, and Doppler shift compensation which keeps the echo frequency within the range of sensitivity of the fovea ([16,17,18]; see Section 2.1), The relatively high frequencies used by HDC bats [19] are suited to the detection and analyses of wing beats of small insects because the glints from small objects are more pronounced the shorter the wavelength (i.e., higher frequency) that impinges and reflects off the object. Overall, this specialized echolocation system allows HDC bats to forage in dense vegetation where distances to background vegetation are short and atmospheric attenuation of the high frequencies used by HDC bats is therefore reduced [20]. Flying in dense vegetation using high frequency pulses with short operational ranges means that these bats have to fly slowly and manoeuvrably. Such flight requires short broad wings that generate lift at low flight speeds [11]. The wings and echolocation system of bats therefore form an adaptive complex by necessity [10,11] which may be constrained by the acoustic and aerodynamic requirements of their specialized ecological niche. If so, this could result in particularly pronounced morphological and echolocation similarity within HDC bats.

There is some evidence that the monogeneric Rhinolophidae are less diverse than bats which use the less specialized low duty cycle echolocation (LDC) [21]. Firstly, both the Rhinolophidae and the Molossidae have highly specialized foraging niches, aerial feeding in confined space and aerial feeding in open spaces, respectively [21]. However, whereas the Rhinolophidae is comprised of a single genus and <80 species [22,23], the Molossidae is comprised of two sub-families, almost 20 genera and more than 100 species [4,24]. Similarly, bats in the genus Myotis, which also does not use HDC echolocation, is comprised of more species (about 100) and has more diverse foraging modes ranging from aerial feeding, gleaning, feeding over water, and trawling [25], despite being a younger genus (approximately 10 Myr compared to 15 Myr for the rhinolophids) [23,26]. The lesser diversity in foraging niche is also reflected in lower phenotypic variation in the rhinolophids compared with the myotids [12]. Furthermore, despite there being 78 species of rhinolophids which display some diversity in body size, their wing shapes are remarkably similar [27]. Although differences in diversity between the Rhinolophidae and LDC echolocators may be the result of adaptive radiation in LDC echolocators, it still begs the question why the younger Myotis, for example, radiated and the older Rhinolophus did not. The specialized echolocation of the Rhinolophus may be at least partly responsible. HDC echolocation as an explanation for convergent phenotypes would be dependent on such convergence being present among functional traits e.g., components of echolocation besides pulse frequency (see Materials and Methods) and wing shape and size. An alternative explanation for convergent phenotypes in HDC bats in the family Rhinolophidae is shared phylogenetic history. Phenotypic analyses have to therefore be phylogenetically informed. Such analyses would also identify cryptic species, an important consideration in a family where species are phenotypically similar enough to be placed in a single genus.

We thus investigated the phenotypes of bat species in southern African representatives of the genus Rhinolophus within the context of potential constraints exerted by evolutionary history and an innovative sensory system in a phylogenetic framework. If convergence stems from stasis as a result of a specialized echolocation system, phenotypes should converge on two levels, within and between species. As outlined above, all horseshoe bat species forage in similar ways using HDC echolocation. Irrespective of phylogenetic history, species may show convergence in traits of their phenotype which are functionally associated with HDC echolocation. We thus compare phenotypic traits among different Rhinolophus species within the context of their phylogenetic relationships. Furthermore, we investigated the within species phenotypic variation between males and females which, again, may show convergence if HDC echolocation constrains the diversification of associated traits. We investigated sex differences in the context of Rensch’s rule [28] which states that sexual size dimorphism (SSD) increases with size in species with male biased SSD and decreases with size in species with female biased SSD. Size differences between the sexes is often attributed to the different roles in reproduction and driven by sexual selection [29] and therefore may contribute to phenotypic differences in bats and might be an opposing force to convergence driven by HDC echolocation.

2. Materials and Methods

2.1. Study Animals

As all rhinolophids use HDC echolocation and thus forage in a similar way, the potential constraints of associated phenotypic traits allow a comparative study among species. The genus Rhinolophus consists of two major phylogenetic clades: an Afro-Palaearctic clade and an Asian clade [23,30,31]. African rhinolophids currently comprise 27 species [32] with many of them having wide geographic distributions [33] over several different biomes. In southern Africa, there are approximately 20 species of rhinolophids and there is therefore at least an opportunity for habitat-mediated phenotypic divergence among African rhinolophids. Rhinolophids are thus ideal for investigating the role of constraints on phenotypic divergence.

In rhinolophid bats, the frequency of the echolocation pulses emitted whilst the bat is at rest on a perch, called the resting frequency, falls within the narrow frequency range of the acoustic fovea. However, the echo of the pulses reflected off of objects returns to the bat’s ears at a higher frequency than the resting frequency when the bat is in flight as a result of the Doppler shift of the returning echo caused by the bats forward motion. This means that the returning echo would fall outside of the frequency range of the bat’s acoustic fovea. The bat compensates for this during flight, hence Doppler Shift Compensation (DSC), by emitting calls at a lower frequency than its resting frequency so that the frequency of the returning echo falls within the range of the acoustic fovea [16,17,18]. The echolocation pulses of rhinolophids are typically of long duration (>30 ms) dominated by a constant-frequency (CF) component with short frequency-modulated (FM) components at the beginning (FMi) and end (FMt) of the CF component [16,17]. These long-duration CF components in combination with DSC allow rhinolophids to receive constant echoes from the background. The bat is then able to detect flying targets by the acoustic glints superimposed on this constant background echo. Acoustic glints are produced when the echolocation pulse is reflected off the fluttering wings of flying insects during the insect’s wing beat cycle [16,17]. These acoustic glints are made possible by the relatively long duration of the echolocation pulses emitted by rhinolophids. When the long CF component of the pulse impinges upon a wing, the amplitude and frequency of the echo depend on the position of the wing and whether the wing is moving towards the source of the pulse or away from it. The amplitude of an echo is positively correlated with the size of the object reflecting the emitted pulse. When the insect wing is perpendicular to the impinging echolocation pulse at the top or bottom of its cycle, all of its surface area reflects the pulse, producing an echo of high amplitude. When the wing is parallel to the direction of the echolocation pulse, only its edge reflects the pulse resulting in a weak echo. Similarly, when the wing is moving down and towards the source of the pulse, during the first part of the power stroke or up and towards the source of the pulse during the recovery stroke, the frequency of the impinging echolocation pulse is Doppler shifted to a higher frequency. In the same way, when the wing is moving away from the origin of the pulse towards the top or bottom of its cycle, the frequency of the echo is Doppler shifted lower than that of the bat’s signal. In this way, acoustic glints, which are amplitude and frequency modulations, are generated and superimposed on the constant echo from the background. All protocols involving animals were approved by the Science Faculty Animal Ethics Committee of the University of Cape Town (Clearance Number 2013/2011/V6/DJ).

2.2. Sampling

Horseshoe bats belonging to 16 species were sampled at several sites (Figure 1) across their geographic ranges in southern Africa. Bats were caught from caves and disused mine shafts using hand nets, harp traps and mist-nets. Harp traps and mist-nets were continuously monitored and captured bats were placed individually in soft cotton bags. Sex and reproductive status were checked immediately following capture so that juveniles, pregnant or lactating bats could be released immediately. Reproductive status was determined by examination of the nipples and palpation of the abdomen of female bats [34]. In addition to this we excluded from analyses all female bats with masses markedly heavier than other females and males of the same species, under the assumption that the heavier females were pregnant. Juvenile bats were distinguished from adults by the presence of cartilaginous epiphyseal plates in their finger bones detected by trans-illuminating the bat’s wings [35]. The forearm length (FA) to the nearest 0.1 mm and body mass to the nearest 0.1 g was measured using dial callipers and a portable electronic balance, respectively. Although body mass is known to vary diurnally and seasonally [36], we included it in the analyses because of its general biological importance and because it was needed to calculate wing loading. Besides excluding pregnant bats, we also minimized variation in mass as a result of gut contents by capturing bats in the afternoon or as they emerged from the roost to forage. This ensured that their guts were empty or close to empty. We however also included FA as a proxy for body size.

Tissue samples for phylogenetic analyses were collected from each captured bat. Skin biopsy punches of 3 mm diameter were taken from the tail membrane using standard skin biopsy punches and stored in 100% ethanol. Liver or muscle tissue were also taken from specimens taken as vouchers and stored in 100% ethanol at −20 °C. All vouchers were collected under research permits issued by the relevant authorities in the respective countries as follows: Botswana (Ministry of Environment, Wildlife and Tourism, EWT 8/36/4 XVI (78)), Malawi (Department of Forestry Licence NO: 1/06/2013/1), Namibia (Ministry of Environment & Tourism 1429/2010), South Africa (Northern Cape Province, Fauna 764/2010; Mupumalanga Tourism & Park Agency, MPB 5253; Cape Nature, 0035-AAA007-00081), and Zimbabwe (Minister of Environment & Natural Resources, 23(1)(C)(11) 25/2011; 19/2012 and 16/2013).

2.3. Phylogenetic Analyses

Whole genomic DNA was extracted with either the DNeasy Blood and Tissue Kit or QIAamp DNA Micro Kit following the manufacturers protocol (Qiagen, Germantown, MD, USA). We used one mitochondrial and two nuclear intron primer pairs as independent markers (mtDNA-R3-F: 5′-TGG CAT GAA AAA TCA CCG TTG T-3′ and mtDNA-F2-R: 5′-ATG GCC CTG AAG AAA GAA CCA GAT G-3′, ACOX2-F1: 5′-CCT SGG CTC DGA GGA GCA GAT-3′ and ACOX2-R1: 5′-GGG CTG TGH AYC ACA AAC TCC T-3′ and COPS7A-F1: 5′-TAC AGC ATY GGR CGR GAC ATC CA-3′ and COPS7A-R1: 5′-TCA CYT GCT CCT CRA TGC CKG ACA-3′) [37,38,39]. The marker mtDNA amplifies a region on the mitochondrial genome coding for cytochrome b, tRNA Treonine, tRNA Proline and a part of the control region [37]. The two intron fragments are non-coding parts of the nuclear genes: “ACOX2” coding for acyl-CoA oxidase 2 and “COPS7A” coding for the signalosome complex subunit 7a respectively [38,39]. These markers were chosen based on the findings by Dool et al. [23] who showed that these markers provide a robust phylogeny of the family Rhinolophidae. Polymerase chain reactions took place in a 25 μL reaction volume. The mastermix (Kapa Biosystems, Wilmington, NC, USA) contained 18.05 μL double distilled water, 2.50 μL buffer with loading dye (10×, conc. 1×), 0.75 μL MgCl (50 mM, conc. 1.5 mM) or 0.50 μL MgCl for the ACOX primer mastermix, 0.50 μL dNTPs (10 mM, conc. 0.2 mM), 0.50 μL of each primer (20 µM, conc. 0.4 μM), 0.20 μL Taq polymerase (5 Unit/µL, conc. 1 U), and 2.00 μL DNA (conc. ~10–80 ng/μL). The touchdown-PCR conditions for all three markers were 95 °C for 10 min; 2 cycles of 95 °C for 15 s, 65 °C for 30 s, 72 °C for 1 min; followed by 2 cycles each at annealing temperature in 2 °C decrements from 65 °C (63 °C–57 °C); 30 cycles of 95 °C for 15 s, 55 °C for 30 s, 72 °C for 1 min; 72 °C for 5 min [23]. Successful marker amplification and general quality of the PCR was checked using a DNA ladder on 1% Ethidium Bromide agarose gels under UV light. Prior to Sanger sequencing (Macrogen Inc., Seoul, South Korea), PCR products were purified using an Exo-FastAP enzymatic reaction (0.125 μL Exo1 20U, 1.0 μL FastAP 1 U, 0.5 μL dH2O; 1.5 µL EXO-FastAP mix (Thermo Scientific, Waltham, MA, USA) were premixed and added to 5.0 µL of the PCR product. The reaction was done at 37 °C for 15 min followed by 85 °C for 15 min and 5 min at 4 °C. Forward and reverse sequences were edited manually in ChromasPro (v1.5, Technelysium PTY Ltd., Brisbane, Australia) and aligned in Mega6 (v6.06; [40]). Previously published alignments of the markers [23] were used as the basis but expanded with species relevant for this study. Alignments were optimised using MAFFT [41] and subsequently analysed for poorly aligned parts using gblocks (v0.91b; [42]).

To reconstruct the phylogenies, the appropriate substitution model for each of the three markers was calculated in jModeltest [43]. The parameters of the model with the best AIC support and lowest ln–Likelihood score, were used in the subsequent analysis. The phylogenetic reconstruction of 16 species based on 451 individuals for the three markers was carried out on the web-based CIPRES platform (v3.3; [44]). Bayesian inference (MrBayes v3.2.3; [45]) was conducted as two independent runs each with 5–10 million generations depending on the result of the convergence diagnostics i.e., the average standard deviation of split frequencies, the potential scale reduction factors as well as ESS values (v1.6.0; [46]) and sampled every 500 generations. The independent runs were combined with Logcombiner (v1.8.4; [47]) thereby discarding values prior to convergence of the runs as burnin assessed via the plot of the log-likelihood through time (Tracer v1.6.0; [46]). A maximum credibility tree including branch lengths and posterior probabilities was generated (TreeAnnotator v1.8.4; [47]) and visualized with FigTree (v1.4.3; [48]).

The phylogeny was then pruned to include only those taxa relevant to our analyses. As R. simulator and R. cf. simulator are phenotypically distinct due to non-overlapping resting frequencies (±80 and ±104 kHz, respectively) but form genetically one clade [23], we had to split the R. simulator the clade by adding very short branch lengths near 0 using the “multi2di” command in the “ape” package [49] to obtain separate OTUs (tips of the tree) for R. simulator and R. cf. simulator in order to match the datasets of phenotypic and phylogenetic datasets for the phylogenetically linear regressions (PGLS). To further match the two datasets, we chose individuals for sequencing from which we collected the phenotypic measurements or if that was not possible we chose specimens which were collected in the same geographic area as the specimens in the phenotypic analysis. The resulting tree topology served as the input phylogeny for the PGLS of phenotypic data.

2.4. Wing Morphology

We took photographs of the right wing of each bat extended in accordance with [50] using a digital camera (Fuji Finepix S1PRO.E, Fuji Photo Film Co. Ltd., Tokyo, Japan) positioned directly above the wing and parallel to the flat table top on which the wing was extended. This minimized angular distortion so that lengths and areas of the wing components could be measured using SigmaScan Pro 5 version 3.20 (SPSS Inc., Cary, NC, USA). The wing was extended on graph paper to allow calibration of the software programme. The size and shape of the wings were assessed through the standard wing measurements of wingspan, wing area, aspect, ratio, wing loading, wingtip length ratio, wingtip area ratio and wingtip shape index, calculated in accordance with Norberg and Rayner [11].

2.5. Echolocation

Echolocation pulses from bats were recorded (for approximately 45 s per bat) from hand-held individuals 30 cm from the microphone of either an Avisoft UltraSoundGate 416 (Avisoft Bioacoustics, Berlin, Germany), using a condenser ultrasound microphone (Avisoft-Bioacoustics CM16/CMPA) or an ultrasound D1000X detector (Pettersson Elektronik AB, Uppsala, Sweden) at sampling frequencies of 500 kHz or 384 kHz. Hand-held calls allow the determination of resting frequency (RF; frequency of maximal energy when at rest) in rhinolophid bats [51] while eliminating variations in frequency as a result of the bats compensating for Doppler shifts during flight [52].

In addition to RF, measured at the centre of the constant frequency component of the pulse, we also measured several other pulse parameters (see [53] and Results Section in this study) from 5 to 30 pulses of hand-held recordings with the best signal to noise ratios for several individuals of each species. Pulses were not chosen from the first 10 pulses of a recording because horseshoe bats are known to tune into their resting frequencies from lower frequencies after periods of silence [51]. The chosen pulses were analysed using Avisoft SASLab Pro automatic measurement function (Avisoft Bioacoustics, Version 4.2, Glienicke, Germany). A Hanning window was used to eliminate the effects of background noise. Measurements were averaged over the pulses of each individual and these average measurements were used in statistical analyses.

2.6. Functional Association between Echolocation and Wing Morphology

In bat guilds, echolocation and wing shape are thought to show adaptations to the foraging mode and habitat [21]. Rhinolophids belong to the guild of narrow space clutter foragers which means they hunt for flying insects close to vegetation. As laid out in the introduction, the HDC echolocation system has evolved to detect and analyse the wing beats of insect flying close to vegetation (for a review, see [54]). This foraging situation requires manoeuvrable flight with narrow turns at low flight speed which is facilitated by shorter and broader wings [11]. We assessed the association between the two functional traits with phylogenetically informed regressions (PGLS; see Section 2.8) across 13 Rhinolophus species. This informed us about variation between species and potential ecological niche separation within this guild [55].

2.7. Sexual Dimorphism

We assessed sexual dimorphism using acoustic parameters of echolocation calls, and indicators of body size mass, forearm length and wingspan. We calculated size dimorphism for each species as the mean measurement for males (M) divided by the mean measurement for females (F; [56]) while controlling for phylogeny through a phylogenetic generalized least squares (PGLS; see Section 2.8) to test for a relationship between ln(M/F) and lnF [56]. Rensch’s rule predicts that in species in which females are bigger than males the slope of this relationship should differ from that for isometry (i.e., 1) and should decrease with size in species with female biased sex ratios (i.e., F > M) and increase with size in species with a male-biased sex ratio (i.e., M > F; [56,57]).

2.8. Statistical Analyses

All phenotypic analyses were done in Statistica (v13.3, TIBCO Software Inc., Palo Alto, CA, USA) and SPSS (v25, IBM Corp., Armonk, NY, USA) unless otherwise stated and all data were log transformed unless otherwise indicated. We used a global level of significance of 5% but report the Bonferroni-adjusted α-level (as α^adj) to allow the reader to assess the effect of multiple comparisons and Type I error. We checked for sexual dimorphism by doing MANOVAS with species and sex as categorical factors and size and wing measurements or echolocation measurements as dependent variables. We kept the numbers of individuals of the two sexes equal as far as possible.

We analysed interspecific variation in size and wing variables as well as echolocation variables through discriminant function analyses (DFA) with 10,000 bootstrap iterations on principal components (PCs) extracted from a principal components analyses (PCA) on the two sets of data. PCAs were conducted to extract uncorrelated principal components in accordance with the assumptions of the DFA. Only PCs with eigenvalues >1 [58] were used in the DFAs.

Closely related species such as the Rhinolophus species studied here share a common ancestry. One can therefore assume that phenotypic traits across species may show similarities as a result of shared ancestry. In studies investigating potential functional adaptations of traits across different species, e.g., wing morphology and its adaptation to different foraging environments, it is necessary to incorporate knowledge about the phylogenetic relationships [59]. Phylogenetic generalised least squares (PGLS) provides a method which takes into account how traits covary due to shared ancestry [60]. We computed phylogenetically informed regressions in R (v3.5.0; [61]) using RStudio (v1.1.453; [62]) and the packages “caper” (v1.0.1; [63]) with “nlme” (v3.1) and “ape” (v5.1; [49]). The PGLS models for the seven regressions were fitted with maximum likelihood estimations of Pagel’s lambda λ [64,65] assuming Brownian motion model of trait evolution (modelfitting was done in the “geiger” package, v2.0.6; [66]). A lambda of 0 indicates no correlation due to phylogenetic relationship i.e., no phylogenetic signal and λ = 1 is consistent with a strong correlation of the traits due to phylogeny and a Brownian Motion model of trait evolution [67].

3. Results

3.1. Phylogeny

The tailored phylogeny is based on a supertree constructed from phylogenies of each of the three markers. The supertree approach was chosen because the individuals sequenced did not completely overlap between the marker sets. The phylogenetic reconstruction of the mitochondrial marker (927-bp alignment) was done under the GTR+I+G model (AIC score: 30,022.621760; −lnL score: 14,102.3109), the intron COPS (621-bp alignment) under the GTR+G model (AIC score: 21,049.566600; −lnL score: 9918.7833) and the intron ACOX (390-bp alignment) under the TVMef+G model (AIC score: 4473.692040; −lnL score: 1780.8460). To obtain a phylogeny matching the species in the phenotypic datasets, we collapsed individual tips within a species-clade (TreeGraph2 v2.14.0-771 beta; [68]). Compared to the latest published Rhinolophus phylogeny [23], we obtained the same topology with good posterior probability support values (Figure 2) but expanded the taxonomic coverage by including R. eloquens and R. deckenii. We subsequently excluded R. hipposideros, R. mehelyi, R euryale, R. alcyone, R. hildebrandti and R. blasii (Figure 2) for the PGLS regressions reported here. As in Dool et al. [25], there are at least two cryptic lineages. R. fumigatus from the western part of southern Africa appears to be a distinct but sister lineage to R. fumigatus from the eastern half of Africa ([23]; this study). We thus refer to them as R. fumigatus_West and R. fumigatus_East, respectively. We used the designation of R. cf. simulator for a lineage that displayed genetic similarity to R. simulator [25] but echolocated at a much higher frequency (see Results Section and [23]).

3.2. Functional Association between Body Size and Wing Morphology

We collected wing and body size data for 13 species and a total of 842 individuals (Table S1). There were significant differences among species (Factorial Anova: F_{108, 5895} = 68.0, p < 0.0001 *, α^adj = 0.0042) and sex (F_{9, 808} = 7.0, p < 0.0001 *, α^adj = 0.0039).

We calculated SSD for mass and forearm for a total of 925 individuals and 13 species and 871 individuals for wingspan (Table S2). Despite significant differences in the sexes across all species (see Anova results above), SSD was not evident within the species of southern African horseshoe bats. Within species, the sex ratios were either equivalent to one, i.e., no gender bias, or only slightly female biased (SSD just below one; Table S2). Phylogenetically informed regressions of mass, forearm length and wingspan against the male:female ratios of these variables were all insignificant (

Table 1. Results of the PGLS regressions comparing echolocation and wing variables between species and between genders (m, male; f, female). An asterisk (*) marks a significant result at the 5% level. The Bonferroni corrected adjusted level of significance is α^adj = 0.0071.

Species Comparison	lnWL × lnMass	logRF × logMass	logRF × logDur	ResidRF × ResidWL
Lambda (ML)	0.00	0.00	0.56	0.00
Intercept b (SE)	0.98 (0.23)	2.41 (0.14)	2.10 (0.41)	6.86 × 10−10 (2.39 × 10−2)
t-test	4.30	17.52	5.18	0.00
p-value	0.001 *	2.2 × 10−9 *	0.000 *	1.00
Slope β (SE)	0.41 (0.11)	−0.51 (0.13)	−0.13 (0.26)	−9.41 × 10−2 (1.79 × 10−1)
t-test	4.18	−3.76	−0.51	−0.53
p-value	0.002 *	0.003 *	0.62	0.61
Adjusted R2	0.58	0.52	−0.07	0.064
F-ratio	17.47	14.16	0.26	0.28
p-value (df 1, 11)	0.002 *	0.003 *	0.62	0.61
Sexual Dimorphism	lnMassf × lnMassmf	lnFAf × lnFAmf	lnWSf × lnWSmf
Lambda (ML)	0.66	0.57	0.58
Intercept b (SE)	2.29 (0.25)	1.47 (0.07)	3.45 (0.06)
t-test	9.28	21.14	61.06
p-value	1.561 × 10−6 *	2.951 × 10−10 *	2.665 × 10−15 *
Slope β (SE)	−1.12 (2.05)	−5.17 (2.25)	−0.72 (0.80)
t-test	−0.60	−2.30	−0.90
p-value	0.56	0.04	0.39
Adjusted R2	−0.06	0.26	−0.02
F-ratio	0.36	5.28	0.82
p-value (df 1, 11)	0.56	0.042 *	0.39

).

The DFA confirmed the results of the MANOVA that there was significant variation in body size and wing morphology amongst species (DFA F_{36, 2444} = 137.4, p < 0.000; Figure 3). The first three principal components (PC) with eigenvalues >1 recovered from the PCA explained 80.2% of the variance. DFA on these three PCs extracted two roots which explained 99.5% of the variance. Root 1 which explained 94% of the variation, was associated with size and manoeuvrability variables including mass, forearm length, wing area, wing span and wing loading (i.e., PC 1,

Table 2. Results of the PCA and DFA on body size and wing variables. Only factors with eigenvalues >1 and only factor loadings >0.7 are shown.

Variable	PCA Factor Loadings
	PC 1	PC 2	PC 3
Mass	−0.97
FA	−0.97
Wing area	−0.91
Wingspan (cm)	−0.90
Wing loading (Nm-2)	−0.79
Wingtip Shape Index		−0.89
Wingtip area ratio		−0.85
Aspect ratio			−0.83
Eigenvalues	4.22	1.86	1.14
Cumulative %	46.84	67.48	80.17
Variable	DFA Results
	Root 1	Root 2	Wilks’ λ	F to Remove	p
PC 1	0.96	0.23	0.08	839.91	<0.0001
PC 3	−0.03	0.63	0.74	21.40	<0.0001
PC 2	0.01	−0.60	0.76	24.07	<0.0001
Eigenvalue	13.10	0.75
Cumulative %	94.1	99.5
Wilks’ λ	0.04	0.54
χ2	2723.80	519.37
df	36	22
p	<0.0001	<0.0001

and Figure 3). Root 2 which explained a further 5% of the variation was associated with shape variables, namely, aspect ratio (PC 3) and the tip shape index (PC 2, Table 2). Species, arranged from the largest to smallest from left to right in Figure 3, formed three major clusters. The cluster on the left consisted of R. fumigatus_West, R. eloquens and R. clivosus which are the largest species with the highest wing loadings. The middle cluster was composed of R. deckenii, R. fumigatus_East, R. capensis and R. damarensis, which have intermediate body sizes and wing loadings (Table S1). The cluster on the right included R. landeri, R. swinnyi, R. denti, R. simulator and R. cf. simulator. These were the smallest species with the lowest wing loadings. Rhinolophus darlingi formed a group on its own, loading between the middle and right clusters along Root 1 and higher than the other species along Root 2. This is probably because it had one of the lowest aspect ratios and the highest tip shape index, suggesting that its separation from the other species along Root 2 was largely based on it having a shorter and broader wing with a more rounded wing tip (Table S1). Even within clusters, the squared Mahalanobis between species pairs were all significant (F_{3, 827} > 4.6, p < 0.003). However, the total classification success was 53% and was under 50% for R. capensis and R. damarensis, in the left cluster and for R. denti, R. simulator in the right cluster of the graph in Figure 2. Two of the species in the left cluster, R. capensis and R. damarensis, are morphologically so similar so that individuals of both species were statistically misclassified and often grouped into either species category. The low classification success among species in the right cluster, R. denti, R. simulator, however is based on morphological similarities to other species of this cluster namely R. cf. simulator for R. denti and R. swinnyi for R. simulator despite marked differences in the resting frequencies of the latter pair (Table S1).

3.3. Functional Association between Body Size and Echolocation among Species

We collected echolocation data for 16 species and a total of 676 individuals (Table S3). There were significant differences among species (Factorial Anova: F_{180, 5996} = 43.1, p < 0.0001) and sex (F_{12, 633} = 1.8, p < 0.05). Sex differences did not vary significantly among species (species*sex: F_{180, 5996} = 1.1, p > 0.28). Univariate results indicated that there were gender differences across species in RF only (F = 11, p < 0.001) and the Tukey post-hoc tests suggests that this difference was not between genders of the same species (p > 0.74, α^adj = 0.0031).

There was significant variation in echolocation amongst species (DFA F_{60, 2566} = 140.0, p < 0.0001; Figure 4,

Table 3. Results of the PCA and DFA on echolocation variables. Only factors with eigenvalues >1 and only factor loadings >0.6 are shown.

Variable	PCA Factor Loadings
	PC 1	PC 2	PC 3	PC 4
Inter-onset interval	0.85
Inter-pulse interval	0.78
Duration of pulse	0.75
Sweep rate of FMt	−0.74
Minimum frequency of FMt		0.91
Resting frequency		0.90
Minimum frequency of FMi		0.86
Duration of FMt			−0.73
Sweep rate of FMi			0.68
Duty cycle				−0.73
Eigenvalues	4.08	3.12	1.97	1.49
Cumulative %	34.01	59.98	76.42	88.82
Variable	DFA Results
	Root 1	Root 2	Wilks’ λ	F to Remove	p
PC 2	−0.33	0.59	0.14	282.60	<0.0001
PC 1	0.07	0.48	0.48	48.50	<0.0001
PC 3	0.05	0.58	0.49	46.53	<0.0001
PC 4	0.01	−0.37	0.72	16.83	<0.0001
Eigenvalue	54.87	1.38
Cumulative %	95.82	98.23
Wilks’ λ	0.003	0.193
χ2	3770.20	1094.90
df	60	42
p	<0.0001	<0.0001

). The first four principal components (PCs) with eigenvalues >1 recovered from the PCA explained 88.8% of the variance. DFA on these four PCs extracted two roots which explained 98.2% of the variance. Root 1, which explained 95.8% of the variation, was dominated by RF and the minimum frequencies of FMt and FMi (i.e., PC 2) and secondarily by the temporal components of the echolocation pulses including inter-onset interval, inter-pulse interval, duration of the whole pulse (PC 1), duration of FMt and sweep rate of FMi (i.e., PC 3, Table 3 and Figure 4). Root 2 which explained a further 2.4% of the variation almost equally correlated with the four PCs including duty cycle (i.e., PC 4, Table 3) which did not contribute to the discriminatory power of Root 1. Species, were arranged from the highest echolocators on the left to the lowest echolocators on the right with three main clusters. Species echolocating at higher than 100 kHz (R. cf. simulator, R denti, R. landeri and R. swinnyi) formed the first cluster on the left. The second cluster was made up of R. clivosus and R. blasii, with mean RFs of 91 kHz, and the third cluster consisted of several species with mean RFs frequencies ranging 80–85 kHz and included R damarensis, R. capensis, R. darlingi, R. simulator and R. deckenii which ranged in mean resting frequencies 74–92 kHz (Table S1). The third cluster on the right was more diverse in their resting frequencies and included R. fumigatus_East and _West echolocating at 55 kHz. Rhinolophus hildebrandtii and R. eloquens formed a pair of species echolocating at 44 kHz with the species on the far right, R. mossambicus, echolocating at 38 kHz. The species are arranged from left to right in increasing order of body size (Figure 4) as in Figure 3. There was much overlap along Root 2.

As with morphological variables, the squared Mahalanobis between species pairs were all significant (F_{4, 657} > 4.2, p < 0.003) with the exception of that between R. fumigatus_East and _West (F = 0.5; p > 0.7). The total classification success was slightly higher than that for morphological variables at 68.6% but was less than 50% for R. capensis because of its overlap with the acoustically similar congenerics R. damarensis and R. blasii (Figure 4 and Table S1).

3.4. Phylogenetically Informed Analyses of Associations between Body Size, Wing Morphology and Echolocation among Species

Regressions involving forearm length did not differ qualitatively from those using mass as a measure of body size, and are thus not reported. Wing loading and resting frequency were positively and negatively associated with mass, respectively (Table 1 and Figure 5 and Figure 6). There was no correlation between resting frequency and call duration nor between resting frequency and wing loading after controlling for body size and phylogeny (Table 1).

4. Discussion and Conclusions

The phenotypes in the Rhinolophidae did indeed converge both between gender and among species. Although sexes differed in size and wing variables this difference was restricted to some species and not large enough to manifest SSD in mass, forearm length and wingspan within species. Similarly, sex differences among echolocation variables were restricted to RF and only among species but not within species. Although there were significant differences among species in size, wing and echolocation variables (Anova results), the correlation between mass and wing and echolocation variables suggest that these differences were largely the result of selection on body size with allometric responses in wing and echolocation parameters. This was supported by the absence of a correlation between wing loading and resting frequency after removal of the influence of body size and phylogeny (Table 1).

There is strong evidence for adaptive complexes between traits associated with body size, wings and echolocation in southern African rhinolophids. Firstly, DFA indicated that with only two exceptions the major difference between southern African rhinolophids was in body size with very little divergence along axes associated with wing shape. However, despite the divergence in body size, these species did form clusters, which is an indication of convergence. Furthermore, the strong allometric relationships between mass and both wing and echolocation variables (Table 1 and Figure 5 and Figure 6) are indicative of adaptive complexes.

Variables related to body size were most responsible for the divergence among species of southern African Rhinolophidae whether size was based on mass or forearm length. Although sexual size dimorphism is prevalent in other mammalian taxa in accordance with Rensch’s rule (as cited in [57]), this was not evident in southern African rhinolophids. Sexual size dimorphism is usually driven by sexual selection on males, favoring larger males as a result of increased fecundity. In bats, however, the few studies done on SSD showed more ambivalent results with little support for a general SSD. SSD in bats is more often than not female biased [69] and this was the case in southern African rhinolophids where >60% of species had slightly female-biased SSDs (this study) and global Rhinolophidae where SSD was largely female biased [70]. However, in both studies on the Rhinolophidae, female body size was only slightly larger than male body size and Rensch’s rule was not supported after controlling for phylogeny. Similarly, in four species of Phyllostomidae, another bat family, SSD was evident among three of the four species across localities but only significant in one species (Anoura caudifer; [71]). Although sexual selection upon which SSD is based, may be evident in bat populations, particularly on traits like echolocation frequency which has been implicated in mate choice [72], the traits on which the general patterns of Rensch’s rule are usually based (e.g., body size) may be subject to only weak selection [70]. Furthermore, natural populations may be subjected to several types of alternative evolutionary processes including drift, habitat-mediated selection and selection to maintain adaptive complexes. All of these may be severe enough to constrain the evolution of sexual dimorphism. Lastly, sexual size dimorphism may not always be evident in traits that may nevertheless be subject to sexual selection [73] because mate choice based on echolocation frequency [72] may result in directional selection in echolocation frequency without sexual dimorphism in bats, especially if changes in echolocation frequency are heritable and expressed in both sexes.

Body size is strongly correlated with many physiological and fitness characters [74,75] and there are a number of factors that could be responsible for its divergence. James’ rule [76], for example, proposes that in cool, dry areas larger body sizes are generally favoured by selection, whereas smaller body sizes are favoured in cool, humid areas and that hot, humid areas would select for smallest body sizes. Such climatic variability associated with altitude has been proposed for the evolution of large body size in a southern African rhinolophid species complex in high altitude refugia and small body size at lower altitude refugia during the Plio-Pleistocene. This divergence in body size resulted in an allometric response in echolocation leading to divergence in pulse frequency [77]. Similarly, in an earlier phylogenetic analyses of echolocation divergence across more than 30 species of Rhinolophidae, it was found that such divergence is best explained by changes in body size with an allometric response in echolocation [30] providing strong support for an adaptive complex between body size and echolocation.

Furthermore, although body size is divergent among species, the allometric relationship between body size and wing variables suggests that southern African rhinolophids are simply scaled up or scaled down versions of the same phenotype. Phenotypic functionality with regards to flight appears therefore to have been maintained across species and across biomes. If so, it may be indicative of the specialization of this family to hunting in vegetation. As far as is known, all species in this family hunt close to or within vegetation [12,78] and their wing morphology is thus specialized for the slow, manoeuvrable flight needed to negotiate confined spaces within clutter [11]. Rhinolophids have longer wingspans and greater wing areas for their mass than most other families and consequently one of the lowest wing loadings too. They also have more rounded wing tips (i.e., high tip shape index (I)). These characteristics allow them to maintain lift at low flight speeds and consequently capable of the slow, manoeuvrable flight needed for flying in confined spaces [11]. The strong allometric relationship between body size and wing loading (all but two of the species are within 95% confidence intervals; Figure 5) suggests strong selection pressure to maintain the allometric relationship across species despite the diverse biomes that each of these species occupy. This may be indicative of the evolution of specialized mode of flight as a result of an adaptive complex between locomotion and a specialized echolocation system.

There is marked variation in echolocation variables among species and resting frequencies is largely responsible for this ranging from a mean of 38 kHz in R. mossambicus to 107 kHz in R. denti (Figure 4 and Table S3). However, as with the wing variables, convergence is evidenced by the formation of several clusters, species of the same size tend to echolocate at more or less the same frequency, larger species echolocating at relatively lower frequencies. This is supported by the strong inverse allometric relationship between body size and resting frequency (Figure 6). Furthermore, the allometric relationships between body size and wing and echolocation variables is suggestive of an adaptive complex between body size and variables associated with flight and echolocation.

The inverse relationship between body size and echolocation frequency reported here for rhinolophids and for several other families of bats [78,79,80,81] can be attributed to the physics of sound production—smaller bats have thinner/shorter vocal chords and smaller resonant chambers producing higher frequency sounds.

Sound radiates from the source and travels through the atmosphere as concentric pressure waves. The generation of these waves dissipates the energy not only due to absorption of the sound energy by air molecules but also by spreading losses as the energy of the sound is spread over an increasingly large wave front. This dissipation of the energy in the sound is called atmospheric attenuation. Atmospheric attenuation is largely dependent on temperature and humidity of the atmosphere and the frequency of the sound being transmitted. All else being equal, higher frequency sound is attenuated more than lower frequency sound [20,82,83] and the operational range of high frequency echolocation at similar intensities is generally lower. This means that bats using short range, high frequency echolocation have to be small with low wing loading to allow the slow manoeuvrable flight necessary to react to a target detected over short distances. Larger bats on the other hand are generally limited to flying in larger spaces because of their size but also because wing loading tends to increase with size (this study; [11]) and they are therefore less manoeuvrable. This means that the operational range of their echolocation has to be longer so that they may detect targets in sufficient time to allow them space within which to manoeuvre. Hence, the correlation between body size and echolocation. Furthermore, it has been proposed that echolocation has constrained the evolution of large body sizes in echolocating bats because, first, echolocation is a short range detection system and bats have to be small to be able to manoeuvre [80] and, second, bats couple wing beat rate with emission of echolocation pulses to save energy [81,84]. Coupling may constrain the evolution of large body size in echolocating bats through limitation of pulse repetition rate because very large bats may not be able to emit pulses at a high enough rate to meet the energy demands of a large body size [81,84]. Larger bats would need to emit at a high rate to get detailed enough sonar images to successfully hunt sufficient prey to sustain their larger bodies. If emission is coupled with wing beat, these large bats would have high wing beat rates which: (1) is energetically costly; and (2) would further impact negatively on their manoeuvrability. Thus, the echolocation system of HDC bats, which is specialized for clutter detection and the emission of several pulses per wing beat, may have constrained their size (they are generally smaller than other echolocating families [80]) and, given the allometric relationship with wing variables, the diversity of their wing morphology as well. This might explain the phenotypic similarity of the Rhinolophidae. Furthermore, the adaptive complex between body size, wing morphology and echolocation was not evident in the genus Myotis, from a different bat family, which uses the less specialized low duty cycle echolocation (LDC) [12].

However, there are some species which deviate from the allometric relationships for both wing and echolocation variables (Figure 5). Nevertheless, these exceptions support the view that similarity in the Rhinolophidae is due to factors other than a lack of genetic variation preventing responses to selection pressures in the diverse habitats occupied by these species. Rhinolophus darlingi and, to a lesser extent, R. fumigatus_West are separated from the rest of the species along Root 2 in Figure 3. Both species fall outside of the 95% confidence intervals of the allometric relationship between body size and wing loading (Figure 5). These deviations are the result of a more rounded wing tip (higher wing tip shape index; Table 1) in R. darlingi and a slightly more pointed wing tip in R. fumigatus_West. In addition, R. darlingi has a higher wing loading and R. fumigatus_West a lower wing loading than expected from their body masses (Figure 5). Rhinolophus darlingi thus couples a rounded wing tip with a relatively high wing loading for its body mass with a resting frequency that scales allometrically with body mass (Figure 6). Similarly, the wing loading of R. fumigatus_West deviates from allometry (Figure 5) but it has a resting frequency that scales allometrically with body size (Figure 6). In contrast, R. clivosus, R. eloquens, R. fumigatus_East and, to a lesser extent, R. simulator have wing loadings which scale allometrically (Figure 5) but resting frequencies which deviate from allometry. R. clivosus has a much higher resting frequency than expected from its mass (see [12]) and its wing loading (Figure 5) whereas the other three species have lower resting frequencies than expected from allometry (Figure 6). These deviations of wing and echolocation variables may involve differential selection pressure on wings and echolocation such that the adaptive complex normally present among these variables, is disrupted. In any case these exceptions suggest that Rhinolophidae have sufficient genetic variability to respond to selection pressure and that their phenotypic similarity is probably due to the constraints of their specialized echolocation system.

It is also unlikely that phenotypic similarity is the result of shared ancestry because the allometric relationships held after controlling for phylogeny (Figure 5 and Figure 6) and similarity in phenotypic characteristics occurred across clades (Figure 2, Figure 3 and Figure 4) and has been reported amongst southern African rhinolophids that are not genetically closely related [8,27]. The wide divergence from the ancestral body size and echolocation frequency (see [27]) and the wide range of body size within extant species further supports the existence of sufficient variation in the Rhinolophidae to respond to other evolutionary forces in their environment. However, with few exceptions (this study; [12]) this response is generally constrained within allometric relationships. This is also true for another family of HDC bats, the Hipposideridae, in which allometry is maintained despite wide ranges in body sizes. In contrast, the allometry between body size and echolocation is absent in at least one genus (Myotis) of LDC bats [12]. It is also absent in the diverse LDC family, Phyllostomidae but present in the Vespertilionidae and the Molossidae.

HDC, Doppler-shift compensated echolocation does relax the constraint on the duration of echolocation pulses, hence the high variation among species and the absence of a correlation with resting frequency. LDC bats prevent damaging their sensitive inner ears during the emission of high intensity echolocation pulses by reducing the transmission of sound to the inner ears. This is done by disrupting the link between the ossicles in the middle ear via contraction of the muscles attached to these tiny bones. LDC echolocators thus have to separate pulse and echo in time to avoid receiving the echo when the transmission of sound to the middle ear is thus disrupted. However, as frequencies increase, detection ranges decrease as a result of the increased atmospheric attenuation of higher frequency sound. LDC bats have to reduce the duration of pulses to avoid pulse/echo overlap over these shorter detection distances. There is therefore a negative relationship between the frequency and duration of pulses in LDC bats which is absent in horseshoe bats (this study, Table 3; [78]). The use of Doppler-shift compensation means that horseshoe bats separate pulse and echo in frequency allowing them to use long duration pulses which enable the detection of insect wingbeats [78].

The phenotype of organisms is the net result of all of the opposing and complimentary evolutionary forces acting upon their lineage over the evolutionary history of that lineage. When an innovative trait arises that confers a substantial advantage in terms of survival and reproduction, the evolution of adaptive complexes between such an innovation with other traits may constrain the evolution of divergence to evolutionary trajectories that do not disrupt such adaptive complexes. HDC echolocation may represent such an innovation.