Listening in the dark: acoustics indices reveal bat species diversity in a tropical savannah

ABSTRACT Surveying biodiversity using bioacoustics has become increasingly common worldwide, although it is mostly concentrated in temperate regions. The variety of automatic recorders, the development of free analytical tools, and several acoustic indices have increased the number of studies worldwide. The bioacoustic approach is essential for application in poorly surveyed regions with the pressure of human activities, such as the Brazilian cerrado. We tested the association of four bat diversity metrics (taxonomic, phylogenetic, and two functional diversity metrics, being one based on morphological and the other on acoustical traits, with five commonly used acoustic indices. We used a dataset of 608.4 h obtained from 30 sampling points in three protected areas in Central Brazil. Using Flexible Discriminant Analysis, we identified 21 bat species used in our subsequent analysis. The Entropy index was the best predictor of taxonomic and phylogenetic diversity, whereas the Acoustic Complexity Index was the best predictor of functional morphological diversity. We concluded that acoustic indices are suitable for estimating the diversity of insectivorous bats in the cerrado. However, we registered only part of the bat community, and bats can vary seasonally masking the real diversity of the study area; thus, this method should be used parsimoniously.


Introduction
Measuring biodiversity is a difficult task, and more than a single metric approach is necessary to characterise it (Willis and Whittaker 2002). Biodiversity can be described in many ways, and three of them include the taxonomic, phylogenetic, and functional dimensions. Such dimensions are essential for capturing the community aspects that cannot be perceived using only species richness (Cadotte et al. 2011;Srivastava et al. 2012;Cisneros et al. 2015). Measurement of taxonomic diversity (species richness, evenness, and diversity) is the most widely used method, which does not show species differences in their ecological roles (Cisneros et al. 2015). Functional diversity has been applied to include species' ecological attributes associated with ecosystem functioning through functional traits (Petchey and Gaston 2006). Functional traits are defined as morphological, physiological, behavioral, or phenological characteristics that reflect the ecological role of a species in an ecosystem (Nock et al. 2016). Although there are indications that phylogenetic diversity tends to correlate with functional diversity , measuring diversity through the relationships between species is a way to incorporate evolutionary aspects in the structuring of communities (Faith 1992;Cadotte et al. 2011;Flynn et al. 2011;Srivastava et al. 2012).
The challenge of calculating the diversity of different animal species is the intrinsic difficulty of knowing which species are present in a specific area. Monitoring biodiversity has inherent challenges for each taxonomic group that depends on data collected from limited spatial and temporal scales (Yoccoz et al. 2001). In general, there are methodological problems related to the effort needed to carry out good inventories, correct identification of species, expenses for personnel, field equipment, and displacement. In countries like Brazil, considering bats, estimates show that it would take around 200 years for a good quality inventory across the country (Bernard et al. 2010). Time is a luxury in regions with vast and constant human pressure, such as the Brazilian cerrado hotspot, where the loss of connectivity is even more critical than habitat loss per se (Grande et al. 2020). One possibility to overcome such difficulties is to apply bioacoustics based on an automated digital recording system, which can simultaneously register vocalising animals in any area or habitat and create a permanent voucher material for the surveyed localities (Sueur et al. 2008a). Automatic recognisers of bat calls can only identify a small number of species (primarily European) (Walters et al. 2012) and cannot be used in the cerrado (Arias-Aguilar et al. 2018;Hintze et al. 2020).
The field effort for biodiversity surveying is considerably smaller than that when using traditional methods for bat inventories, such as mist-nets for capturing species. The amount of data generated by recordings are typically large. The use of acoustic indices (AIs) is a promising approach for overcoming this problem. AIs are mathematical representations of the acoustic diversity of a community that can be applied to quickly characterise and compare areas for biodiversity assessments over time (Sueur et al. 2014). These indices can reflect alpha diversity (Machado et al. 2017) or beta diversity when the intention is to compare how acoustically similar (or dissimilar) two or more areas are (Mammides et al. 2017;Hayashi et al. 2020). The AIs approach is promising, mainly for registering acoustic communities composed of species that are difficult to capture with mist nets, such as aerial and gleaning insectivores that forage in clutter and edge spaces. Such an approach is even critical in regions where quick decisions must be made to prevent environmental impacts and promote biodiversity conservation. However, it is important to know what acoustic data tell us about the diversity of a given taxonomic group, or the association level of acoustic diversity and its biological diversity.
Several studies have been conducted in Brazil to evaluate the potential application of AIs for biodiversity studies (Machado et al. 2017;Duarte et al. 2021;Campos et al. 2021). However, no research was conducted on the application of AIs for bats. We believe that it is worthwhile to apply the indices to study bats, as they represent a significant part of the mammalian fauna in the Neotropics (Simmons and Voss 1998;Aguiar et al. 2020) and a key group for ecosystem functioning (Muscarella and Fleming 2007;Aguiar et al. 2021). Therefore, we aimed to verify whether AIs reflect the different dimensions of the diversity of aerial insectivorous bats in the cerrado. To do this, we tested the following hypotheses: (1) AIs will be positively related to bat diversity (taxonomic, phylogenetic, and functional) in the cerrado, that is, the higher the diversity, the higher the index value, (2) Indices that have frequency measurement adjustments, such as Acoustic Diversity Index (ADI) and Acoustic Evenness Index (AEI), respond better to the diversity of cerrado bats in larger frequency bandwidths because insectivorous bats present a wide range of frequencies; and (3) indices such as the Acoustic Complexity Index (ACI), which allows adjusting the size of the temporal bands, respond better to the diversity of cerrado bats in bands with a shorter duration because bats usually emit several calls of short duration. Data from the literature show for several regions a positive relationship between species diversity and AIs (Machado et al. 2017;Mammides et al. 2017;Jorge et al. 2018;Retamosa et al. 2018), that is, areas with high species diversity tend to present high index values. Therefore, we expect to identify which biodiversity aspects can be captured by the indices and how feasible they are when applied to bat acoustic monitoring because AIs are primarily developed for species with low frequency and distinct temporal structures (e.g., birds and frogs).

Materials and methods
We conducted the study in three protected areas in the cerrado biome of central Brazil  S − 47°43' Long W), and the Environmental Protection Area Gama e Cabeça de Veado (15°55' Lat S − 47°52' Long W). Details about each studied area are provided in the Supplementary Information (S1). The cerrado biome is a Neotropical savannah-like vegetation characterised by a mosaic of habitats ranging from evergreen forests along rivers and streams, dry forests, scrublands, and open fields (Eiten 1972;Ratter et al. 1997). The climate is seasonally marked by a rainy season from October to March, and a dry season from April to September. The average annual rainfall is 1500 mm, with mean temperatures ranging from 22°C to 27°C (INMET, https://tempo.inmet.gov.br).

Data collection
A total of 30 sample points were set, with 10 points in each protected area. Bat calls were recorded with Song Meter SM2Bat (Wildlife Acoustics, Maynard, MA, USA) ultrasonic detectors (referred to as Autonomous Recording Unit or simply ARU) in wave format with a sample rate of 384 kHz and 16 bits. We installed one ARU per sample point at each studied site, with a minimum distance of 1 km between them. Recordings were conducted for three minutes at two-minute intervals, covering 12 h from sunset to sunrise, always on new moon nights. Each site was sampled for two nights, with an interval of one to three months between sample nights. The samples were collected between February and September 2014, with 480 recording hours per study site, totaling 1,440 h or 17,280 files. Such efforts are efficient for surveying a significant portion of the existing local biodiversity of bats (Pieretti et al. 2015;López-Baucells et al. 2021) and birds (Metcalf et al. 2021).

Bat species identification
We analyzed the sound files with Avisoft-SASLab Pro version 5.2 (Avisoft Bioacoustics) software, with sonograms and spectrograms constructed from 1,024 Fast Fourier Transform (FFT) points and a Hamming window function with 93.75% overlap. We considered a bat present at a site when at least one pass, defined by us as a sequence of three or more sequential pulses in the searching phase (echolocation calls that bats produce while searching for food). We identified the passes by spectrogram analysis, measuring acoustic traits such as the start frequency (SF), end frequency (EF), frequency of maximum energy (FME), call duration (CD), and pulse interval (PI) (Supplementary Information, Table S2). Species calls characteristics were identified according to the methods described by Arias-Aguilar et al. (2018), Barataud et al. (2013), and Jung et al. (2014). After species identification, we conducted Flexible Discriminant Analysis (FDA). The FDA tests whether predetermined bat species can be correctly classified according to selected metrics (the acoustic traits, in our case). The analysis can be performed without assuming data normality (Hastie et al. 1994). As Russo and Jones (2002) suggested, we used a discrimination threshold of 0.8 to consider the passes satisfactorily classified as the species previously identified and used them in later analysis, where we quantified the diversity metrics response to AIs. Molossus currentium and Pteronotus gymnonotus were recorded only once and were not included in the discriminant analysis. However, their spectral characteristics, such as pulse shape and frequency bandwidth, allow them to be safely identified. FDA was performed using the FDA function of the mda package version 0.5-2 (Hastie and Tibshirani 2020) in R software (R Core Team 2020).

Acoustic indices
We divided each three-minute file into a one-minute file using the Kaleidoscope software (Wildlife Acoustics). We used five of the most commonly used acoustic diversity indices related to alpha diversity (Sueur et al. 2014;Bradfer-Lawrence et al. 2019) to analyze the data: the Bioacoustic Index (BI) (Boelman et al. 2007), Acoustic Entropy Index (H) (Sueur et al. 2008b), Acoustic Diversity Index (ADI), Acoustic Evenness Index (AEI) (Villanueva-Rivera, Pijanowski et al. 2011) and Acoustic Complexity Index (ACI) (Pieretti et al. 2011). The BI and H indices are related to the amplitude and evenness of the biophonic sound frequencies, whereas the ADI and AEI are related to evenness among the pre-established frequency bandwidths. The ACI is based on the difference in the biophonic amplitude between the time intervals. For a more detailed description of the use of AIs, it is suggested to refer Bradfer-Lawrence et al. (2019). AIs calculations were conducted with soundecology version 1.3.3 (Villanueva-Rivera and Pijanowski 2016) and seewave version 2.1.8 (Sueur et al. 2008b). Before calculating the indices, we cut the frequency spectrum of the sound files according to the frequency range of the bat echolocation calls from 8 to 125 kHz. We calculated the indices with distinct bandwidth values for the AIs that use predefined frequency (ADI and AEI) or time (ACI) intervals. ADI and AEI were calculated with 8-, 10-, 15, and 20 kHz frequencies bandwidth sizes, and ACI was calculated using 1, 2, 5, and 10 seconds time intervals. We summarised the indices by calculating the mean value per index per sampling point across the two survey months (Supplementary Information, Table S3). Due to short time intervals during sampling campaigns, we assumed no influence of seasonality in our data.

Diversity metrics
We calculated one diversity metric based on Taxonomic Diversity (TD), one metric for Phylogenetic Diversity (PD), and three metrics based on Functional Diversity (FD). The first FD metric is based on Functional Diversity -morphology traits (FDm). The second FD is based on acoustic traits (FDa), while the third is based on a combination of morphological and acoustic traits (FDam), and the last one is based on Phylogenetic Diversity (PD). In calculating all diversity metrics, we considered a bat species only if the discriminant analysis successfully separated it. The species richness of the sample sites was defined as the number of distinct bat species identified at each sampling point. Functional morphological diversity was based on the dissimilarities between the cranial and wing measures of the species found in the sample sites. These morphological traits reflect the use of different food resources and foraging niches of aerial insectivorous bats (Bogdanowicz et al. 1999), while functional diversity indices based on the dissimilarity of morphological traits indicated a variation in the ecological function of a species in a community (Reiss et al. 2009). Morphological functional diversity was based on the dissimilarity of the species' cranial and wing measurements obtained from the literature (Supplementary Information, Table S4).
Cranial measures, diet, and echolocation calls are functionally related (Arbour et al. 2019), and wing traits, such as aspect ratio and wing load, are associated with foraging modes (Amador et al. 2020). Functional diversity based on acoustic trait dissimilarities can also reflect the foraging diversity of insectivorous bats, as traits such as frequency bandwidth, duration, and pulse interval are related to foraging mode and environment use (Jones and Holderied 2007;Luo et al. 2019). The acoustic dissimilarities between species were based on the bioacoustic traits of the bat calls recorded in our samples ( Supplementary Information, Table S3).
Phylogenetic diversity was calculated using the phylogenetic distances available from Shi and Rabosky (2015) ( Supplementary Information, Figure S5). To standardise the diversity metrics, Rao's quadratic entropy (Botta-dukát 2005) was used as a diversity index of the functional and phylogenetic diversities, calculated with the functions dbFD and Road in the packages FD version 1.0-12 (Laliberte and Legendre 2010) and picante package version 1.8.2 (Kembel et al. 2010). All metrics are shown in the Supplemental Material (Table S4).

Statistical analyses
We used generalised linear models to test and compare the efficiency of individual AIs in predicting the bat diversity metrics. The general structure of all models is the same: glm (diversity metric: acoustic index). As species richness data are discrete, we used Poisson error distribution for species richness models and Gaussian error distribution for functional and phylogenetic models (Zuur et al. 2009). We used the small-sample correct Akaike Information Criterion (AICc) values to compare the model fit, where the same diversity metric responds to distinct AIs, distinct frequencies, and time scales used in the calculating indices. Models with ∆AICc ≤2 were considered similarly adjusted (Burnham and Anderson 2002). Finally, we performed a likelihood ratio test for a more parsimonious model of each diversity metric. The general linear models were performed with the GLM function, and models AICc were calculated using the AICc function of the package AICcmodavg version 2.3-1 (Mazerolle and Mazerolle 2017). The spatial autocorrelation in the models was tested by calculating the Moran's I autocorrelation index of the model residuals using the function 'Moran.I' of the package ape version 5.5 (Paradis and Schliep 2019). We used the package vegan version 2.5-7 (Oksanen et al. et al. 2020) to perform a Mantel correlation test between bat morphological and acoustic normalised trait matrices with phylogenetic distances. We also calculated the Pearson correlation between the AIs and between the diversity metrics. All analyses were performed using R software (R Core Team 2020). We considered p < .05 as significant differences in all tests.

Results
The FDA satisfactorily grouped 21 previously identified bat calls (DI ≥ .8; Table 1), which were considered valid species (or genus) records for subsequent analyses. There was no correlation between phylogenetic distances and the distances between species based on morphological traits (Mantel r = > −.04, p = .387). However, there were significant correlations between phylogenetic distance and acoustic traits (Mantel r = .20, p = .01), and between morphological and acoustic traits (Mantel r = .40, p = .021). Table 1. Model selection of the bat diversity metrics predicted by distinct acoustic indices calculated with distinct sound frequency (Acoustic Diversity Index -ADI and Acoustic Entropy Index -H) and time scale (Acoustic Complexity Index -ACI). The calculations of the H and Bioacoustic Index -BI do not use frequency and time band size. The ∆aicc is relative to the models' comparison within the same index calculated in distinct scales. for all models, we used General Linear Models with the structure: diversity metric ~ acoustic index, family = "Gaussian"; except for taxonomic diversity, we used the "Poisson" family. There was no high correlation (here considered as r ≥ .7) between the diversity metrics (rTD × PD = .68, rTD x FDm = .16, rTD x FDa = .29, rPD x FDm = .33, rPD x FDa = .36, rFDm x FDa = .40), which indicates a non-correlation of these metrics.
The ADI and AEI were highly correlated (r = > −.86) as well and the ACI was highly correlated with H (r = .87) and BI (r = .74). Based on the AICc values, the H index was the most parsimonious predictor of species richness and phylogenetic diversity (Table 1).
Although the H index also significantly influenced the functional morphological diversity, this diversity metric was better explained by the ACI (Figure 2). No index could significantly predict acoustic functional diversity (Table 2), although the models with the H index as a predicted variable were the most parsimonious, with lower AICc values. None of the best-supported models showed significant spatial autocorrelation  The frequency bandwidth did not influence the response of the diversity metrics to ADI, AEI, or the time interval to ACI. Generally, models that used the same indices calculated with different frequency bandwidths (ADI and AEI) or different time intervals (ACI) showed a ∆AICc ≤2 (Table 1).

Discussion
Our results indicate AIs as a valid approach for the rapid assessment of bats in regions with a low level of biodiversity inventories and where native vegetation is rapidly converted to anthropic areas, as is the case of the Brazilian cerrado. As hypothesised, the diversity metrics were positively related to at least one AI. The H index and ACI had the best predictive potential. Contrary to our expectations, the predictive potential of the indices did not increase with changes in the frequency or time windows. Therefore, AIs are potentially suitable proxies for species richness, phylogenetic diversity, and functional diversity of bats in the savannahs of central Brazil. We accurately discriminated the presence of 21 taxa of insectivorous bats, a group challenging to catch with mist nets and often underrepresented in studies using traditional sampling methods (Hintze et al. 2020;Jung et al. 2014). Previous studies have shown that primarily insectivorous species account for 59.3% of the 118 bat species in the cerrado (Aguiar et al. 2016). Thus, there is a potential application of bioacoustics for surveying under-sampled areas in the cerrado (Bernard et al. 2010;Aguiar et al. 2020).
The index with the best predictive performance for bat diversity metrics is acoustic entropy (H), which measures the distribution of energy across acoustic space (Sueur et al. 2008a). The higher the H index, the higher the species richness (TD) and phylogenetic diversity (PD), with no association with functional acoustic diversity (FDa). Higher values of the H index are expected to be associated with species-rich environments because they tend to produce sounds that occupy more of the soundscape (Sueur et al.008a). This may be the case for our sampling points. The points found in the smaller protected area (AGCV), which has more anthropogenic disturbance, showed H values lower than the more extensive and conserved areas (PNB and PNCV), irrespective of the type of environment (grasslands, dense cerrado, or gallery forest).
The best predictor of functional morphological diversity (FDm) was ACI. However, the differences between the models using the ACI and H index as predictors were small (<2), indicating that both models were equally parsimonious (Burnham and Anderson 2002). The ACI is an index reflecting vocal activity and might indicate the abundance of individuals more than species richness (Retamosa et al. 2018). The strong association of ACI with FDm suggests that sampling points with high ACI values had more distinct species in cranial and wing measures, but not necessarily more species. Cranial and wing morphology, size, and echolocation call characteristics are predictive of their diets and play a significant role in determining the partitioning of food resources (Findley and Black 1983;Aldridge and Rautenbach 1987;Ramírez-Fráncel et al. 2021). Insectivorous bats differ in prey size, exoskeleton hardness, flight shape, body size, vocalisation and echolocation, dentition, and skull shape. Taxa that resemble each other most strongly morphologically are also the most similar in dietary intake. Considering the FDm calculation using wing and skull measurements, the range of values could reflect the coexistence of species that exploit a diverse variety of insects. A typical community of bats in the cerrado is composed of insectivorous bats of different sizes, ranging from 5 to 100 g (Paglia et al. 2012). The relationship between FDa and H index was strongly nonlinear, although it was not significant. Therefore, we considered that none of the tested indices could predict local variation in bat acoustic traits. The association of H and other indices with species richness has already been described, particularly for birds (Machado et al. 2017;Mammides et al. 2017;Bradfer-Lawrence et al. 2020). However, the bestperforming indices and response patterns vary considerably (Fairbrass et al. 2017;Ferreira et al. 2018;Moreno-Gómez et al. 2019).
To the best of our knowledge, this is the first study to implement AIs with distinct dimensions of bat diversity in the neotropical region. Many studies have also tested the effectiveness of AIs in representing species richness in several regions and found a significant relationship between some AIs (e.g., ACI, H, and BI) and predicted species richness. However, practically all studies conducted in Neotropics used recorders which record in the audible frequency range (e.g., Machado et al. 2017;Jorge et al. 2018;Ferreira et al. 2018;Retamosa Izaguirre et al. 2018;Campos et al. 2021;Duarte et al. 2021;Oliveira et al. 2021). As we showed, the study of bats requires ultrasound devices to properly record them, especially due to their wide echolocation range of 19.7 to 120.9 kHz Frequency of Maximum Energy (in the case of our study).
Considering our results, we highlight that applying AIs as an indicator of bat diversity is viable but must be carefully considered. The relationship between AI and species richness can be influenced by vegetation structure, abiotic conditions and background noise (Sueur et al. 2014;Buxton et al. 2018), which may lead to a significant variation in the correlation between the AIs and observed diversity (Mammides et al. 2017;Moreno-Gómez et al. 2019). This variation is reflected in the low precision of the models with a marked dispersion of H and ACI at sites with high diversity. Considering our sampling points were located inside protected areas, and data collection was conducted in 'normal' weather conditions (no rain or wind), we assumed that our analysis was not affected by these aforementioned factors. Therefore, we recommend that index comparisons control the environment and seasonal variation, the range of time and frequency of the faunal interest group, and the sampling effort, as signal masking and total record time affect the precision of the indices (Bradfer-Lawrence et al. 2019; Metcalf et al. 2021). Furthermore, the survey quality can influence the relationship between the AIs and biodiversity measures (taxonomic, phylogenetic, morphological, or acoustic functional diversities) since a small sampling effort can misrepresent the species. In our study, we registered near 30% of the 69 species of insectivorous bats and close to 40% of the genera expected for the cerrado (Paglia et al. 2012). Even considering that our recordings did not capture insectivores bat species from the Phyllostomidae family, we believe that the relationship between the AIs and the phylogenetic or morphological diversity would not change because we captured most of the lineages and a wide variety of morphological traits.
The size of the frequency bands for ADI and AEI and the size of the time band for ACI did not influence the model performance. Although the change in the band size influenced the parsimony in some model sets (∆AIC >2), there was no general tendency of model improvements with a systematic increase or decrease in the band size. The majority of AI studies use the default settings of the R package functions (Sueur et al. 2008b; Villanueva-Rivera and Pijanowski 2016) of 1000 Hz as the size of frequency bands (e.g., Fairbrass et al. 2017;Machado et al. 2017;Ferreira et al. 2018), which we consider narrow for bat species, considering several bat species have a large frequency bandwidth (Jones and Holderied 2007;Arias-Aguilar et al. 2018). Most species recorded at our site had a frequency bandwidth greater than 1 kHz and exceeded 30 kHz (Table 2). A narrow frequency band could overestimate local acoustic diversity, assigning distinct spectral signatures to the same species. Unlike birds or frogs, bats can emit different calls or vary the calls according to the type of environment in which they navigate, or even when they are alone or flying with other individuals (Habersetzer 1981;Obrist 1995). Thus, if the band frequency is too narrow, the AIs will be inflated, because the same species will appear in different bands. However, if the frequency band is too wide, several species would be considered at the same band frequency, underestimating the acoustic diversity. Although our data do not show a direct tendency in this direction, we recommend using the mean bandwidth frequency as the band frequency window of the species that possibly occurs in the study area. The mean bandwidth of the recorded species at our study sites was approximately 8 kHz (8.4), which proved to be a good frequency band for our models.
Then, the spectral structure of echolocation calls was phylogenetically structured. The phylogenetically closer species showed more similarity in their acoustic traits. Acoustic distances also correlated with species distances concerning morphological structures, reinforcing that echolocation traits are phylogenetically conserved and closely related to cranial morphology (Arbour et al. 2019;Luo et al. 2019;Roemer et al. 2019). Therefore, the sites with high phylogenetic and morphological functional diversity are expected to present high acoustic complexity and high values of acoustic entropy because both indices increase with biophony heterogeneity (Bradfer-Lawrence et al. 2019). However, acoustic functional diversity showed great variation at sites with a high H-index. Among the acoustic traits, the interval between pulses, which had a greater coefficient of variation and consequently had great weight in acoustic functional diversity, showed little influence on the H index, which was more influenced by pulse frequency variation (Sueur et al. 2008a). This discrepancy in the weights of the variables that influence both indices is responsible for their nonlinear relationship.
The fundamentals of AIs as indicators of biodiversity are the subdivision of acoustic space by distinct species (Krause 1993;Sueur et al. 2014). The partitioning of acoustic space by bat species with their echolocation calls in frequency bands and temporal patterns and the phylogenetic signal trends of these traits lead the indices H and ACI to predict, with decreasing precision, the taxonomic, phylogenetic and functional bat diversity. We conclude that certain indices, especially H, can be used to predict bat diversity, but should be used cautiously due to their relatively low precision and non-linearity. Furthermore, we call attention to properly set the necessary parameters to calculate the indices when using R packages. R packages, such as seewave (Sueur et al. 2008b), or soundecology (Villanueva-Rivera and Pijanowski 2016), have their default parameters defined for audible sounds, and the use on bats will require some specific tunings. Thus, we advise conducting preliminary analysis regarding the characterization of bat species calls present in a study area to define the frequency bandwidth or time window sizes adequately.