Environmental correlates of geographic divergence in a phenotypic trait: A case study using bat echolocation

Abstract Divergence in phenotypic traits may arise from the interaction of different evolutionary forces, including different kinds of selection (e.g., ecological), genetic drift, and phenotypic plasticity. Sensory systems play an important role in survival and reproduction, and divergent selection on such systems may result in lineage diversification. Such diversification could be largely influenced by selection in different environments as a result of isolation by environment (IbE). We investigated this process using geographic variation in the resting echolocation frequency of the horseshoe bat species, Rhinolophus damarensis, as a test case. Bats were sampled along a latitudinal gradient ranging from 16°S to 32°S in the arid western half of southern Africa. We measured body size and peak resting frequencies (RF) from handheld individual bats. Three hypotheses for the divergence in RF were tested: (1) James’ Rule, (2) IbE, and (3) genetic drift through isolation by distance (IbD) to isolate the effects of body size, local climatic conditions, and geographic distance, respectively, on the resting frequency of R. damarensis. Our results did not support genetic drift because there was no correlation between RF variation and geographic distance. Our results also did not support James' Rule because there was no significant relationship between (1) geographic distances and RF, (2) body size and RF, or (3) body size and climatic variables. Instead, we found support for IbE in the form of a correlation between RF and both region and annual mean temperature, suggesting that RF variation may be the result of environmental discontinuities. The environmental discontinuities coincided with previously reported genetic divergence. Climatic gradients in conjunction with environmental discontinuities could lead to local adaptation in sensory signals and directed dispersal such that gene flow is restricted, allowing lineages to diverge. However, our study cannot exclude the role of processes like phenotypic plasticity in phenotypic variation.

