Functional diversity of Himalayan bat communities declines at high elevation without the loss of phylogenetic diversity
Study area and sampling locationsWe conducted this study in Kedarnath Wildlife Sanctuary (30° 25′–30° 41′ N, 78° 55′–79° 22′ E), located in Uttarakhand state in the western Himalayas of India. This sanctuary covers a broad elevational gradient from 1400 to 4000 m above sea level (asl) (Fig. 1), with corresponding changes in habitat types: from Himalayan moist temperate forests dominated by Quercus spp. at low elevations, to sub-alpine forests dominated by Rhododendron spp. and alpine meadows at high elevations34. This sanctuary is known to harbour 26 species of bats35.Figure 1Map of India showing the location of the study area, Kedarnath Wildlife Sanctuary, and the sampling locations within the study area. Elevation is in m asl. The map was created using QGIS (v 3.6.3-Noosa) (QGIS Geographical Information System, www.qgis.org). Please note that the geographical boundaries represented in the map may contain areas considered disputed.Full size imageWe sampled at four locations spanning an elevational gradient of 2200 m. Sampling points within each location were spread across the elevations mentioned in parentheses: Mandal (1500–1800 m), Ansuya (2000–2200 m), Chopta (2700–3000 m) and Tungnath (3300–3700 m) (Fig. 1). Sampling was conducted between late-March and mid-May in 2018 and 2019, starting at lower elevations and then moving to higher elevations. This sampling duration coincides with summer in the Himalaya. To comprehensively sample the bat diversity, we employed a combination of automated ultrasonic recorders and capture sampling using mist-netting. Fieldwork was approved by the Internal Committee for Ethics and Animal Welfare, Institute for Zoo and Wildlife Research (approval no. 2018-06-01), and conducted under a permit issued by the Uttarakhand State Forest Department, Government of India (permit no. 2261/5-6).Sampling strategyFor acoustic sampling, we placed full spectrum passive ultrasonic recorders (SongMeter SM4BAT, Wildlife Acoustics, Maynard, MA, USA) in different habitat types (open, forest edge, and forest) at each elevation (hereafter, “passive recordings”). The recorders were programmed to record bat calls for two consecutive nights at each sampling point, from dusk to dawn (9–10 h/night), using a sample rate of 500 kHz/s, an amplitude threshold of 16 dB and a frequency threshold of 5 kHz. The dominant habitats at Ansuya and Tungnath are montane evergreen forests and alpine meadows respectively, therefore only these habitats were sampled at these elevations. The exact number of sampling points per habitat for each elevation is given in Table S1. On separate days after completing acoustic sampling at a site, we set up nylon and monofilament mist nets of 4, 6 or 9 m length, 16 × 16 and 19 × 19 mesh sizes (Ecotone GOC, Sopot, Poland) for four hours following dusk (starting between 18.30 h in early summer and 19.30 h in late summer). The captured bats were handled and measured following the guidelines of the American Society of Mammalogists36. To further refine identification in light of the paucity of taxonomic knowledge in the region, we collected only one specimen each of taxonomically-challenging species in accordance with our field research permit. We measured body mass (accuracy 0.1 g) using a spring balance (Pesola, Schindellegi, Switzerland), and forearm length (accuracy 0.01 mm) with vernier calipers (Swiss Precision Instruments SPI Inc., Melville, NY, USA). Next, we gently stretched the left wing and placed the live animal perpendicular to the background of a graph sheet of 1 × 1 cm grids. We photographed the outstretched wing using a Nikon D3400 DSLR camera at 55 mm zoom from a distance of about 90 cm. Subsequently, we released the bats and recorded their echolocation calls at a distance of 5 to 10 m using a handheld ultrasonic detector (Anabat Walkabout, Titley Scientific, Brendale, QLD, Australia) and saved them as audio files of .wav format. These recordings (henceforth referred to as “reference recordings”) formed the dataset used to develop a call library for identification.Call classifier and analysis of passive recordingsReference recordings from 2018 and 2019 were labelled using Raven Pro 1.5 (Cornell Lab of Ornithology, Ithaca, NY, USA) to generate a dataset of acoustic parameters for identification. We visualized calls using a spectrogram with Hanning window, size 1024 samples with 95% overlap. From each recording, we selected 10 clear pulses and measured the following parameters: average peak frequency, maximum peak frequency, centre frequency, minimum peak frequency, peak frequency at the start and end of the call, bandwidth at 90% peak amplitude, average entropy, and call duration. All frequency variables were measured in Hz and time variables in ms. We used the peak frequency contour to determine start and end frequencies and also used bandwidth at 90% peak amplitude because higher frequencies attenuate quickly with distance from the emitting bat (causing changes to the bandwidth), and these measures are therefore more reliable in field circumstances. Using this labelled call library as a training dataset, we trained a fine K-nearest neighbours classifier using supervised learning within the ‘Classification Learner’ app in MATLAB (Mathworks, Inc., Natick, MA, USA). We further employed fivefold cross-validation to obtain estimates of the accuracy of each classifier in assigning calls to species. Using these pairwise values of relative accuracy (%), we generated confusion matrices for these classifiers where the species identities were represented in the columns and rows as ‘True’ and ‘Predicted’ classes, respectively. Any species with classification accuracy below 85% was clubbed with possible confusion species into a “sonotype”, to improve accuracy of the classifier in the most conservative way possible (Fig. S4). The complete list of sonotypes and their mean echolocation call parameters is presented in Table 1. The classifier identified these sonotypes with > 80% accuracy, with the exception of Miniopterus and the Plecotus type B call (which, however, we could manually identify because of their call structures and frequencies). For all subsequent analyses on functional diversity and phylogenetic diversity, we used these sonotypes to ensure accurate identification.Table 1 Trait matrix of the sonotypes in our assemblage (FA in mm; fmaxe, pfc.min, and pfc.max in kHz; Duration in ms).Full size tableNext, we analysed the passive recordings manually in Raven Pro. We labelled calls in subsets of 15 min per hour of the passive recordings. For each hour, the 15-min subsets were in the time windows 0–5 min, 20–25 min and 40–45 min, so as to spread out our sampling window across the hour. Following labelling, we obtained sonotype IDs using the classifier, and then verified them manually by visual comparison to the call library to improve discrimination. For every 5-min interval, we made a presence-absence matrix where 1 indicated the presence of a sonotype and 0 indicated its absence. The number of 5-min intervals in which a sonotype was detected (hereby “acoustic detections”) was summed up for each sampling point. We measured the relative abundance of sonotypes as the proportion of its total number of acoustic detections relative to the total number of acoustic detections of all sonotypes in a given elevational location. The use of such a presence-absence framework is akin to ‘Acoustic Activity Index’37 which represents a relatively less biased index of activity that is less affected by differences in vocal behaviour and echolocation frequencies of different species of bats.Assessing detectabilityTo assess the completeness of our species inventory, we estimated the species richness of each sampling point using the first-order Jackknife Estimator (Jack 1)38. Jack 1 is a nonparametric procedure for estimating species richness using presence or absence of a species in a given plot rather than its abundance39. Mean species detectability was calculated as the ratio of the observed to estimated species richness for different sampling point-year combinations40,41. We then assessed whether this mean species detectability depended on the habitat type, year, and location by fitting a linear model with the above-mentioned variables as fixed factor predictors and the mean detectability as a response. We also determined species-level detectability by following the approach of Kéry and Plattner42. If a sonotype was detected by mistnetting or acoustic sampling in sampling event i, we modelled its probability to be detected in sampling event i + 1. For each sonotype, we fitted a generalized linear mixed-effects model (logit link and binomial error distribution) with detection/non-detection as the response variable, and habitat type, location, and year as the fixed factor predictors. Site and species were included as random intercepts. The significance of the fixed effects was assessed with the Likelihood Ratio Test. This test allows one to choose the best of two nested models by assessing the ratio of their likelihoods. The significance of the random effect (species) was assessed by applying a parametric bootstrap (number simulations = 100) to the model with and without the random effect, using the function bootMer of ‘lme4’ package. In short, a parametric bootstrap consists of fitting the model to the data and bootstrapping the obtained residuals. For these and other statistical analyses we used R version 4.0.2 (R Core Team 2020).Taxonomic diversityWe calculated rarefied incidence-based species richness (SR) and Simpson diversity extrapolated to 50 sampling events (the number of sampling events in Mandal) using the ‘iNEXT’ R package43. The calculations were performed on a sonotype-by-sampling point presence-absence matrix with detections from both acoustic sampling and mistnetting pooled together. In the matrix, columns represented sampling units (Night 1, Night 2 and so on) and rows represented sonotype. By using sonotypes instead of species, we likely underestimated the SR, but this underestimation was uniform across elevations and is unlikely to change the pattern of SR with elevation.Functional diversityOur functional trait matrix (Table 1) comprised seven morphological and acoustic traits involved in guild classification, foraging and micro-habitat preferences (abbreviation followed by units): forearm length (FA, mm), aspect ratio (AR), wing loading (WL, N/m2), tip-shape index (I), echolocation peak frequency/frequency of maximum energy (FmaxE, kHz), minimum and maximum frequencies of the peak frequency contour (pfc.min and pfc.max, kHz) and call duration (D, ms). FA was measured in the field using vernier calipers. We used ImageJ (National Institutes of Health, Bethesda, MD, USA)44 to measure total wing area, areas of hand and arm wings and the wingspan from the standardised wing photos that were taken in the field. We calculated AR, WL, and I from these measurements following the equations given in Norberg and Rayner45. AR and WL both represent parameters that are correlated with flight aerodynamics and behaviour. I is influenced by the shape of the wing tip where values of 1 and above indicate broad, triangular tips, while those below 1 indicate acute wing tips. The four acoustic traits represent the shape of the echolocation call and they were measured from the reference recordings using Raven Pro, as described above.We first calculated the means for each of the seven traits across all species within a sonotype (thus obtaining one average trait value for each sonotype) (Table 1) and then used those to compute four multivariate functional diversity (FD) indices: functional richness (FRic), divergence (FDiv), evenness (FEve)46, and dispersion (FDis)47, using the function dbFD() in the ‘FD’ R package47. Our FD measures are unlikely to be underestimated due to the pooling of species into sonotypes because these species were similar in acoustic and morphological traits. FRic is the convex hull volume of the traits of species present in a community, measured in the multidimensional trait space. This measurement is not weighted by abundance, relative abundance or biomass of the species in the community, but, it is standardised such that it ranges from 0 to 1. FDiv reflects the distribution of abundance across taxa (sonotypes in our case) in the functional space. High FDiv means the taxa with extreme trait values are more abundant in a community whereas low FDiv means that those with the trait values close to the centre of the functional space are more abundant48. FEve, on the other hand, measures the evenness in the abundance distribution of taxa in the functional space. FEve is high when all taxa have similar abundances, and it is low when some functional groups are abundant while others are rare48. Lastly, FDis is measured as the mean distance of all taxa to the abundance-weighted trait community centroid. We performed two sets of analyses: one using the number of mistnet captures as a proxy for relative abundance, and another using the number of detections of different sonotypes in 5-min intervals in the passive recordings as a proxy of relative abundance. We did not pool acoustic detections and mistnet captures as they have inherently different detection probabilities and measure different entities (relative number of detections vs. number of captured individuals). Owing to rhinolophid bats at lower elevations being taxonomically and functionally different from the remaining species pool, we performed another set of FD calculations, excluding the four rhinolophid species and using acoustic detections as relative abundance. One species, Tadarida teniotis was commonly detected at all elevations on acoustic recorders, but we were unable to capture it as it foraged high above the ground, and thus were unable to collect morphological trait data. Additionally, in using acoustic detections as a measure of relative abundance, we had to exclude the non-echolocating pteropodid bat Sphaerias blanfordi which was caught only once at Chopta. Therefore, our FD values are likely systematically underestimated across all elevational communities, which does not affect the comparison of community composition across elevations.Phylogenetic diversityUsing the nexus file of a published phylogeny49, we pruned the tree to represent species in the 14 sonotypes. For each of these types, we chose the species most commonly mist-netted as representative of its group. Published DNA sequences are lacking for some of the species in this region, so we chose their closest relatives from the phylogeny instead. Thus, we made the following replacements: (a) Nyctalus leisleri represented the AMN sonotype, (b) Eptesicus serotinus represented the EH sonotype, (c) Murina aurata for Murina sonotype, (d) Myotis longipes for MS sonotype, (e) Pipistrellus javanicus for MP sonotype, and (f) Plecotus turkmenicus for Plecotus sonotype. After pruning the tree, we calculated three indices of phylogenetic diversity using the ‘picante’ R package50: Faith’s phylogenetic diversity (PD), Mean pairwise distance (MPD) and Mean nearest-taxon distance (MNTD). Faith’s PD is a measure of phylogenetic richness which is obtained by summing the branch lengths of the tree connecting the species in the community. MPD and MNTD measure phylogenetic dispersion of communities; whereas MPD measures the average phylogenetic distance among all the taxa in a community, MNTD measures the same for the nearest neighbouring taxa51. We weighted MPD and MNTD by relative abundance of the sonotypes in each community (like FD, the number of detections in five-minute intervals in the passive recordings was used as a proxy of relative abundance).Null model testingAs FD and PD are strongly correlated to species richness52, we used a null model to assess whether the observed was significantly different than expected due to chance alone. We produced the null distribution of each FD and PD index by randomizing the community matrix 999 times using the ‘independent swap’ method53,54, so as to preserve the species richness at each site and the number of sites in which each species can be found. Our randomization was further constrained by elevation, so that the abundances were randomized among the sampling points within each elevation. The null model allows for calculation of an effect size (difference between the observed value and mean of the null distribution). Given the range of FD and PD values, the effect sizes are not comparable across communities with vastly different species richness55. Therefore, standardized effect sizes (SES) of each index were calculated at each site as the difference between the observed value and the mean of the null distribution, divided by the standard deviation of the null distribution. SES > 1 and SES More
