Phenotypic traits evolution and morphological traits associated with echolocation calls in cryptic horseshoe bats (Rhinolophidae)

Bats provide an excellent casestudy for studying evolution due to their remarkable flight and echolocation capabilities. In this study, we sought to understand the phenotypic evolution of key traits in Rhinolophidae (horseshoe bats) using phylogenetic comparative methods. We aimed to test the phylogenetic signals of traits, and evaluated the best‐fit evolutionary models given the data for each trait considering different traits may evolve under different models (i.e., Brownian Motion [BM], Ornstein‐Uhlenbeck [OU], and Early Burst [EB]) and reconstruct ancestral character states. We examined how phenotypic characters are associated with echolocation calls and minimum detectable prey size. We measured 34 traits of 10 Asian rhinolophids species (187 individuals). We found that the majority of traits showed a high phylogenetic signal based on Blomberg′s K and Pagel′s λ, but each trait may evolve under different evolutionary models. Sella traits were shown to evolve under stabilizing selection based on OU models, indicating sella traits have the tendency to move forward along the branches toward some medial value in equilibrium. Our findings highlight the importance of sella characters in association with echolocation call emissions in Rhinolophidae, as calls are important for spatial cognition and also influence dietary preferences. Minimum detectable prey size in Rhinolophidae was associated with call frequency, bandwidth, call duration, wingspan, and wing surface area. Ultimately, understanding trait evolution requires sensitivity due to the differential selective pressures which may apply to different characteristics.


