Global patterns of functional trait variation along aridity gradients in bats

Aim: Our understanding of the biological strategies employed by species to cope with challenges posed by aridity is still limited. Despite being sensitive to water loss, bats successfully inhabit a wide range of arid lands. We here investigated how functional traits of bat assemblages vary along the global aridity gradient to identify traits that favour their persistence in arid environments. Bats. Methods: We mapped the assemblage- level averages of four key bat traits describing wing morphology, echolocation and body size, based on a grid of 100- km resolution and a pool of 915 bat species, and modelled them against aridity values. To support our results, we conducted analyses also at the species level to control for phylogenetic autocorrelation. Results: At the assemblage level, we detected a rise in values of aspect ratio, wing loading and forearm length, and a decrease in echolocation frequency with increasing aridity. These patterns were consistent with trends detected at the species level for all traits. Main conclusions: Our findings show that trait variation in bats is associated with the aridity gradient and suggest that greater mobility and larger body size are advantageous features in arid environments. Greater mobility favours bats’ ability to track patchy and temporary resources, while the reduced surface- to- volume ratio associated with a larger body size is likely to reduce water stress by limiting cutaneous


| INTRODUC TI ON
Arid and semi-arid regions cover c. one third of the planet's land surface, encompassing key habitats for a wide array of species (Safriel et al., 2005). Additionally, the extent of such regions is predicted to increase by 11-23% by the end of the century, thus affecting a greater number of biological communities (Huang et al., 2016). Arid environments are characterized by water limitations and unpredictability in precipitation, which result in low primary productivity and large spatio-temporal variation in food resources (Maestre et al., 2016). These features, often in combination with high temperatures and other abiotic attributes associated with aridity, represent major stressors for local species, shaping assemblage composition (O'Connell & Hallett, 2019;Rymer et al., 2016;Schwimmer & Haim, 2009). A range of behavioural, morphological and physiological strategies and adaptations are employed by species to cope with these stressors. For example, vertebrates may reduce water stress in several manners, e.g. via increased skin resistance to cutaneous evaporative water loss, burrowing habits and high mobility to access scattered resources (Fuller et al., 2014;Oliveira & Scheffers, 2019;Rymer et al., 2016;Williams & Tieleman, 2005). These adaptations have been mostly studied in species inhabiting deserts, and we still have a poor understanding of which traits favour the presence of species with increasing aridity across broad geographical extents.
Over large geographical scales, average functional traits of species assemblages are associated with environmental conditions along gradients as a result of the combination of environmental filtering and ecological adaptation (Barnagaud et al., 2017;Holt et al., 2018;Lawing et al., 2017;Violle et al., 2014). Describing variation of species traits along aridity gradients can, therefore, provide insights on which features contribute to the ability of species to withstand arid conditions (Rymer et al., 2016).
Bats are a diverse and ubiquitous order of mammals, being found across all continents, except the polar regions. Their widespread distribution is accompanied by a great morphological and ecological diversification, with the ability to occupy a wide range of environmental conditions, from temperate and tropical rain forests to hyper-arid deserts (Simmons, 2003). Bats are considered particularly sensitive to water loss due to the presence of naked wing membranes (Minnaar et al., 2014) and to the exposure to air convection during flight (Carpenter, 1986). Nonetheless, they are a successful group in arid environments (Carpenter, 1969). The nocturnal habits of bats represent an advantage to face aridity constraints, and species-specific adaptations to extreme aridity, such as reduced active thermoregulation and modified kidney structure, further contribute to their high diversity in many arid ecosystems (Bondarenco et al., 2014;Carpenter, 1969;Munoz-Garcia et al., 2016). Thanks to these features, bats constitute a good candidate group for investigating functional traits' responses to increasing aridity.
Along the global spectrum of moist to arid environmental conditions inhabited by bats (hereafter 'aridity gradient'), a progressive reduction in primary productivity and vegetation cover occurs, alongside the increased pressure of the stressors mentioned above (Safriel et al., 2005). In this context, bats and other wildlife are gradually exposed to structurally more open environments. In echolocating bats, the primary functional traits that mediate their interaction with the structure of the environment relate to features of wing morphology and echolocation calls (Aldridge & Rautenbach, 1987).
These features are associated, respectively, with flight abilities and sensory perception in relation to background objects. A high wing loading (high weight to wing surface ratio) confers a fast but poorly manoeuvrable flight, while high values of aspect ratio (narrow wings) are associated with energy-efficient flight (Norberg & Rayner, 1987).
When it comes to echolocation, low frequency calls are less subject to atmospheric attenuation and can travel farther, being thus well suited to open spaces, whereas high frequency calls provide more detailed information about nearby objects and are associated with more cluttered environments (Denzinger & Schnitzler, 2013).
Species foraging in open spaces (open space foragers) are typically adapted to fast flight, with high values of wing loading and aspect ratio, and tend to exhibit low echolocation frequencies (Norberg & Rayner, 1987). In contrast, bats flying within the vegetation (narrow space foragers) are characterized by manoeuvrable flight with low values of wing loading and aspect ratio and high echolocation frequencies (Denzinger & Schnitzler, 2013). With trait values covering a continuum between these two extremes, morphological and echolocation trait variation has been applied to identify the effects of an increase in open space following forest fragmentation, reporting the partial loss of narrow space foragers from deforested areas (Estrada-Villegas et al., 2010;Farneda et al., 2015;Hayes et al., 2007).
Similarly, the progressive decrease in the structural complexity of vegetation with aridity might drive a shift in assemblage traits towards values associated with greater adaptation to open spaces.
The increase in aridity itself also exerts a greater water stress on species. Several adaptations to preserve water have been identified in bats, such as modified kidneys to concentrate urine, or torpor, which contributes to reducing respiratory evaporative water losses (Carpenter, 1969;Munoz-Garcia et al., 2016). Cutaneous evaporative evaporation. These findings highlight the importance of extending attention from species-specific adaptations to broad scale and multispecies variation in traits when investigating the ability of species to withstand arid conditions.