and drift (isolation by distance, IbD; Wright, 1943) in combination with selection for environment-specific traits (isolation by environment, IbE a subset of isolation by adaptation; Safran et al., 2016). Both isolation and selection are likely to vary spatially along a continuum of environmental gradients (e.g., distance and climate). However, traits that are heritable and have a profound impact on fitness (e.g., sensory traits) are more likely to be impacted by natural selection than drift, resulting in adaptation to novel environments (Boughman, 2002;Campbell et al., 2010;Endler, 1992). Nevertheless, the role of phenotypic plasticity in phenotypic variation cannot be excluded (Ghalambor, McKay, Carroll, & Reznick, 2007;Thibert-Plante & Hendry, 2011;Via et al., 1995). Phenotypic plasticity can influence variation if it promotes successful dispersal and survival in different environments (Pfennig et al., 2010). Such dispersal and resultant gene flow may favor the evolution of increased phenotypic plasticity, over adaptive genetic divergence, because it would promote adaptation to new environments within a few generations (Pfennig et al., 2010), thereby promoting phenotypic variation.
Sensory signals (e.g., acoustic signals) are important in the context of lineage diversification (Mutumi, Jacobs, & Winker, 2016;Odendaal, Jacobs, & Bishop, 2014;Slabberkoorn & Smith, 2002) because they are essential to the survival and reproduction of animals that rely on such signals for orientation and foraging as well as in mate choice and assortative mating (Bolnick & Kirkpatrick, 2012;Coyne & Orr, 2004).
These signals are propagated through the environment, and specific environmental conditions can influence the evolutionary trajectory of the signaling system. Acoustic signals, in particular, are part of a sensory system that relies on audition and the propagation of sound through the atmosphere.
Atmospheric conditions (i.e., climate) can therefore exert a strong influence on geographic variation of complex signals, such as bird song (Lengagne & Slater, 2002), frog mating calls (Lingnau & Bastos, 2007), and the echolocation calls of bats (Lawrence & Simmons, 1982;Luo, Kosel, Zsebok, Siemers, & Goerlitz, 2014). These acoustic signals may diverge along climatic gradients as a result of variation in atmospheric attenuation of sound. Atmospheric attenuation, caused by the scattering and absorption of sound by the atmosphere, is the result of a complex interaction between the humidity and temperature of the air as well as the frequency of the sound (Hartley, 1989; Lawrence & Simmons, 1982;Luo et al., 2014;Mutumi et al., 2016).
Climate could therefore play a pivotal role in driving the evolution of signaling systems through its effect on atmospheric sound attenuation (Griffin, 1971;Richards & Wiley, 1980;). For example, climatic conditions were found to have influenced the echolocation call frequency of Hiposideros ruber (Guillén, Juste, & Ibáñez, 2000) and some horseshoe bat species (Mutumi et al., 2016) through atmospheric attenuation of acoustic signals. Furthermore, other studies have suggested that bat echolocation call frequency diverges along climatic gradients (Snell-Rood, 2012).
Such environmentally mediated differences in sensory systems leading to lineage diversification can be facilitated by geographic or environmental isolation of populations (Schluter, 2001). Isolation by geographic barriers and distance can influence trait variation by restricting gene flow between populations. This is usually manifested by an association between genetic or phenotypic differences and geographic distances and referred to as isolation by distance (IbD; Wright, 1943).
Alternatively, when a species is distributed across several biomes, discontinuities in habitat can also result in environmentally mediated ecological isolation of local populations and a restriction of gene flow among them. Gene flow can be restricted in such situations either by selection against dispersers moving between populations in different habitats or by individual preference for remaining in a particular habitat (Nosil, 2004;Nosil, Vines, & Funk, 2005). This is usually manifested by an association between phenotypic differences and environmental differences (IbE, Crispo, Bentzen, Reznick, Kinnison, & Hendry, 2006;Nosil, 2012;Rundle & Nosil, 2005;Shafer & Wolf, 2013).
Bat echolocation has evolved primarily for orientation (Schnitzler et al., 2003;Schuchmann & Siemers, 2010;Simmons & Stein, 1980) and food acquisition (Heller & Von Helversen, 1989). There is ample evidence that the frequency of bat echolocation is adapted to the habitat and foraging mode of bats (Aldridge & Rautenbach, 1987;Jones & Holderied, 2007). There is also some, but not conclusive evidence, that it is secondarily involved in communication (Bastian & Jacobs, 2015;Knörnschild, Jung, Nagy, Metz, & Kalko, 2012;Siemers, Beedholm, Dietz, Dietz, & Ivanova, 2005). Echolocation thus provides an opportunity to investigate the role of evolutionary processes in the geographic variation of a phenotypic trait that is essential to survival and reproduction. However, only a few studies have investigated geographic variation in resting frequency and the environmental factors (humidity and temperature) responsible for it (Luo et al., 2014;Mutumi et al., 2016).
A potentially confounding factor in understanding the effect of climate on geographic variation of resting frequency is the inverse correlation between body size and echolocation frequency in bats (Jones, 1996(Jones, , 1999. Larger bats tend to have lower frequencies, because they have longer and thicker vocal chords as well as larger resonant chambers (Jacobs, Barclay, & Walker, 2007). James' Rule (JR) proposes that animals occurring in hot humid environments generally have smaller body sizes than animals of the same species that occur in cooler, dryer environments, and the largest animals are expected to occur in cool, dry areas (James, 1970). This is because selection at cooler environments would favor larger individuals with lower surface area-tovolume ratio, and large body size enables individuals to conserve heat in cooler environments (Meiri & Dayan, 2003). Any variation in echolocation frequency might therefore simply be the result of the response of body size to different temperatures and humidity.
The southern Africa horseshoe bat, Rhinolophus damarensis Roberts, 1946, has a wide distribution in southern Africa which stretches from the more arid regions in northwestern South Africa and southern Namibia to the more mesic regions of northern Namibia and southern Angola (Jacobs et al., 2013). Genetic analyses based on mitochondrial cytochrome b and the nuclear thyrotropin beta chain precursor indicated that this species consists of two lineages, a northern lineage restricted to the more mesic regions of northern Namibia and a southern lineage restricted to the more arid regions of central and northwestern South Africa, extending into central Namibia. The phylogenetic analyses on cyt b revealed that the two lineages split ~4.8 MYA. There was a 4.1% sequence divergence between the northern and southern lineages with much lower within lineage divergence of 0.7% and 0.5%, respectively (Jacobs et al., 2013). Rhinolophus damarensis therefore provides an ideal opportunity for testing the roles of spatial separation and environment on the diversification of an acoustic signal. If IbE plays a dominant role in echolocation frequency divergence then: (1) call frequency variation should be correlated with climatic variables across populations; (2) call frequency variation should be partitioned in accordance with regional groupings; that is, there should be a strong signal of IBE between the northern and southern lineages; and (3) if call frequency is the result of selection rather than lineage-specific differences (e.g., arising by genetic drift), the correlation between call frequency and climate should also be evident within a particular region and in the same direction as the correlation across all populations. Furthermore, call frequency variation should not be correlated with geographic distance; that is, there should be no signal for IbD. Under James' Rule, we predicted that there should be a correlation between body size and climatic factors and between body size and call frequency. We currently lack data on gene flow and dispersal to adequately test the role of phenotypic plasticity and other evolutionary processes, but we discuss their potential influence in promoting phenotypic variation.

| Study animal
Rhinolophus damarenesis (Jacobs et al., 2013) is a small insectivorous bat with a relatively high echolocation frequency (85.4 ± 1.4 kHz) and a mean forearm length of 49.5 ± 1.7 mm (Jacobs et al., 2013). It has a wide distribution in the western half of southern Africa stretching from western South Africa through Namibia to southwest Angola (Jacobs et al., 2013). This region is characterized by a wide range of arid and mesic climatic conditions. Like other horseshoe bats (Rhinolophidae), R. damarensis uses high duty cycle (HDC) echolocation. Duty cycle is the ratio of call duration to call period and is usually expressed as a percentage. Rhinolophids have several harmonics in their echolocation calls but invariable concentrate call energy in the second harmonic (Fenton, 1994;Hartley & Suthers, 1988; Figure 1). Their calls ( Figure 1) have a prominent constant frequency (CF) component preceded and followed by a frequency modulated (FM) component (Neuweiler, 1984). Horseshoe bats are perch hunters that are able to echolocate from a resting position while scanning the environment for prey (Jacobs et al., 2007). The peak frequency (frequency of maximum amplitude) of the CF component of a call emitted while at rest is called the resting frequency (RF).
This RF is 100-300 Hz lower than the reference frequency, a narrow range of frequencies to which the acoustic fovea (a region of overrepresentation of neurons in the auditory cortex) is tuned (Schuller & Pollak, 1979). Resting frequencies, which are the calls used when foraging from a perch, can therefore be recorded from handheld bats. This eliminates variance in frequency caused by horseshoe bats compensating for Doppler shifts in the returning echo during flight (Schnitzler & Denzinger, 2011).

| Ethics statement
Methods used in this research were carried out in accordance with guidelines of the American Society of Mammalogists (Gannon & Sikes, 2007). We followed the sampling techniques/guidelines outlined by Kunz and Parsons (2009), and these were approved by the Science 1197/2008). Lands accessed to sample bats were publicly or privately owned and sometimes protected. At all times, the necessary permission was obtained. We did not sample protected or endangered species.

| Sampling
Data were collected across the known range of R. damarensis, in southern Africa from seven localities along a latitudinal gradient ranging from 16°S to 32°S (Figure 2). The geographic coordinates (latitude and longitude) of each locality at which we sampled R. damarensis were F I G U R E 1 A typical call of Rhinolophus damarensis. Top = oscillogram; Bottom = spectrogram recorded using a Garmin GPS (model Colorado, Garmin International Inc, Kansas). Bats were captured from their roosting areas such as caves, mines, abandoned buildings, old mine shafts, hollow trees, and culverts under roads using hand nets during the day, or mist nets and harp traps as bats emerged from the roosts at dusk. Mist nets and harp traps were checked regularly throughout the trapping period to ensure that bats were not trapped in these for too long. All captured bats were transferred directly into soft cotton holding bags.

| Morphology
The body mass of each bat was measured using a portable electronic balance (to the nearest 0.01 g), and forearm length (FA) was measured to the nearest 0.01 mm using digital calipers. Juveniles were excluded from this study and were identified by the presence of epiphyseal plates in their finger bones. The plates were detected by transilluminating the bats' wings (Anthony, 1988). Forearm length instead of mass was used as an indicator of body size because mass varies diurnally, seasonally and with life-history stage (Jacobs et al., 2013;Mutumi et al., 2016). The sex of each bat was also recorded.

| Echolocation
Echolocation calls were recorded from handheld bats positioned 10 cm in front of an Avisoft Ultrasound Gate 416 (Avisoft Bioacoustics, Berlin, Germany) microphone connected to a HP Compaq nx7010 notebook computer running Avisoft SasLab Pro software (Avisoft Bioacoustics Schönfließer, Germany) with the sampling rate set at 500 kHz. We measured 10 high-quality calls (high signal-to-noise ratio) that occurred after the first 10 calls in each recording because horseshoe bats tune into their peak resting frequency after a period of silence (Schuller & Suga, 1976;Siemers et al., 2005). Calls were analyzed in BatSound Pro (version 3.20; Pettersson Elektronik AB, Uppsala, Sweden) using a sampling rate of 44,100 Hz (16 bits, mono) and a threshold of 15. We measured the peak resting frequency (kHz) from the fast Fourier transformation (FFT) power spectrum (size set at 1,024 samples) of the dominant 2nd harmonic in each call for 10 high-quality (high signal-to-noise ratio) calls. The average RF over the 10 calls for each bat was used in the analyses.

| Environmental variables
ArcGIS shape files for climate data were downloaded from BIOCLIM (http://www.worldclim.org/bioclim) and OEI (www.en.openei.org) websites. The shape files were analyzed in ArcGIS v.10 to extract data on relative humidity (RH) and annual mean temperature (AMT).
Relative humidity was based on data taken at 10 m above the surface of the earth by NASA Surface meteorology and solar energy (SSE Release 6.0, Data Set, November 2007), a 22-year (1983-2005) monthly and annual average dataset (http://eosweb.larc.nasa.gov/ sse/). Annual mean temperature shape files were based on monthly temperature values that were averaged over the period of 50 years . All environmental data were extracted for a radius of 10 km around each locality at which we sampled R. damarensis. A radius of 10 km was used because the home range of rhinolophids of similar size to R. damarensis (e.g., R. euryale and R. ferremequinum) has been measured at within 10 km (Bontadina, Schofield, & Naef-Daenzer, 2002;Flanders & Jones, 2009;Goiti, Aihartza, Garin, & Zabala, 2003;Jones & Morton, 1992). Currently, there are no data on the home range of R. damarensis.

| Atmospheric attenuation and detection range
We calculated the mean atmospheric attenuation for each locality and then used these to calculate the mean detection ranges of prey (PDR) for each population of R. damarensis using the web calculator developed and described in Stilz and Schnitzler (2012). Prey detection ranges were calculated to determine the impacts of ecological constraints F I G U R E 2 Major biomes of southern Africa (Rutherford, 1997)
The calculations required the following information: (1) climatic conditions (e.g., AMT, RH, and atmospheric pressure) at each sampled site, where atmospheric pressure was kept at that for normal atmospheric conditions, taken as 101.325 pascal; (2) resting frequency (Hz) of each individual bat; (3) the dynamic range, which is the difference between peak intensity (dB SPL) measured at 1 m and the auditory threshold of the bat (assumed to be 0 dB SPL for horseshoe bats (Holderied & von Helversen, 2003;Long & Schnitzler, 1975); (4) reflection loss, C 1 , which accounts for the fraction of the energy reflected, and (5) geometric spreading, C 2 , which quantifies the loss of energy due to spreading multiplied by 2 for both outgoing emitted call and the returning echo.
The values of the latter two factors are dependent on the geometry of the reflected wave. The web calculator generates C 1 and C 2 depending on the target selected. We chose the point function reflector which resembles the insect prey of bats. We used a call intensity of 123.7 dB SPL for (c) above, measured by Fawcett, Jacobs, Surlykke, and Ratcliffe (2015) for a horseshoe bat of similar size (R. capensis). The actual call intensity of R. damarensis is currently unknown. We compared the calculated PDRs for the different localities using the Kruskal-Wallis test followed by multiple comparisons.

| Predictor variables
To investigate the plausibility of the James Rule, IbE, and IbD, we evaluated candidate models with different combinations of environmental, spatial, and biological predictor variables to determine their influence on RF across the distributional range of R. damarensis.
Biological predictor variables included forearm length (FA), which was used as a proxy for body size (James' Rules) and gender which was used to determine whether there was sexual dimorphism in RF within

R. damarensis.
The four alternative environmental variables considered were annual mean temperature (AMT), relative humidity (RH), atmospheric attenuation (AA), and prey detection range (PDR). Temperature and relative humidity have previously been suggested as factors largely responsible for atmospheric attenuation of sound (Hartley, 1989;Lawrence & Simmons, 1982;Luo et al., 2014 Similarly, the predictors region (Reg) and latitude (Lat) were considered as alternative covariates in separate candidate models. The factor Reg comprised two levels representing two genetically differentiated lineages within R. damarensis (Jacobs et al., 2013). One lineage included populations that were located north of latitude −24ºS and was designated the northern lineage. The other lineage included populations further south than latitude −24ºS and was designated the southern lineage ( Figure 2). By contrast, the covariate Lat represented a linear predictor, which would imply a continuous geographic divergence in RF instead of a geographic break. The split between the two R. damarensis lineages was considered to provide insight into potential spatial divergence of RF in R. damarensis.

| Statistical procedure
Geographic variation in RF may be the result of stochastic processes (e.g., genetic drift). One way in which such stochastic processes would be manifested is in a positive relationship between the RF of populations and the geographic distance between them. We therefore also determined whether differences in RF were associated with geographic distances by calculating a dissimilarity matrix of RF (kHz) differences among localities and Euclidean distances among populations from their geographic coordinates (longitude and latitude). A simple pairwise Mantel test was employed to analyze associations between RF differences and geographic distances (Legendre & Fortin, 2010).
We used 10,000 permutations based on Monte Carlo simulation tests.
We applied linear mixed-effects models (LMEs) to analyze the relationship between the response variable (RF) and the environmental (AMT or RH or AA or PDR), spatial (Lat or Reg) and the biological (FA, sex) predictor variables. Linear mixed-effects models incorporate both fixed and random effects (Verbeke & Lesaffre, 1996;Zuur, Ieno, Walker, Saveliev, & Smith, 2009), where the fixed effects quantify the overall effects across the different localities and the random effects quantify variation of the fixed effect parameters across localities (Bolker et al., 2008). The random effect for sampling location was included to account for the nested sampling design as a result of sampling several individuals from a single location. In addition, LMEs can incorporate autocorrelation in the data, which was considered here to account for potential spatial dependencies among sampled localities (Bjørnstad & Falck, 2001) to ensure that IbE is not being driven by spatial autocorrelation (Shafer & Wolf, 2013). It was also for this reason that we investigated IbD (Shafer & Wolf, 2013).
We further explored whether selection or stochastic processes were responsible for RF variation by running separate LMEs on just the southern populations. If RF variation is the result of selection under different climates, the correlation between climatic variables and RF among all populations should also exist and be in the same direction among the southern populations.
We first used a reasonably complex set of variables (AMT, Reg, and Sex with resting frequency as the response variable) to determine the best mixed-effects structure (Zuur et al., 2009). Llinear mixed-effects models were fitted with and without spatial correlation structures and a random effect for sampling locality. The candidate LMEs were fitted with various spatial autocorrelation functions on longitude and latitude coordinates, including corGaus, corExp, corRatio, corSpher, and AR1. However, the most parsimonious error structure, as judged by the Akaike's information criterion (AIC), was a linear mixed model (LME) that only included a random effect for locality. This was also supported by nonparametric spatial spline correlograms based on the model residuals, which showed that spatial autocorrelation, evident for a fixed-effects linear model, could be effectively accounted for by including the random effect for locality (Figure 3). Inspection of residuals from the random-effects LME showed that the residuals closely approximated a normal distribution and provided no evidence for violation of the assumption of homogenous variance (Figure 4). The LME with a random effect for locality was then selected as the most parsimonious error model for analyzing the variation in RF.
To determine the optimal combination of covariates, competing models were ranked based on their AICc, which is the AIC corrected for a small sample size (n) relative to the number of parameters being estimated (K). We ran a total of 40 model permutations to examine which model subsets best explained the variation in our data. Akaike information criterion corrected differences (Δ i ) and Akaike weights (w i ) were used to identify the best candidate models (Amar, Koeslag, Malan, Brown, & Wreford, 2014;Odendaal et al., 2014). Models with the lowest AICc and highest Akaike values were considered the most parsimonious, and the differences in AIC scores (Δ i ) were calculated to determine the likelihood that a given model was the best approximating model relative to other candidate models (Amar et al., 2014;Odendaal et al., 2014;Symonds & Moussalli, 2011). A model with a (Δ i ) value of zero was considered the best approximating model, and those with values of up to two were regarded as good as the best model (Symonds & Moussalli, 2011). Evidence ratios of the best fitting models were calculated relative to the other subsequent candidate models to determine the relative evidence for the best approximating model in relation to the other candidate models (Amar et al., 2014).
The predictions of James' Rule were tested by incorporating forearm in our models because body size can influence resting frequency variation in echolocating bats due to the allometric relationship between body size (e.g., forearm) and resting frequency (Jacobs et al., 2007;Jones, 1996Jones, , 1999. All our analyses were conducted in R version 3.1.2, and the following packages were used: "AICcmodavg" for model selection and multimodel inference to compare and rank multiple competing models and to estimate those that best approximate the true processes underlying the biological phenomenon under study (Symonds & Moussalli, 2011); "lme" for fitting the linear mixed-effects model and to incorporate both fixed and random effects in a linear predictor expression from which the conditional mean of the response was evaluated (Bates, Mächler, Bolker, & Walker, 2015); "effects" for displaying the effect sizes of linear, generalized linear, and other models (Fox & Andersen, 2006); "car" a companion to applied regression for regression diagnostics and other regression-related tasks (Fox, 2002); "ncf" a spatial nonparametric covariance function to make correlograms and to check for spatial autocorrelation (Bjørnstad & Falck, 2001); and "Ade4" for estimating geographic distances between sites (Legendre & Fortin, 2010;Thioulouse, Chessel, Dole′Dec, & Olivier, 1997).

| RESULTS
We recorded and analyzed calls from 106 adult R. damarensis individuals from seven localities (Table 1). This species had an average RF of 85.43 ± 1.3 kHz and average call duration of 31.14 ± 5.9 ms. Mean RF ranged from 84.4 ± 0.7 kHz to 87.6 ± 1.1 kHz across localities.
Differences in RF were not associated with geographic distances.
A simple pairwise Mantel test revealed that there was no significant association between resting frequency differences and geographic distances in R. damarensis (Monte Carlo test observation: 0.084, 10,000 replicates, p = .5474).
The LME model selection suggested that the three most supported models all included AMT as an environmental predictor (Table 2), which together resulted in an accumulative weight (Cum.wt) of 77%.
In the order of highest ranking based on Akaike information criterion weight (highest AICwt), these three models were Reg + AMT + Sex, Reg + AMT + Sex + FA, and Lat + AMT + Sex. However, the evidence ratio indicated that the first model (Reg + AMT + Sex) was 3 and 9 times stronger than the second and third models, respectively ( Table 2). Each of the three variables in the most parsimonious model had a p-value of <.05 (Table 3).
All LMEs models indicated that forearm (FA) as a proxy for body size did not strongly influence resting frequency (RF) variation in R. damarensis. This was further supported by regression analysis which showed no correlation between FA and RF (R 2 = 0.02, p = .21).
RF of R. damarensis varied with AMT, exhibiting a positive correlation between RF and AMT (Figure 5a). Comparison of spatial (LAT and REG) and environmental (RH and AMT) factors using linear F I G U R E 3 Spatial autocorrelation due to strong spatial dependencies in observed data (a) and (b) a graph with no spatial autocorrelation after spatial dependencies had been accounted for by the simple random-effects model  mixed-effects models revealed that AMT and REG had the lowest AIC value. This, and the fact that LAT was not a significant predictor of RF, indicated that environmental discontinuities rather than continuous geographic divergence in RF strongly influenced resting frequency variation in R. damarensis. Prey detection range, AA, and RH were also not significant predictors of RF variation (Figure 5b-d). There was also T A B L E 2 Results from linear mixed-effects models (LMEs) testing for association between the resting frequency of Rhinolophus damarensis and environmental variables sexual dimorphism in the RF of R. damarensis, and females had higher RF than males (Figure 5e).
The inclusion of Reg in the model indicated that the northern and southern regional groups, based on genetic lineages uncovered by Jacobs et al., 2013; had distinct RFs; the southern region had higher RF than the northern region ( Figure 5f). However, the RFs of the two regions overlapped; the southern region had RFs which ranged from 82 to 89 kHz and that of the northern region ranged from 82 to 86 kHz.
Our separate linear mixed-effects models on just the southern populations revealed that the model with the highest ranking was RF ~ AMT + Sex. None of the other variables (PDR, AA and RH) in the lower ranked models had a p-value of <.05. Similar to the analyses which included all populations (Table 3), the association between AMT and RF was also positive (slope = 0.62279, p < .002).
The mean detection range and atmospheric attenuation of echolocation signals across populations of R. damarensis were 9.68 ± 0.63 m and 2.29 ± 0.19 dB/m, respectively. Across populations from the northern lineage, the mean atmospheric attenuation (2.43 ± 0.12 dB/m) was higher and mean detection range

| DISCUSSION
We found no evidence for James' Rule or for random genetic drift.
Body size was not correlated with RF nor climatic variables, suggesting that variation in RF was not the result of concomitant variation in body size as proposed by James' Rule. Similarly, the Mantel test showed no IbD effect and there was therefore no evidence that genetic drift was responsible for the variation in RFs. In contrast, the LMEs indicated that there was support for IbE in the form of an association between RF and region, which was based on two geographically separated genetic lineages. Furthermore, RF variation was also associated with the climatic variable, AMT across all populations and a similar association was found across the southern populations. This suggests that regional differences in RF were not simply due to lineage-specific differences but that local adaptation had exerted an influence. However, the relationship we found between RF and AMT is correlative rather than causative. Alternative evolutionary processes, for example, phenotypic plasticity, cannot therefore be excluded at this stage.
Most studies on geographic variation have focused on humidity as the main predictor of such variation (for a review see Jiang, Wu, & Feng, 2015) providing evidence that echolocating bats that are found in humid areas will experience severe atmospheric attenuation resulting in lower RF than those that are found in dry areas (Guillén et al., 2000;Heller & Von Helversen, 1989;Lawrence & Simmons, 1982).
For example, geographic variation in the RF of R. pusillus was positively correlated with mean annual rainfall (Jiang et al., 2010). However, it is known that the propagation of echolocation calls and acoustic signals, in general, are influenced by both temperature and relative humidity (Lawrence & Simmons, 1982;Luo et al., 2014) through their effect on atmospheric attenuation. Atmospheric attenuation is the result of a complex interaction between temperature, humidity, and the frequency of the sound being propagated (Luo et al., 2014). For example, in two species of horseshoe bats, R. simulator and R. swinnyi, distributed in the more mesic eastern half of southern Africa, acoustic divergence in RF was influenced by the interaction between temperature and humidity and the degree of influence was higher in R. swinnyi than in R. simulator because R. swinnyi echolocated at higher frequencies (R. swinnyi = 103.77 ± 1.70 kHz; R. simulator = 80.32 ± 2.20 kHz; Mutumi et al., 2016). However, our results for R. damarensis indicate that temperature (AMT) was the predominant climatic factor responsible for the divergence in echolocation RF. This was also the case for R. ferrumequinum (Jiang et al., 2015) in Asia and for R. capensis in southern Africa. Similar to R. damarensis, R. capensis has a distribution that extends into the more arid western and northwestern regions of South Africa (Neumann & Bamford, 2015;Odendaal et al., 2014) and its RF was also not correlated with RH (Odendaal et al., 2014).
This may be because both R. damarensis ( Figure 5b) and R. capensis (Odendaal et al., 2014) occupy more arid regions (Table 1; compare with Mutumi et al., 2016) and in such areas temperature is a better predictor of RF than RH.
The effect of temperature is mediated by its interaction with RH (see Luo et al., 2014;Mutumi et al., 2016), and this interaction may vary from one region to another such that one or the other may exert the predominant influence on atmospheric attenuation and the propagation of acoustic signals. Our results and those of Mutumi et al. (2016) suggest that in warm mesic regions, both temperature and relative humidity are likely to influence the propagation of acoustic signals, but in warm arid regions, the predominant climatic factor that defines atmospheric attenuation is more likely to be temperature. The complex nonlinear interaction between temperature and humidity may also explain the paradoxical results obtained by studies on the climatic influence on geographic variation in acoustic signals. For example, geographic variation in RF was positively correlated with mean annual rainfall in R. pusillus (Jiang et al., 2010) but negatively correlated with mean annual rainfall in Hipposideros ruber (Guillén et al., 2000). These contrary findings are often attributed to the complexity of natural selection (e.g., Jiang et al., 2015) but may also be a consequence of the nonlinear effects and potential collinearity of temperature and humidity on atmospheric attenuation (Luo et al., 2014). Neither Jiang et al. (2010) nor Guillén et al. (2000) considered the interactive effects of temperature and humidity on atmospheric attenuation.
Divergence in acoustic signals may arise from the action of different evolutionary forces, among which environmental selection is the most common (Podos & Warren, 2007;Wilkins, Seddon, & Safran, 2013). In the case of R. damarensis, divergence in RFs may be the result of selection for optimal detection ranges within the respective habitats of the populations in the two regions. Populations of R. damarensis in the northern region have RF on the lower end of the observed range for populations in the southern region. The lower RF of populations from the northern region could be the result of the more mesic conditions in northern Namibia. Acoustic divergence within species can arise when signals undergo selection for optimal propagation of acoustic signals in the local acoustic environment (Wiley & Richards, 1982) leading to populations occupying different ecological niches (Nosil, 2012).
Although we found significant differences in estimated PDRs between populations, there was no clear distinction between the PDRs of the two regions (Tables 1 and 4), suggesting that the ecological niches of the bats occupying these two regions, if different, may be defined by more than just atmospheric attenuation and its effect on PDR.
Nevertheless, the two regions have different climates. The northern part of Namibia is a more mesic region and characterized by woodland vegetation (Hoetzel, Dupont, & Wefer, 2015;Simmons, Griffin, Griffin, Marais, & Kolberg, 1998;White, 1983), while central Namibia and the Karoo ecoregions (Succulent and Nama) of South Africa are semi-arid (Neumann & Bamford, 2015;Okitsu, 2005;Thuiller et al., 2006; Figure 2). The importance of region as a predictor of RF variation may be indicative of these environmental discontinuities.
If so, the coincidence of phenotypic difference with environmental discontinuities strongly suggests that divergence in this species is the result of IbE possibly facilitated by adaptation of sensory signals to local environmental conditions. Such environmental discontinuities in the context of the propagation of acoustic signals may be partly defined by atmospheric conditions such as temperature and humidity (this study; Jiang et al., 2015;Mutumi et al., 2016;) but in all probability also by other unobserved environmental factors (i.e., latent variables). The association of latitude and altitude with RF variation in two species of horseshoe bats (Mutumi et al., 2016) and with region in R. damarensis suggests that in addition to temperature and humidity, other as-yet-unidentified environmental variables influence RF and this too may be suggestive of the influence of environmental discontinuities on phenotypic traits.
Geographic variation in RF between populations of R. damarensis in the two regions may have evolved as a consequence of adaptation to local environmental conditions. Such local adaptation could impede gene flow between the northern and southern regions occupied by R. damarensis (Gillam & McCracken, 2007;Rundle & Nosil, 2005) through several processes including selection against migrants (Hendry, 2004;Thibert-Plante & Hendry, 2009) (Fitzpatrick, 2012;Pfennig et al., 2010). One potential limitation to plasticity is the tight coupling between RF and the acoustic fovea in high duty cycle bats, like R. damarensis, that use Doppler shift compensation (Neuweiler, 1984;Schnitzler & Denzinger, 2011). However, HDC bats are able to shift their RFs (up to 3.9 kHz) in response to both neighboring conspecifics and different ambient noise conditions (Hage, Jiang, Berquist, Feng, & Metzner, 2013;Hiryu et al., 2006).
Although small, such shifts in frequency encompass the differences between the mean RFs of neighboring populations of R. damarensis and may allow bats to optimize their detection range under the local climatic conditions without a change in allele frequency. However, an adequate test of phenotypic plasticity requires greater knowledge of the genetic basis and heritability of phenotypic traits, of gene flow between populations and of interhabitat differences in foraging behavior of these bats.
Sexual dimorphism within a species may suggest that social rather than ecological selection is responsible for any divergence. Sex in R. damarensis was a significant predictor of RF with females generally echolocating at slightly higher RF than males ( Figure. 5d). This is also true for several other horseshoe bat species (Chen, Jones, & Rossiter, 2009;Siemers et al., 2005;Yoshino et al., 2006). Bats exhibit a wide diversity of mating systems, but little is known about the mechanisms involved in courtship and mating (Grilliot, Burnett, & Mendonc, 2009).
Recent evidence suggests that echolocation calls may be a source of information for individual recognition (Siemers & Kerth, 2006;Siemers et al., 2005). Sexual dimorphism in RF may have an essential social function by signaling the sex of the caller (Kazial & Masters, 2004), promoting mate recognition (Guillén et al., 2000) and reproductive success (Grilliot et al., 2009). It has been shown for some species that call frequencies of males and females do not overlap, and frequency may encode the sex of the signaler reliably (Neuweiler et al., 1987 Recent studies suggest that climate change may have a severe negative impact on the sensory ecology of sound-mediated behaviors (Luo et al., 2014). An increase in temperature and/or relative humidity, combined with the high RF used by R. damerensis, is likely to result in a decrease in the detection volume of R. damarensis under a global warming scenario (Luo et al., 2014). Global warming may therefore directly impact the sensory ecology of organisms reliant on sound-mediated behaviors (Luo et al., 2014). If changing ambient temperature increases atmospheric attenuation, then echolocating species of bats may have reduced prey detection distances (Luo et al., 2014). Reduction in prey detection distances could lead to ineffective foraging and an increase in the costsbenefits ratio of foraging, which could lead to local extinction of populations.
Climate-induced variation in acoustic signals may not be restricted to bats. There is also some support for climate-induced geographic variation in bird song (Snell-Rood, 2012), suggesting that lineage diversification in general may be driven by habitat-mediated differences in acoustic signals, especially when such differences are accompanied by genetic differentiation. Climatic gradients in conjunction with other ecological discontinuities could lead to local adaptation and/or phenotypic plasticity with directed dispersal (Edelaar et al., 2008) such that gene flow is restricted, allowing lineages to diverge.