Introduction
The adaptive radiation of flying mammals is characterized by the success of flying and echolocation, which enabled them to colonize a variety of habitats and by navigating in the dark access a wide variety of niches occupied by birds during the day. This evolutionary success is also amplified by morphological innovations such as wing shape and echolocation emission structures (Jones & Teeling, 2006;Arbour et al., 2019;Amador et al., 2020). Understanding the mechanisms that underlie phenotypic diversity is one major question in evolutionary biology because phenotypic traits of extant species are shaped by their evolutionary history (Losos & Glor, 2003;Münkemüller et al., 2012;Garamszegi, 2014). In this respect, numerous facial traits have been observed in bat species, yet these have been relatively little studied. Echolocation in bats is a complicated task because the call is not only used to locate targets but also to analyze the features of targets, and distinguish prey from nonprey (Neuweiler, 1990). Calls may be optimized to reflect habitat clutter, in addition to temperature and humidity. Horseshoe bats (Rhinolophidae) emit calls with constant frequency signals (CF), and call duration is relatively long (>30 ms) compared to other bat families with CF calls (such as Hipposideridae), and are able to compensate for the effect of Doppler-shift (Schnitzler & Denzinger, 2011). Over the course of evolution, the echolocation behavior and hearing systems of bats have adapted to foraging systems, what they eat and how they acquire food (Koselj et al., 2011;Lazure & Fenton, 2011;Jacobs et al., 2014;Cvikel et al., 2015;Jones et al., 2016;Arbour et al., 2019;Brokaw & Smotherman, 2020;López-Cuamatzi et al., 2020;Giacomini et al., 2021).
Variations in signal emissions often require modification and specialization of structures that produce them (Jones, 1997). It has been proposed that the structure of the noseleaf in nasal-emitters bats could partially explain variations in the echolocation calls (Zhuang & Müller, 2006;Vanderelst et al., 2011). However, while previous studies have focused on the associations between bat skull size and morphology related to echolocation within an evolutionary context (Bogdanowicz et al., 1999;Arbour et al., 2019;Giacomini et al., 2021), little attention to other morphological traits. Body size has been proposed to have constrained echolocation calls in bats, as aerial-feeding bats must be lighter to fly while echolocating, which is energy costly (López-Cuamatzi et al., 2020). Higher frequency species generally have a smaller body size (allometric relationship) to facilitate feeding, locomotion, and reproduction success strategies (Dumont, 2007;Roff & Fairbairn, 2007;Jacobs & Bastian, 2018;Moyers-Arévalo et al., 2020). Forearm and body mass have been used as proximate measures of body size in bats, but some species deviate from the expected allometric relationship (Jacobs et al., 2007;Wu et al., 2015). In these species, the forearm may not be the best predictor of echolocation call frequency (Chornelia et al., 2022). Previous studies show that other features such as ear pinnae are also involved in scanning spatial information using echolocation calls (Mogdans et al., 1988;Kuc, 2010;Yin & Müller, 2019;Leiser-Miller & Santana, 2020), and noseleaf structures are actively involved in shaping emission beams (Arita, 1990;Zhuang & Müller, 2007;Kugler & Wiegrebe, 2017;Zhang et al., 2019). However, most of these studies have focused on Phyllostomidae, and Old-World leaf-nosed and horseshoe bats have been largely neglected, despite the convergence of various characteristics between the two groups. The correlation between ears pinnae and echolocation has been investigated for some Rhinolophidae species in previous studies (Zhao et al., 2003;Wu et al., 2015;Yin & Müller, 2019), and experimental studies highlight the important acoustic functions of various parts of the noseleaf in Rhinolophidae (Zhuang & Müller, 2006, 2007Vanderelst et al., 2010;Feng et al., 2012;Vanderelst et al., 2013), however, these studies have not yet accounted for evolutionary history.
The noseleaf acts like a baffle to redirect and focus the echolocation calls emissions from nostrils and may vibrate to emit acoustic energy more efficiently (Kuc, 2010;Vanderelst et al., 2010). The noseleaf is a spear-with leaf-like structures surrounding nostrils, and may vary in shape and structure, from simple shapes (i.e., Phyllostomidae, Megadermatidae) to more complex structures (i.e., Rhinolophidae) (Vanderelst et al., 2013;Chornelia et al., 2022). The spear-like in posterior components of the Rhinolophidae noseleaf is called the lancet and include multiple furrows (with some species has basal lappets, i.e., additional features in the extension base of sella), and are expected to be adapted to the frequent modulated (FM) component of Rhinolophidae calls (Zhuang & Müller, 2006;Vanderelst et al., 2013). Noseleaf characters have been shown to be effectively used as traits for the identification of potential cryptic species of Rhinolophidae (Chornelia et al., 2022). Even though multiple experimental studies have shown the function of noseleaf in call emission, these characters are still underappreciated for quantitative analysis.
In addition, many species of Rhinolophidae are cryptic, which denotes the species that are morphologically similar and thereby challenging to differentiate based on traditional taxonomic approaches (Mayer & Von Helversen, 2001;Bickford et al., 2006;Stuart et al., 2006;Clare, 2011;Roux et al., 2016;Struck et al., 2018). Cryptic species may be particularly common for taxa that rely on nonvisual cues for mate identification, such as species that rely on olfactory or auditory cues and therefore nocturnal species (especially those which like the Rhinolophidae have a highly conserved call structure). The value of calls as indicators of undescribed species, where echolocation calls may play important role in secondary sexual signals is known from various examples such as Pipistrelles in Europe (Pipistrellus pipistrellus and Pipistrellus pygmaeus) (Jones, 1997), and phonic types in Hipposideros bicolor complex (Kingston et al., 2001). Conversely, some species show regional dialects, such as Hipposideros armiger (Sun et al., 2020) highlighting the need for integrative approaches to use calls as indicators, but not as absolute characters to be used in isolation while verifying species identity. Therefore, the understanding of interrelationships between phylogenetics and bioacoustics in bats is complex, but crucial to study, as bat echolocation is not only constrained by physical features but also reflects the environment and community structure. To the best of our knowledge, apart from single species studies, no study has been conducted more broadly on the evolutionary traits of Rhinolophidae noseleaf characters, although such similar studies have been done largely in the Phyllostomidae family (Leiser-Miller & Santana, 2020;Hall et al., 2021).
Species diversity and phenotypic diversity observed in nature are one of the most challenging tasks to understand from the evolutionary perspective (Fenton, 2010;Garamszegi, 2014). Major evolutionary changes have occurred over extended timescales which limit the sources to estimate the ancestral forms of organisms, therefore, using phylogenetic comparative approaches will enable to draw evolutionary inferences of phenotypic characters in the past (Pagel, 1992;Losos & Glor, 2003;Freckleton, 2009;Münkemüller et al., 2012;Revell, 2012;Garamszegi, 2014;Soul & Wright, 2021). In this study, we sought to understand the phenotypic evolution in Rhinolophidae using phylogenetic comparative methods. Noseleaf characters are known to correlate with acoustic traits, however, how these traits evolve and the relationship between these phenotypic traits and acoustic traits within an evolutionary context remains unresolved. Therefore, our study aims to test the phylogenetic signal, evaluate evolutionary models of each trait given the data, and use the phylogenetic relationship to estimate the expected covariance of one or several predictor variables on a single response variable (acoustic and minimum prey size detectable [MPSD]) while controlling for potential phylogenetic signals in response. Therefore, in this study, we aim to (1) quantify the phylogenetic signals in noseleaf traits, and to analyze the correlation of noseleaf traits, echolocation calls parameters, and external morphological traits including wing. We use two common indices to measure the strength of the phylogenetic signals, Bloomberg′s K (Blomberg et al., 2003) and Pagel′s λ (Pagel, 1992) (as reviewed in Münkemüller et al., [2012]). (2) Considering the traits may not evolve under similar models, therefore, we test the best-fit evolutionary models given to the data for each trait including BM, OU, and Early Burst (EB). (3) We investigated how noseleaf, wing, and external morphological measurements are associated with evolutionary characteristics of echolocation calls. (4) Finally, we explore the relationship between echolocation calls with MPSD within Rhinolophidae. Bat species using shorter wavelengths generally feed on smaller insects (Houston et al., 1999) predicted using Rayleigh scattering calculations to estimate the sizes of insects detectable by different echolocation call frequencies (Jones & Holderied, 2007). Specialization in sensory performance in bats is related to access to food resources such as prey size and to facilitate species coexistence (Safi & Siemers, 2010). Therefore, bats with high frequencies (short wavelengths) are predicted to be able to detect smaller insects, for instance in Vespertillionidae, which supports the hypothesis that call frequency and prey size can be functionally linked (Thomas et al., 2004). We expect call parameters may be associated with MPSD in Rhinolophidae. We also tested other morphological traits parameters that may be linked to minimum detectable prey size. To the best of our knowledge, this study is the first to assess the evolution of noseleaf characters in Rhinolophidae related to acoustic traits and can be used as a basis to explore the phenotypic evolution of complex noseleaf structures for future studies.