K E Y W O R D S
aridity gradient, assemblage level, bats, body size, echolocation, functional trait variation, global patterns, wing morphology water loss, one of the significant routes of water loss in bats, is dependent on body surface area (Minnaar et al., 2014). It is possible that larger body sizes, characterized by a lower surface-to-volume ratio, are comparatively less exposed to water loss and thus favoured as aridity increases (Munoz-Garcia et al., 2016). This could be particularly true for bats, since larger bodied bats have relatively smaller wings when compared to smaller species (Findley et al., 1972), thus reducing exposure via the naked membranes of the wings. Alternatively, since aridity is often associated with warm temperatures, smaller sizes could be associated with high aridity as a feature favouring heat exchange, as stated by Bergmann's rule (Bergmann, 1847;Blackburn & Gaston, 1996;sensu Meiri, 2011;Torres-Romero et al., 2016).
However, several studies suggested that this might not be generalizable across all taxonomic groups and geographical scales (Medina et al., 2007;Rodríguez et al., 2008). In echolocating bats (hereafter only 'bats'), the overall small body size across species might allow for sufficient heat dissipation in hot environments, while water retention might be a more critical constraint than temperature in highly arid regions. However, trends in traits identified along aridity and temperature gradients are likely not entirely comparable as a number of arid areas are located in cold regions (e.g. cold deserts).
In this study, we aimed to investigate the set of functional traits that enable bats to cope with arid conditions. We examined how a bat assemblage's measures of wing morphology (aspect ratio and wing loading), echolocation frequency and body size vary along the aridity gradient. Due to the decreasing vegetation structural complexity, we hypothesized that increasing aridity will be associated with bat assemblages with high wing loading and aspect ratio and low echolocation frequency. Due to the need to mitigate evaporative water loss, we also hypothesized that aridity will select for assemblages of bats with larger body sizes. To test our hypotheses, we conducted analyses at the assemblage level (with an assemblage defined as the set of echolocating bat species coexisting in the same grid cell). While a grid-based approach not accounting for phylogenetic relationships among species is appropriate to test ecological hypotheses (e.g. environmental filtering, de Bello et al. 2015), the inclusion of models incorporating phylogenetic signal can provide further information and disclose evolutionary patterns (de Bello et al. 2015). Therefore, we complemented our analysis by testing our hypotheses also at the species level via a phylogenetic comparative approach (Figure 1, see e.g. Olalla-Tárraga et al., 2010). In addition, this species-level F I G U R E 1 Diagram summarizing the main methodological steps. FMAXE = frequency of maximum energy; GAM = generalized additive model; CWM = community weighted mean; PGLS = phylogenetic generalized least squares EchoBank, , as well as with targeted literature searches to reduce gaps in the dataset. A list of the data sources is found in the Appendix. We collected information on the main traits of interest, namely forearm length (as a proxy for body size), wing loading, aspect ratio and frequency of maximum energy (FMAXE) of the echolocation call ( Figure 1a). Since it is recommended to impute data gaps using a larger set of correlated traits than those of interest to the study (Van Buuren, 2018), we also collated data on body mass and echolocation call duration. Wing loading and aspect ratio were defined following Norberg and Rayner (1987; see Supporting Information Appendix S1 for further details), while FMAXE was obtained for the harmonic with maximum energy (Monadjem et al., 2010). When collecting data from secondary sources (large databases) we discarded data obtained via extrapolation or other indirect methods.

| Imputations
We imputed data gaps using a multiple imputation technique (multi- . We further checked imputation accuracy by introducing an additional artificial 5% of missing data to the original data and by verifying that imputed data approximates observed data. To compare observed and imputed data we employed the percent bias, which measures the average tendency of the imputed values to be larger or smaller than the observed ones (Zambrano-Bigiarini, 2017). For satisfactory performance, the upper accepted limit for the percent bias is 5% (Demirtas et al., 2008;Van Buuren, 2018). We conducted this for 10 randomly selected trees of the ones in use, and simulated data were compared to observed data after obtaining the median across the 25 imputed datasets. All 2,500 datasets were retained and used to assess the influence of imputation uncertainty on the results of both the assemblage-and species-level analyses.

| Spatial data processing
To estimate assemblage-level average traits, we generated a global grid of 100 km × 100 km resolution using the equal-area Mollweide projection, and for each cell determined species presence and cell areas covered using IUCN range maps (Figure 1d). This resolution is often considered a good compromise to investigate broad scale geographical patterns while accounting for the intrinsic uncertainty in the range maps (Buckley & Jetz, 2008;Santini et al., 2017). For each imputed dataset, we then calculated the assemblage mean of each trait per cell, where each species value is weighted by the area of the cell that its range covers (often referred to as community weighted mean, hereafter only 'CWM'). CWMs were calculated using the function 'ddply' from the R package 'plyr' (Wickham, 2011). We ex-  Figure S3.9), which was rescaled at the 100 km × 100 km resolution in order to assign a value of aridity to each assemblage. It is important to note that the Global Aridity Index (hereafter only 'AI') is expressed as the ratio between mean annual precipitation and mean annual evapotranspiration (Trabucco & Zomer, 2019). Therefore, the greater is the AI, the less arid is the environment. Since our focus is on the influence of aridity on traits, we will discuss trends referring to 'increase in aridity' and 'decrease in aridity', which correspond, respectively, to smaller and greater AI values. Data for the species-level analysis were obtained by extracting the value of the 5th percentile of the AI (high aridity) within the species geographical ranges (Figure 1d). This allowed us to focus on the putative highest aridity level experienced by each species in its distribution.