Echolocation calls
Echolocation calls were collected from captured individuals in the field and were recorded from hand-released. Calls were recorded using a Pettersson M500-384 (Pettersson Electronic AB; www.batsound.com). A total of 12 210 calls sequences were analyzed using BatSound ver4 (Pettersson Electronic AB) at a sampling rate of 44.1 kHz and the spectrogram was set at 1024 sampling sites in Fast Fourier Transform (FFT) in conjunction with the Hanning window. We measured 5-10 search pulses from each high-quality call sequence, and then averages were taken following methods outlined in Vaughan et al. (1997). Parameters measured included maximum frequency (Fmax), minimum frequency (Fmin) measured from the spectrogram, frequency at maximum energy (FmaxE), or the most intense frequency measured using a power spectrum created using an interpolated (95%) 1024point FFT in conjunction with Hamming windows, bandwidth (BW) (Fmax-Fmin), duration (Dur) was measured from the oscillogram, and interpulse interval (IPI) was measured from the beginning of the call being measured to the beginning of the next call.

Statistical analysis
In this analysis, we focused on 10 species sensu lato, with a total of 187 individuals with potential cryptic species: Rhinolophus rex ( In this study, R. macrotis were treated as one sensu lato species considering that taxonomic status of R. macrotis is complicated and not fully resolved. Previous studies on the R. macrotis complex have split it into three species (Rhinolophus episcopus, Rhinolophus siamensis, Rhinolophus osgoodi) , however, their approach to delineate species is mainly based on genetic data and associated morphological and acoustic variation has not been explored. The two phonic types of R. macrotis in this study are closely related with low mtDNA divergence. R. macrotis 58 kHz size is smaller than R. episcopus (FA = 45.80-48.98 mm/48.8 kHz) and R. macrotis 77 kHz has a similar body size to R. siamensis (FA/calls = 38-42 mm/68-74 kHz) but higher call frequencies and to the best of our knowledge, 74-77 kHz is the highest echolocation calls recorded within R. macrotis complex and has not previously been reported. Clear species boundaries remain challenging particularly when limited genetic data are available, and sampling showed limited geographical coverage. Therefore, in this study, we treated R. macrotis as one species sensu-lato (though we tested the impact of a split on the modified tree topology). Considering all species within the R. macrotis complex diverged around 1-2 Ma, with potential incomplete lineage sorting, mtDNA introgression and hybridization may have occurred, thus can be assumed these species were from a direct common ancestor. The potential cryptic species were defined based on analysis incorporating phylogeny, species delimitation based on multispecies coalescent models, and integrative taxonomy using a combination of acoustic and morphological data based on Chornelia et al. (2022).
First, we tested the selected characters for phylogenetic signal using phylosig function in the R package "phytools" v1.0-1. The time-calibrated tree was generated in BEAST v2.6.3 (Bouckaert et al., 2019), the Maximum Clade Credibility tree was obtained from 45 000 trees sampled in the posterior distribution from BEAST, and summarized in TreeAnnotator v2.6.3, and visualized in FigTree v1.4.4 and was used for further analysis (see details in Appendix S3, tree file is supplemented in Supplementary Material). The tree topology, in general, is consistent with recent studies using large-scale supermatrix genomic tree (i.e., Álvarez-Carretero et al., 2021), however, further increases in the number of species and genomes may improve the systematic stability of the group. The tree branches were pruned prior to the analysis to represent 13 potential cryptic species. A total of 34 characters were included in the analysis, comprising seven external characters, 14 noseleaf characters (lancet, sella, anterior and posterior noseleaf), 11 wings parameters, and two echolocation call parameters. Prior to comparative analyses, the data averages were taken from multiple individuals per species (Appendix S2). We calculated both Blomberg′s K (Blomberg et al., 2003) and Pagel′s λ (Pagel, 1992) indices under the assumption of Brownian Motion (BM) evolutionary models. We hypothesized the similarity of traits in phylogenetically related species is due to its inheritance from a common ancestor, and therefore, contains strong phylogenetic signals. The traits which show strong phylogenetic signals should indicate higher species resemblance to each other compared to other species drawn at random from the tree. Furthermore, we fitted alternative univariate models of trait evolution including the Early-Burst (EB) model and the Ornstein-Uhlenbeck (OU) model using the R package "geiger". Brownian Motion, EB, and OU models were compared using Akaike Information Criterion (AIC), small-sample corrected AIC (AICc), and AIC weighted (AICw). We transformed the tree branches for the traits that showed the non-BM model (i.e., OU) using rescale function in the R package "geiger", before ancestral state reconstruction. The ancestral state characters are then mapped into the pruned tree using fastAnc function in the R package "phytools" 1.0-1 (Revell, 2012).
Second, to test whether echolocation calls (FmaxE/call frequency at maximum energy) are linked with external body characters and noseleaf, we ran a phylogenetic generalized least squares (PGLS) regression using the function pgls from R package "caper". We use echolocation call parameters as a response variable, body size + noseleaf parameters as predictor variables, and pruned time tree generated in BEAST. The variables were log-transformed before the analyses. Finally, we identify the factors associated with minimum prey size with other echolocation calls parameters and wing shape. To estimate the MPSD based on calls frequency, we used the following equation adapted from (Yang, 2010): where λ is the wavelength, v is the speed of sounds in the air (v = 340 m/s), and f is the call frequency. We hypothesized that minimum prey detectable size is correlated with call durations and BW, and also associated with the wing shape of the species.

Phylogenetic signal of selected characters in cryptic Rhinolophidae
We measured the resemblance tendency of closely related species, including the potential cryptic species. In total 20 characters showed relatively strong phylogenetic signals with significantly lower P values (P ≤ 0.05). In external characters, we found that forearm length, body mass, tibia length, headbody length, and tail length contain substantially high phylogenetic signals. In noseleaf characters, we found noseleaf total area, noseleaf length, noseleaf width, lancet height, lancet-tip height, and lancet tip to connecting process height contain relatively high phylogenetic signals. In wing characters, we found wingspan, wing tip index ratio, wing loading, armwing length, arm-wing-area, hand-wing area, and total wing surface area contain higher phylogenetic signals. Additionally, call frequency at maximum energy and BW exhibit relatively higher phylogenetic signals under BM models using Blomberg′ s K and Pagel′s λ indices (Table 1; Fig. 1). The BM′s continuous trait evolution models indicate the changes in the traits are gradual over time, thus, the amount of change of the traits along the branch length of the phylogenetic tree is proportional to the branch length. Bloomberg′s K compares the variance of Phylogenetic Independence Contrast (PIC) under the BM model, K = 1 means that species resemble one another as much as we expect under BM, K < 1 means there is less phylogenetic signal expected under BM and K > 1 means that there is more phylogenetic signal (Blomberg et al., 2003). Moreover, when λ is close to zero meant there is no phylogenetic signal expected in the data under the BM model, λ = 1 correspondent BM and the traits contain phylogenetic signal, and 0 < λ < 1 indicated weaker phylogenetic signal (Pagel, 1992).
However, some characters also show higher values for Blomberg′s K and Pagel′s λ (K > 1; λ > 1; P ≤ 0.05) which indicates stronger trait similarity between related species than expected under the BM model of continuous trait evolution. The characters are body mass, tibia length, tail length, noseleaf length, lancet height, lancet tip height, lancet tip-connecting process height, wing tip index ratio, wing loading, arm-wing length, and arm wing area. The higher indices values (>1.5) indicate that traits did not evolve under BM trait evolution models, and BM is not the best fit model given the data, instead, EB or Accelerated (AC) trait evolution models are the best fit model chosen based on AIC (AICc and AICw). Four characters were shown to evolve under the EB trait evolution model, for instances, body mass (K = 1.886, λ = 1.018, AICw = 0.796), lancet height (K = 2.151, λ = 1.018, AICw = 0.921), lancet tip-connecting process length (K = 2.052, λ = 1.0181, AICw = 0.881), and wing tip ratio I (K = 1.753, λ = 1.076, AICw = 0.644) (Table 2; Fig. 2).
The results suggest different morphological traits evolve under different trait evolution models. The majority of traits evolved under BM models (17 characters), which indicates traits diversification or evolutionary changes of phenotypic along the branches are random, in both direction and distance over any time interval, and OU models (13 characters), suggesting that traits evolve toward adaptive peaks where the intermediate forms which lack adaptive fitness were removed by natural selection. Additionally, only four traits evolved under EB (accelerated) evolutionary models, which indicates the greatest level of early phenotypic divergences which correspond to adaptive radiation. To assess the robustness of the results when we split Rhinolophus macrotis into two tips, we ran an additional test with a modified tree topology (see Section 2). The result of the simulation test demonstrated the consistency of results presented here, including the phylogenetic signals, best-fit trait evolutionary models given the data, as well as ancestral trait reconstructions and correlation between acoustic and phenotypic traits despite the slight modification in the tree topology (results are not presented here). Therefore, in this case, a slight modification in the tree tip did not affect the overall results, however, it is important to determine the uncertainty associated with the tree topology by repeating the analysis using the alternative topologies.
We found relatively weaker relationships between echolocation calls and wing parameters, with R 2 = 0.5 in wingspan and wing surface area (PGLS; P = 0.002284, R 2 = 0.5861; P = 0.003, R 2 = 0.5506, respectively) ( Table 3). Wing tip index, wing tip ratio, and wing loading were not significant predictors for call frequency, including other morphological traits such as tail length, lancet height, and lancet tip-connecting process length.

Phylogenetic signal and trait evolution of selected characters in cryptic Rhinolophidae
Our results demonstrate different patterns of trait evolutionary models in Rhinolophidae. Among 34 traits measured and included in this study, 20 characters were shown to contain strong phylogenetic signals based on Bloomberg′s K and Pagel′s λ indices, which is expected under BM models (0 > K > 1; 0>λ > 1; p ≤ 0.05). However, some traits are a better fit with other evolutionary trait models, such as OU and EB, indicating the traits may evolve in different ways, which may be expected due to the various selective pressures that each trait may be subject to. OU is a model of adaptive evolution that includes the maintenance of traits by stabilizing selection as the dominant evolutionary force with a constant pull toward an optimum value (Hansen, 1997;Beaulieu et al., 2012). EB is a model with adaptive radiation scenarios, which initially show rapid morphological evolution in the early evolutionary stage followed by a slow evolutionary rate and relatively stasis/stable adaptive peak (Harmon et al., 2010;Ingram et al., 2012;Puttick, 2017). The two indices are assumed to show continuous trait evolution, BM process (random evolutionary walk), however, for some cases, BM may not be the best model when there is consistent selection toward a single optimum trait value, that is, the OU model (O′Meara et al., 2006;Beaulieu et al., 2012). Alternative models can incorporate a range of processes other than simple random walks BM, such as punctuational and speciational evolution, accelerating and decelerating rates of evolution, and adaptive radiation (Pagel, 1992;Hansen, 1997;Freckleton, 2009).
Our results suggest that body mass, lancet height, lancet tip-connecting process, and wing tip index evolved under EB evolutionary models, suggesting these traits diversify quickly as the species enter new niches and the rate of evolution slows through time (Arbour et al., 2019). EB or Accelerated (EB/AC) model evolution is rarely observed in evolutionary processes, in contrast, BM includes constant and random-walk evolutionary processes and is the most common evolutionary model of morphological traits. Although, previous studies found that maximum morphological disparities tend to establish early in evolutionary history (Hughes et al., 2013), for instance, skull phenotypic diversification in bats (Giacomini et al., 2021).
Additionally, other than BM models, most of the characters show some degree of fit with OU trait evolutionary models, particularly for traits that contained lower phylogenetic signals measured from Blomberg′s K and Pagel′s λ. Interestingly, the majority of characters that better fit with the OU trait evolutionary model are sella characters including internarial cup, sella area, sella height, lancet base width, lancet angle, and sella base width among the species we included in the study. Meanwhile, other noseleaf characters such as noseleaf area, noseleaf length and width, and noseleaf tip height were better fit with BM models. Our results suggested that sella traits have a tendency to move forward along the branches toward some medial value in an equilibrium (Hansen, 1997), in other words, other intermediate forms and sella trait diversification may be eliminated by natural selection, as the shape of the sella is likely to relate to its function. Previous studies showed evidence of natural selection in the fiber-type composition of the locomotor muscle of lizards based on the OU model (Scales et al., 2009), and indicated a common selection pressure across the clades. This emphasizes the importance of understanding trait evolution in bats and optimizing the selection of traits that have often been neglected despite their importance. The majority of the previous evolutionary trait studies on bats focus on skulls and dentition (Jacobs et al., 2014;Clavel & Morlon, 2020;Esquivel et al., 2021;Giacomini et al., 2021), but few studies have been focussed on measuring the noseleaf traits, particularly in Old-World bats.

The echolocation calls evolutionary associated with phenotypic traits
By incorporating phylogenetics into the analysis, we consistently showed most noseleaf and sella characters were significantly associated with the evolution of call fmaxe. This result is supported by previous findings in (Chornelia et al., 2022). The most common parameters used to scale the relationship between body size and echolocation calls among bats are forearm and body mass (Bogdanowicz et al., 1999;Jacobs & Bastian, 2018), however, many species groups including Rhinolophidae showed deviation from the allometric rules between echolocation calls and forearm size (Jacobs et al., 2007;Wu et al., 2015;Chornelia et al., 2022). Furthermore, lower coefficients in previous studies suggest that factors other than body mass and body size may influence the call frequency of maximum energy of bats (López-Cuamatzi et al., 2020), and we found noseleaf-metric and sella characters are better predictors of the evolution of echolocation calls. In this study, we found forearm and body mass (PGLS; P << 0.001, R 2 = 0.6639; P = 0.01, R 2 = 0.4649) showed relatively lower associations than other noseleaf and sella traits (PGLS, R 2 ≥ 0.72). The combination of all traits in sella and noseleaf also showed a relatively higher association with the evolution of call frequency (PGLS; P = 0.001, R 2 = 0.9965), and combined with ears size and forearm almost give a similar result (PGLS; P = 0.03, R 2 = 0.9998). Our results also suggest that wing shape such as wingspan and total area of wings are related to call frequency (PGLS; P = 0.002, R 2 = 0.5861; P = 0.003, R 2 = 0.5506), however, the values are relatively low compared to noseleaf and sella traits. However, it is also worth noting that lower R 2 does not necessarily mean the models are not the best, but there are several possible reasons why R 2 was low, for instance, random noise in the data due to measurement errors (Freckleton, 2009), thus, here we only use R 2 as a diagnostic to assess lack of model fit.
In bats, a recent study proposed "nested-constraint hypothesis" which suggested echolocation and body mass are strongly linked, and flight in conjunction with echolocation constrains the body mass in echolocating bats (Moyers Arévalo et al., 2020). This hypothesis suggests that negative relationships between frequency and body mass may be associated with efficient echolocation for small bats since the early evolution of the clade, in contrast with some larger bats that lost their echolocation ability, such as Pteropodidae. Echolocation might have subsequently been lost in Pteropodidae (with tongue-clicking or wing-clapping evolving secondarily), which share common ancestry with Rhinolophidae, though Rhinolophids have sophisticated echolocation and acoustic features (i.e., Doppler-shift effect compensation, High-Duty Cycle calls). However, whether echolocation evolved once or twice, and if echolocation arose independently in rhinolophoids remains an open question (Jones & Teeling, 2006).
While this hypothesis may be applied in the Chiroptera as a whole, there is a deviation from the rule in family groups such as Rhinolophidae, where more complex mechanisms enable greater modification of call structure. Deviations from the expected allometric relationship have been recorded in a number of species, such as H. larvatus (body mass~19 g, calls frequency~86.5 kHz) (Furey et al., 2009) and R. rex (body mass~10 g, calls frequency~25 kHz) (Chornelia et al., 2022), including other rhinolophids such as R. clivosus, R. fumigatus, and R. simulator (Jacobs et al., 2007). Deviation from the allometric hypothesis has also been found in Phyllostomidae. Previous studies also suggest that other morphological features, such as vocal tract geometry (Neuweiler, 1990) or noseleaf morphology (Reijniers et al., 2010;Vanderelst et al., 2010Vanderelst et al., , 2013 may be better predictors of emitted frequency than body size. For example, nasal capsules may also act as better predictors for calls than body size for nasalemitting species as it negatively correlated with calls (Jacobs et al., 2014) even though nasal capsule volume alone cannot explain the anomaly of high echolocation than expected in R. clivosus. Here, we propose the sella area is strongly associated with echolocation calls, suggesting that the sella size progressively wider and taller as call frequency decreases.

Predicting the MPSD associated with echolocation calls
Traits associated with feeding, locomotion, and reproduction in bats have been shown to be explicitly linked with call structure (Hall et al., 2021). A metaheuristic algorithm inspired by bats' acoustic behavior, namely Bat Algorithm (BA) generated a simple mathematical formula to estimate bat prey sizes, which is equal to the speed of sound in air (v = 340 m/s) divided by maximum frequency (FmaxE). The algorithm was generated by considering that each ultrasonic emission may last from 5 to 20 ms, therefore, around 10-20 pulses are emitted every second. When hunting, calls may be sped up to 200 pulses per second, which Note: Lambda measured the covariances between species, λ = 0 correspondence to nonphylogenetic regression, λ closer to 1 indicates the evolution of residual error is Brownian; PGLS, phylogenetic generalized least squared. Note: Lambda measured the covariances between species, λ = 0 correspondence to nonphylogenetic regression, λ closer to 1 indicate the evolution of residual error is Brownian; MPSD, minimum prey size detectable; PGLS, phylogenetic generalized least squared.
demonstrates the incredible ability of signal processing about 300-400 µs in bats' ears (Yang, 2010). Our result suggested the wavelength of Rhinolophidae calls is in the range of 3-14 mm, which is expected in the same order as prey size. Other than the frequency at maximum energy (FmaxE), we determined that prey size was associated with ears size, body mass, frequency maximum, frequency minimum, BW, and call duration. This is in agreement with previous experimental studies in Rhinolophidae (i.e., Vanderelst et al., 2011). Rhinolophidae ears are important in capturing the echoes from the emitted calls by actively moving the ear pinnae during echolocation (Mogdans et al., 1988;Firzlaff & Schuller, 2004;Vanderelst et al., 2011;Yin & Müller, 2019). Moving ear pinnae can sometimes be able to detect glints (frequency and amplitude produce by fluttering insect wings) (Siemers & Ivanova, 2004;Koselj et al., 2011;Lazure & Fenton, 2011;Schnitzler & Denzinger, 2011), because narrow bandwidth of CF call frequency prevents Rhinolophidae from using spectral cues to localize the reflectors (Firzlaff & Schuller, 2004). Therefore, moving ears pinnae received the incoming echoes from insect glints as amplitude cues which provide stable localization information and estimate the origin of the echoes (Vanderelst et al., 2011). This also indicates that Rhinolophidae only hunts insects that introduce glints in their echoes, that is, fluttering insects such as moths (Schnitzler et al., 2003;Lazure & Fenton, 2011;Schnitzler & Denzinger, 2011). The long-duration calls in Rhinolophidae can distinguish smaller insects from larger ones because smaller insects produce more acoustic glints (small modulations in amplitude and frequency) (Jones & Teeling, 2006). Longer duration of calls may be able to detect larger prey. The association of MPSD with wing was only shown in wingspan and wing area, but no significant relationships were observed with wing-tip index, aspect ratio, and wing loading. This result suggests smaller species with lower call frequencies such as R. rex may be able to detect larger prey compared to medium-sized rhinolophids with high echolocation calls, however, further experimental studies are needed to confirm this pattern. In addition, highfrequency sound attenuates more rapidly, thus, highfrequency calls attenuate too fast to be useful over long distances (Lawrence & Simmons, 1982). This trade-off means that species with high calls have to detect prey and their reaction time would be shortened, increasing energetic cost, and meaning that such species may have more specific prey requirements as their range of detectability will be smaller.

Future directions
Research into the evolution of phenotypic traits in bats is clearly in the preliminary stages. Nevertheless, previous studies have shown interesting findings on the evolutionary relationship among morphological traits, echolocation, ecological demands, functional traits, and bat behavior (i.e., Jacobs et al., 2014;Brokaw & Smotherman, 2020;López-Cuamatzi et al., 2020;Esquivel et al., 2021). Yet, gathering complex data for greater systematic sampling to achieve robust results on these inter-relationships remains a challenge. Some caveats should be acknowledged, first, we use a smaller subset of taxon sampling within Rhinolophidae due to limited noseleaf data available from previous studies. Therefore, we recommend increasing systematic sampling including a breadth of traits (phenotypic, acoustics, and genetics) in future studies. We did not attempt to fit more complex evolutionary models on the continuous traits, as fitting complex evolutionary models require bigger character dataset within larger phylogenetic trees and is usually computationally intensive. However, the results presented here provide a baseline for future research. In addition, unstable systematics and tree topologies potentially influence the phylogenetics comparative analysis in general especially if it involves species that are moved across the root of the tree (Blomberg et al., 2012). However, changes that happen closer to the tip, as in the case of Rhinolophus macrotis here, cause less drastic effects, nevertheless, controlling the uncertainties in phylogenetics is needed to assess the robustness of the results (Martins & Housworth, 2002;Symonds, 2002;Garamszegi & Gonzalez-Voyer, 2014). Thus, future work may incorporate larger datasets with a complete time-calibrated phylogeny. Finally, fitting evolutionary models to data and model selection based on AIC does not necessarily mean that it is a perfect model for the data, but it does highlight the relationship between these traits. Further, we should note that the possible errors in tree topology are likely to bias the interpretation of phylogenetic comparison (Pagel, 1992;Blomberg et al., 2003;Freckleton, 2009). Phylogenetic reconstructions for bat clades are still generally derived from a limited number of genes, and therefore, have a degree of uncertainty. However, generating super-matrix time-trees based on genomic analysis is becoming increasingly possible to overcome these constraints (Shi & Rabosky, 2015;Amador et al., 2018;Álvarez-Carretero et al., 2021). Furthermore, such data will make it possible to test larger morphological and acoustic datasets, and in improve the accuracy of time-tree estimation.

Conclusion
Our study highlights the importance of understanding the phenotypic evolutionary processes in bats. The results suggest that different morphological traits may have evolved under different evolutionary models, likely due to different selection pressures on different characteristics. Interestingly, most of the sella characters measured from representative species in this study evolved under OU models, which suggests that sella characters evolved toward equilibrium or stabilizing selection, and any transition or intermediate forms of sella were eliminated by natural selection or other pressures. Noseleaf metrics are better predictors of echolocation calls than body size or forearm length in Rhinolophidae when phylogenetics are incorporated into the analysis. As echolocation-related traits are associated with the foraging strategy, we found that MPSD is associated with ear size and other echolocation parameters such Fmax, Fmin, BW, and duration. Ultimately call structure relates to its multifunctional use, as well as being impacted by the physical structures responsible for its production. Due to high levels of conservatism in some call elements, calls can provide an initial clue to the possible existence of cryptic species, which may be differentiable on the basis of various noseleaf traits.