| Assemblage-level analysis
To examine variation in community functional traits along the global aridity gradient while controlling for spatial autocorrelation using a trend surface analysis, we employed generalized additive models (GAMs, Wood, 2017). While aridity was included as a linear predictor, spatial patterns were modelled by including a two-dimensional smoother of the coordinates of the centroid of each cell (see e.g. Fletcher & Fortin, 2018;Safi et al., 2011). This also accounted for other unobserved variables potentially affecting trait variation spatially. We verified the efficacy of the GAMs in controlling the spatial autocorrelation of the data with the use of variograms (see Supporting Information Figure S4.10). Models were built independently for each trait of interest and they were run across all datasets produced via the imputations, for a total of 10,000 models, 2,500 for each trait (Figure 1e). To exclude potential outliers, each model was re-run after excluding the 0.1% most influential observations as identified by the function 'influence.gam' from the 'mgcv' R package (Wood, 2011). The models were specified as follows: where CWM represents the community weighted mean for a certain trait, AI is the AI value for the assemblage and s(X, Y) specifies a twodimensional spline smoother for the spatial position of the cells. Before fitting the models, the AI was raised to the power of 0.4 to improve normality. We chose a 0.4 power transformation as it proved to be the most effective in improving normality of the AI. GAMs were fit using the function 'gam' (Wood, 2011). We tested the significance of aridity for modelling each trait's CWM by comparing the full models as described above to null models retaining only the spatial component.
If the difference in the Akaike information criterion (ΔAIC) between the null and full model was ≥ 2, aridity was considered significant in explaining trait variation. We repeated all validation procedures for 25 randomly selected datasets out of 2,500. Results across models were summarized for each trait by providing medians and 95% confidence intervals of models' parameters.
To ensure that the trends observed were not an artefact of unincluded factors, we calculated the Pearson's r for the correlations between the AI and other environmental variables that could potentially affect CWM, namely annual mean temperature, temperature seasonality and precipitation seasonality. We obtained these environmental layers from WorldClim version 2.1 (Fick & Hijmans, 2017) and processed them analogously to the aridity layer. In particular, we projected and rescaled them at the 100 km × 100 km resolution, extracted values for the cells included in the analysis and then, calculated Pearson's r of the correlations with the AI.

| Species-level analysis
As an independent test of the robustness of our results, we also tested our hypotheses using species as the unit of the analysis and accounting for phylogenetic autocorrelation. We employed phylogenetic generalized least squares (PGLS) models to assess the variation of traits across species as a function of the aridity level tolerated by each species (Figure 1f). Using the 5th quantile instead of the mean or the median of the AI value allowed us to focus on the putative extreme aridity experienced by the species within their current distribution while excluding possible errors due to the coarseness of these maps. We also tested weather using the 2.5th quantile or the median would result in qualitatively different results (see Supporting   Information Table S4.3). To facilitate model convergence, trait data were log-transformed and the AI was raised to the power of 0.4. Analogously to the assemblage-level analysis, we produced a total of 10,000 PGLS models, 2,500 for each trait, and each model was re-run after excluding model-specific outliers. Outliers were identified as observations with studentized residuals > ±3, following Jones and Purvis (1997). The models were fitted via the function 'gls' from the R package 'nlme' (Pinheiro et al., 2019) by specifying a Pagel's λ correlation structure available from the package 'ape' (Pagel, 1999;Paradis & Schliep, 2019), and they included 904 species for which phylogenetic information was available (Figure 1f). Results across models were summarized for each trait by providing medians and 95% confidence intervals of models' parameters. We provide a CWM ∼ AI + s (X, Y) measure of the phylogenetic signal detected for the relationship of each trait with aridity via the Pagel's λ as estimated by the models.
All analyses were conducted in the R environment, version 3.6.1 (R Core Team, 2019). for diagnostic graphs for a sample tree). The comparison between imputed and observed data for artificially added missing records always produced a percent bias below the threshold of 5% (Supporting Information Table S2.1), confirming the suitability of multiple imputation by chained equation to fill in data gaps in our dataset. From the geographical representation of the CWMs, we can identify major patterns in trait variation and regions retaining extreme values (Figure 2a-d). Globally, aspect ratio and wing loading peak within arid regions (Figure 2a,b, see also Supporting Information Figure S3.9 for representation of aridity classes). The highest FMAXE are observed in the humid Southeast Asia, with some high values also located in the arid Arabian Peninsula, while a low peak is observed both in the arid zone of Australia and South America and in the Nearctic and Palaearctic zones (Figure 2c).

| Assemblage-level analysis
Forearm length, instead, is greater in the semi-arid and arid regions of India and in the arid and hyper-arid zones surrounding the Sahara, with low values widespread in the Nearctic and Palaearctic zones (Figure 2d).
Model selection indicated that the models including the AI (full models) performed substantially better than null models in explaining the trait data (ΔAIC always well above 2, see Supporting Information Table S4.2). As predicted, aspect ratio and wing loading (addressing mobility) showed a positive relationship with increasing aridity [across 2,500 models for each trait, aspect ratio: median β aridity (95% Confidence Intervals (CIs)

| D ISCUSS I ON
In this study, we provide novel insights into the functional responses of bats to aridity. In accordance with our hypotheses, we observed that increased aridity is associated with wing morphology and echolocation attributes suitable to move and forage in open spaces (higher aspect ratio and wing loading and lower FMAXE values) and with larger body sizes (greater forearm length). Despite the moderate levels of phylogenetic conservatism for the traits considered, these patterns were consistent between the two different methodological approaches, indicating that our results are valid independently from the phylogenetic relationships among species.

| Patterns of wing and echolocation traits along the aridity gradient
As aridity increases, assemblages are increasingly characterized by long narrow wings and low echolocation frequencies. These trends reflect the changes in vegetation structure associated with aridity and likely result from environmental filters favouring species better adapted to structurally more open habitats (Salinas-Ramos et al., 2020). The prevalence of these features in deserts and semideserts is likely linked to the relaxation of the environmental pressures typically associated with flight in highly cluttered spaces.
Our results suggest that, as openness of the habitat increases, manoeuvrability requirements to avoid background obstacles are reduced, while speed and cost-efficiency of flight to cover longer distances are favoured. Accordingly, there is an increased need to receive information about prey and physical obstacles further away in space to allow turning, and this is made possible by the use of lower echolocation frequencies (Denzinger & Schnitzler, 2013). These trends were confirmed by the species-level approach.
However, the trend for FMAXE received less support, potentially due to the complex and plastic nature of echolocation that is poorly represented by a single parameter (Denzinger & Schnitzler, 2013).
As arid lands are expected to expand under climate change (Huang et al., 2016), the trends observed in this study might gain further importance in the future. In fact, previous studies on bat functional traits have associated low values of wing aspect ratio with a higher extinction risk , likely due to the limited dispersal ability and smaller range sizes that are associated with low values for this trait  2019, but F I G U R E 3 Summary of the results of (a) the assemblage-level analysis and (b) the species-level analysis. (a) Distribution (median and 95% confidence intervals) of the coefficients for the aridity index (AI) across the 2,500 models per trait. (b) Distribution (median and 95% confidence intervals) of the coefficients for the 5th percentile of aridity estimated across the 2,500 models per trait. Lower values of the AI correspond to higher aridity levels, and therefore for example, a negative coefficient indicates an increase in trait values with increasing aridity. The detected signal is considered highly probable if the 95% confidence intervals do not include zero, indicating that effects were consistent across models. FMAXE = frequency of maximum energy see also Varzinczak et al., 2020). The expansion of arid lands could result in more open habitats and, according to our findings, this could create less suitable conditions for clutter foragers and thus potentially increase their extinction risk.
In addition to a greater suitability to open spaces, greater mobility can also be advantageous for other aspects of bat ecology in arid environments. Accessing scattered water bodies and ephemeral and patchy food sources, conditions typical of desert and semi-desert habitats, likely leads to increased daily movements, as shown for example, for the desert-adapted species Rhinopoma microphyllum (Bell, 1980;Conenna et al., 2019;Egert-Berg et al., 2018;Korine & Pinshow, 2004). Higher wing loading and aspect ratio, as observed in our study for bat assemblages in more arid areas, allow for fastand energy-efficient flights, representing an advantage when facing long nightly commutes (Norberg & Rayner, 1987 (Norberg & Rayner, 1987). All the regions mentioned above also hold corresponding trends in FMAXE, with progressively lower values as aridity increases, supporting the idea of these assemblages being dominated by open space species.
However, our results represent a simplification of complex patterns, and traits might respond to thresholds in aridity in a nonlinear matter, or diverge from main trends under extreme conditions also due to low species richness. For example, we also identify arid areas that display a prevalence of clutter-type values across both wing and echolocation traits, for example, the Arabian Peninsula ( Figure 2). It is possible that, at extremely low environmental productivity, prey is occasionally more tightly associated with vegetation patches, favouring species capable of foraging closer to the vegetation or gleaning prey from plants and the ground. Indeed, several desert-adapted species display such hunting strategies (e.g. Plecotus christii, Asellia tridens, Anthrozous pallidus; Bell, 1982;Korine & Pinshow, 2004). Further research is needed to explain the prevalence of narrow space foragers in some of the highest-aridity areas.

| Patterns of body size along the aridity gradient
The body size in bat assemblages increased with aridity. This pattern was strongly supported by the phylogenetic-controlled approach, indicating that species larger than expected by their phylogenetic position were more frequent in drier environments. Our results suggest that, given the overall small sizes and susceptibility of bats to cutaneous evaporative water loss, it might be advantageous to retain a large body size due to its relative lower surface-to-volume ratio. This is in accordance with James (1970), who stated that, when considering the gradient from humid to dry, larger sizes are expected to be associated with drier conditions. This trend might gain more relevance when seen in conjunction with the increase, at higher aridity, in wing loading, signalling a decrease in wing surface relative to body size.
Similarly, Olalla-Tárraga et al. (2009) showed that anuran body size positively correlated with water deficit in the Neotropical Cerrado biome, and labelled this pattern as responding to the 'water availability hypothesis', according to which anurans from drier areas exhibit larger body sizes to reduce water loss and prevent desiccation.
These patterns have also been previously associated with the 'water conservation hypothesis', which states that species are expected to be larger in more arid environments compared to mesic ones due to the advantage conferred by a larger body size in reducing desiccation rates (Nevo, 1973). While the relationship between bat body size and aridity has received limited attention in bats, the relationship with temperature has been extensively investigated, with various outcomes either confirming, refuting or contrasting Bergmann's rule (e.g. intraspecific: Jiang et al., 2019;Penone et al., 2018;interspecific: Safi et al., 2013). The divergence from Bergmann's rule in endotherms has been associated with size-independent adaptations to cold or warm climates (Medina et al., 2007), such as the ability to hibernate in bats. However, an interesting future avenue would be investigating whether these various identified trends between body size and temperature arise because of the combination of warm temperatures with various aridity levels, thus determining different outcomes depending on whether arid conditions are associated with either cold or warm temperatures (James, 1970).

| Limitations of the study
The results of this study are consistent across the different methodological approaches employed. Nevertheless, we would like to draw the reader's attention towards potential shortcomings of the study.
First, the effects of aridity on bat traits detected across the analyses are relatively weak. Many other factors, including other environmental variables, historical and local components are likely to combine to determine the trait values observed. Second, Hawkins et al. (2017) have recently highlighted that the repeated co-occurrence of species across grid cells, as is the case for the assemblage-level analysis in this study, can potentially increase the risk of Type I errors. While we encourage the reader to consider these caveats when evaluating the results of this study, we believe the validity of the study is not hindered since our aim was that of detecting overall effects of aridity on traits, rather than predicting trait values. Additionally, the results of the species-level analysis, which is insensitive to the bias caused by repeated species co-occurrences, support the overall robustness of conclusions.

| CON CLUS IONS
Taken together, our analyses revealed patterns of variation along the global aridity gradient in several important bat traits. Importantly, as aridity increases, we observed a prevalence of traits favourable to hunting in open spaces, as well as larger body sizes. These features might be advantageous in dry climates, allowing bats to exploit habitats with open structure while coping with scattered resources, as well as reduce evaporative water loss.
Environmental aridity poses a multitude of challenges to species and can be a driver of changes in local assemblage composition.
Describing variation in functional traits along the aridity gradient, as presented in this study, can help us identify these changes and suggest how different groups of species might respond to the expected expansion in arid lands.

ACK N OWLED G M ENTS
The authors thank Enrico Bernard for contributing field data on wing morphology and Elina Numminen for insightful methodological discussions. IC was supported by the Ella and Georg Ehrnrooth Foundation (grants #8-6349-18 and #9-7968-5) and by the Nordenskiöld-Samfundet. RR was supported by ARDITI-Madeira's Regional Agency for the Development of Research, Technology and Innovation (grant M1420-09-5369-FSE-000002).

AUTH O R CO NTR I B UTI O N S
IC, LS and DR conceived the ideas and designed methodology; IC collected the data; IC and LS analysed the data and MC, RR and AM provided important input into data interpretation. IC led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

DATA AVA I L A B I L I T Y S TAT E M E N T
The