Bat echolocation plasticity in allopatry: a call for caution in acoustic identification of Pipistrellus sp.

Animals modify their behaviours and interactions in response to changing environments. In bats, environmental adaptations are reflected in echolocation signalling that is used for navigation, foraging and communication. However, the extent and drivers of echolocation plasticity are not fully understood, hindering our identification of bat species with ultrasonic detectors, particularly for cryptic species with similar echolocation calls. We used a combination of DNA barcoding, intensive trapping, roost and emergence surveys and acoustic recording to study a widespread European cryptic species complex (Pipistrellus pipistrellus and Pipistrellus pygmaeus) to investigate whether sibling bat species could exhibit extreme echolocation plasticity in response to certain environmental conditions or behaviours. We found that P. pygmaeus occupied the acoustic niche of their absent congeneric species, producing calls with P. pipistrellus’ characteristic structure and peak frequencies and resulting in false positive acoustic records of that species. Echolocation frequency was significantly affected by the density of bats and by maternity rearing stage, with lower frequency calls emitted when there was a high density of flying bats, and by mothers while juveniles were non-volant. During roost emergence, 29% of calls had peak frequencies typical of P. pipistrellus, with calls as low as 44 kHz, lower than ever documented. We show that automatic and manual call classifiers fail to account for echolocation plasticity, misidentifying P. pygmaeus as P. pipistrellus. Our study raises a vital limitation of using only acoustic sampling in areas with high densities of a single species of a cryptic species pair, with important implications for bat monitoring. Ultrasonic acoustic detectors are widely used in bat research to establish species inventories and monitor species activity through identification of echolocation calls, enabling new methods to study and understand this elusive understudied group of nocturnal mammals. However, echolocation call signalling in bats is intrinsically different to that of other taxa, serving a main function of navigation and foraging. This study demonstrates an extreme level of plasticity, showing large variation in call frequency and structure in different situations. We showcase the difficulty and limitation in using acoustic sampling alone for bat monitoring and the complications of setting parameters for species identification for manual and automatic call classifiers. Our observations of call frequency variation correlated with density and absence of congenerics provide novel insights of behavioural echolocation plasticity in bats.


Introduction
In the face of huge biodiversity losses, accurate cataloguing of species and understanding their ecology is key to setting effective conservation priorities (Burgin et al. 2018;Frick et al. 2020). To identify species, we frequently rely on morphology and behavioural data that encompass complex interactions and responses of animals to their changing environments. Accurate species identification of bats can be particularly challenging due to their behaviour as nocturnal flying mammals, and their limited study compared to other taxa. Many bat species have gone undocumented over time (Jones 1997;Kunz and Parsons 2009;Altringham 2011), especially cryptic species-defined as two or more closely related species that could be easily misidentified due to their morphological similarity (Bickford et al. 2007). Where these species are not distinguished, erroneous records of intraspecies variation, distribution ranges, and niche width are recorded (Jones 1997). Therefore, for cryptic species, it is important to understand their unique morphology, behaviour and ecology and the limitations to differentiating between them. Cryptic species are widespread in bats (Kiefer et al. 2002;Jones and Barlow 2004;Srinivasulu et al. 2019) and have been increasingly revealed through DNA metabarcoding and acoustic technologies (Mayer et al. 2007;Simmons and Cirranello 2020).
Ultrasound detectors are widely used to record bat echolocation signals and provide a tool to investigate cryptic diversity where morphological differences are not easily recognised. Cryptic bat species have been discovered using acoustics by uncovering different phonic types within what was thought to be a single species. This happened with Scotophilus dinganii in Southern Africa (Jacobs et al. 2006), Hipposideros bicolor in Southeast Asia (Kingston et al. 2001) and Pteronotus cf. parnellii in South America (De Thoisy et al. 2014), amongst others examples (Mayer and von H e l v e r s e n 2 0 0 1 ; R a m a s i n d r a z a n a e t a l . 2 0 1 1 ; Chattopadhyay et al. 2012;Hintze et al. 2016). In Europe, one of the most researched and widespread bats, Pipistrellus pipistrellus, was divided taxonomically following acoustic, behavioural and genetic studies (Barlow 1997;Barlow et al. 1997; Barratt et al. 1997;Vaughan et al. 1997a). Primarily, the split was driven by the discovery of two phonic types with echolocation peak frequencies of~45 kHz and~55 kHz (Jones and Van Parijs 1993), later attributed to P. pipistrellus and P. pygmaeus respectively (Jones and Barratt 1999). This frequency partitioning allows occupation of separate acoustic niches, enabling coexistence (Jones and Holderied 2007).
However, echolocation calls are highly variable, and we continue to uncover new environmental and ecological factors that shape bat echolocation behaviour. Bat acoustic signals have been shaped through natural selection to favour optimal echoes from their prey  and to facilitate intraspecific communication (Jones 1997). Call structure and frequency is affected by habitat type and foraging strategies (Kalko and Schnitzler 1993;Kalko 1995). Bats flying in cluttered habitats tend to emit short broadband calls at high repetition rates (e.g. Myotis spp.) or constant frequency calls (e.g. Rhinolophus spp.), and bats in open habitats use quasiconstant frequency long narrow-band calls (e.g. Nyctalus spp.) (Kalko 1995). Echolocation adaptations to habitat have evolved convergently, where even across genera, different species use similar call structures when foraging in similar habitats (e.g. Siemers and Kalko 2001). Furthermore, bats modify their echolocation frequency in the presence of other bats to avoid overlap with each other's calls, which is known as echolocation jam avoidance behaviour (Ulanovsky et al. 2004;Gillam et al. 2007).
Ecological plasticity-the ability of an organism to alter its behaviour, physiology or morphology in response to an environment (Price et al. 2003)-is broadly documented in bats in the form of echolocation frequency modulation (Russo et al. 2018). Bats exhibit plasticity in echolocation signals across different seasons and geographic areas (Barclay and Brigham 2004;Snell-Rood 2012;Mutumi et al. 2016;López-Baucells et al. 2018), in response to temperature and humidity (Guillen et al. 2000;Chaverri and Quirós 2017), and in relation to sex, age and anatomy (Siemers and Kalko 2001). In fact, Lameira et al. (2010) conducted a review of acoustic signals in terrestrial mammals and found that bats were the taxon with the most studies showing geographic variation. Understanding this capacity of bats to modify their echolocation calls is crucial for bat studies, which increasingly rely on acoustic methods.
Technological advancements in bioacoustics have improved our capacity to record and store large volumes of high-quality acoustic data. Accessibility to bioacoustics is also increasing with the release of low-cost and compact full spectrum detectors (e.g. AudioMoth from Open Acoustic Devices; Hill et al. 2019). In parallel, new software developments now enable the use of artificial intelligence for automated species identification of bat echolocation data (e.g. Kaleidoscope by Wildlife Acoustics Inc., USA; Tadarida by Bas et al. (2017), France; and Bat Classify by Chris Scott and John Altringham, UK), including regionally specific automated classifiers (e.g. Bertran et al. 2019). However, automated classifiers still have significant error rates in identification, often heightened in areas with sympatric species with similar echolocation calls (Russo and Voigt 2016;Rydell et al. 2017). The limitations of classifier algorithms in accounting for acoustic signal plasticity are not fully understood-and may be particularly flawed in distinguishing between the echolocation calls of cryptic species.
Given that bats exhibit highly plastic echolocation behaviour, oversimplified and overconfident approaches to acoustic identification could be resulting in inaccurate species monitoring. We aimed to examine this issue focusing on a pair of cryptic species: P. pipistrellus and P. pygmaeus. The echolocation calls of these sibling species are generally classified following parameters aligned with Russ (2012), with peak frequencies of 41.6-50.6 kHz for P. pipistrellus, and a higher range of 50.2-64.1 kHz for P. pygmaeus. While this classification is accurate for the species in some areas of their distribution, there are very few studies that investigate echolocation plasticity to determine how call frequencies vary in areas where these species do not occur together or how echolocation is affected by population density. Previous research has shown intercolony variability in echolocation calls of P. pipistrellus (Fornuskova et al. 2014) and plasticity in echolocation of P. pygmaeus while foraging in different habitats (Bartonička and Řehák 2005;Bartonička et al. 2007). Moreover, studies in a similar species, Pipistrellus kuhlii, showed echolocation variability in different behavioural contexts (Berger-Tal et al. 2008). We aimed to build on previous studies of echolocation plasticity by investigating the extent of variability and its causes, focusing on P. pygmaeus in a study area where these bats occur at very high densities and have been monitored since 2000.
In the Ebro Delta in Catalonia, bat boxes were installed in the 1990s with huge success in occupancy, with a total of over 3500 P. pygmaeus using the boxes by 2006 (Flaquer et al. 2006). To date, two PhD, numerous masters projects and national monitoring programs have been focused on bats in the delta, encompassing box checks, captures, acoustics, radiotracking and pest-control ecosystem service studies (Flaquer et al. 2005(Flaquer et al. , 2006(Flaquer et al. , 2009(Flaquer et al. , 2014Flaquer and Puig-Montserrat 2012;López-Baucells et al. 2013;Puig-Montserrat et al. 2015;Bideguren et al. 2019). Acoustic monitoring consistently shows abundant records with peak frequencies corresponding to P. pipistrellus, as well as P. pygmaeus, according to the parameters in Russ (2012) and Barataud (2015). However, in 20 years of extensive monitoring, no individuals of P. pipistrellus have ever been identified, from roosts or trapping surveys, using morphology and DNA barcoding (CF, pers. observ.). This has raised important scientific questions regarding whether the lower frequency calls are in fact P. pipistrellus or if they are being misidentified due to behavioural echolocation plasticity of its sibling species and, if so, which factors underlay this behaviour.
This study investigates bat echolocation plasticity and to what extent it can compromise acoustic species identification of cryptic bat species. We test the hypothesis that P. pygmaeus exhibits extreme and previously undocumented echolocation plasticity and overlaps with the call parameters of P. pipistrellus. We predict that:

i)
Based on extensive data from the past 20 years, all sampled bats of this cryptic species pair will be confirmed morphologically and genetically as P. pygmaeus.

ii)
Bats confirmed as P. pygmaeus will be recorded echolocating at frequencies typical of P. pipistrellus, and the presence of a single species will be supported by a unimodal frequency distribution of echolocation calls recorded through passive acoustic monitoring. In contrast, a bimodal frequency distribution would suggest the presence of both cryptic species. iii) Given that the lower frequency echolocation calls (typical of P. pipistrellus) are expected to be emitted by P. pygmaeus which are abundant and widely distributed in the Ebro Delta, these calls are expected in all habitat types and throughout all hours of the night. In contrast, a patchy geographical and temporal distribution of these lower frequency calls would suggest they are emitted by low numbers of P. pipistrellus that have avoided detection during sampling. iv) The lower frequency echolocation calls will comprise a higher proportion of the echolocation calls in areas of high bat activity compared to areas of low bat activity, as bats modify their echolocation frequency in the presence of other bats to avoid overlap with each other's calls.
Finally, we also investigate the echolocation of adult bats only compared to adult and juvenile bats together and the accuracy of bat species identification from acoustic recordings using an automatic classifier and manual classification via an online questionnaire.

Study area
The study was carried out in the Ebro Delta in the north-east of the Iberian Peninsula (0°43′ 3.69″ E 40°42′ 21.61″ N). This is one of the largest European deltas and consists of wetlands (65 km 2 ), crops (240 km 2 , of which~85% are rice paddies), urban areas (16 km 2 ) and 50 km of coastline with salt marsh areas (Puig-Montserrat et al. 2015). It has a semi-arid Mediterranean climate with a dry summer and a mean annual temperature of 16.8°C (Curcó 2006). The urban areas have high densities of houses and multi-storey buildings, and only few scattered buildings occur outside of the urban areas. There are no underground roosts, and trees are scarce, with scattered eucalyptus trees (Eucalyptus globulus), palm trees (Washingtonia filifera) and poplars (Populus alba). Bat boxes were first set up in 1999, starting with 69 bat boxes (Flaquer et al. 2006) and increasing to~500 in 2019 (X. Porres, pers. commun.). Over 15 bat box models are used (Bideguren et al. 2019), set up on electric poles, trees and buildings.

Trapping and box check surveys
Surveys were conducted between June and August 2019. Sampling sites ( Fig. 1) covered four habitat types: urban (n = 5), rice paddies (n = 5), lagoons (n = 5) and salt marshes (n = 5) across the Ebro Delta. Trapping was carried out twice at each site, a month apart, over 4 h starting at sunset, using traditional sampling methods-ground-level mist nets and harp traps (Kunz and Parsons 2009). Mist net numbers and sizes were selected based on habitat suitability and access limitations, with 6 to 14 mist nets and up to three harp traps set each night. Acoustic lures (one Apodemus BatLure and one Binary Acoustic Technology AT100 lure) were used to broadcast social calls of Pipistrellus sp. (Avisoft-Bioacoustics 2009; Barataud and Guido 2019) to attract the target species and increase sampling efficiency (e.g. Russ et al. 1998;Hill and Greenaway 2005;Lintott et al. 2013). Traps were checked every 15 min. Species identification was done following field guides Flaquer and Puig-Montserrat (2012) and Dietz and Kiefer (2016). A wing tissue sample (2-mm biopsy punch) was collected from adult bats and preserved in 96% ethanol for molecular species identification. Buccal saliva swabs of all bats were collected as part of a separate epidemiology project and were used for molecular species identification of juvenile bats. Over the two months, 390 bat boxes were also inspected across the Ebro Delta to check for presence of P. pipistrellus. As this study involved focal animals in the field, blind data recording was not possible.

Static passive acoustic monitoring
Full spectrum acoustic data was collected in 2016 and 2019 at 20 trapping sites ( Fig. 1) using Song Meter passive acoustic detectors (Wildlife Acoustics Inc., USA) and AudioMoths v.1.1.0 (Hill et al. 2019), as summarised in the supplementary material S1. The Song Meter detectors were programmed with a frequency range of 12-196 kHz, a trigger level of 12 dB and a high pass filter at 12 kHz and were set out for six nights in 2016 to record from 30 min before sunset to 30 min after sunrise at each site. AudioMoths were programmed for continuous recording at a sample rate of 256 kHz and medium gain and used to record for 4 h per night at each site during the 2019 trapping surveys.

Roost emergence surveys
Emergence surveys were conducted at 12 P. pygmaeus maternity roosts to investigate plasticity in echolocation calls as bats Fig. 1 Map of the Ebro Delta study site, showing the 20 trapping and acoustic sampling points distributed across four habitat types. Habitat data was downloaded from land cover maps from the Catalonian Department of Territory and Sustainability (http://territori.gencat.cat/ca/ inici/) and adapted using ArcGIS version 10.6 emerged. Two surveys were carried out at each roost, starting 15 min before sunset until~40 min after sunset. Pipistrelles are early emerging species (Jones and Rydell 1994;Duvergé et al. 2000), and during this time most bats left the roost. The end time was adjusted based on visual observations at each site, to avoid recording bats arriving from other locations. Surveyors stood~5 m away from the roost and recorded calls using an EMT2 Pro (Wildlife Acoustics Inc., USA). Faecal samples were collected from each roost for DNA species confirmation.

Data analysis
DNA high-throughput metabarcoding and bioinformatics DNA was extracted from 97 wing punches following a hot sodium hydroxide and tris (HotSHOT) genomic DNA preparation method (Truett et al. 2000). DNA was also extracted from 24 faecal samples collected from 12 bat roosts using a Stool DNA isolation kit (Norgen Biotek).
Sequences were obtained using a massive-parallel barcoding strategy (Ershova et al. 2019). The polymerase chain reaction (PCR) was conducted using the universal Leray-XT primer set (Wangensteen et al. 2018) which amplifies a 313 base pair (bp) fragment of the mitochondrial cytochrome c oxidase subunit I (COI) gene from a broad taxonomic range of eukaryotic species, including most metazoans. The primers used were forward primer mlCOIint-XT (5′-GGWACWRGWTGRACWITITAYCCYCC-3′) and r e v e r s e p r i m e r j g H C O 2 1 9 8 ( 5 ′ -TAIACYTCIGGRTGICCRAARAAYCA-3′). DNA amplification followed a 1-step PCR protocol, and the products were pooled into two libraries and prepared for sequencing in a MiSeq™ Reagent Nano Kit v2 (2 × 250 cycles).
Barcoding sequences from FASTQ files were analysed using a metabarcoding pipeline based on OBITools (Boyer et al. 2016). Taxonomic assignment was done using the OBITools ecotag algorithm against a reference database of COI Leray sequences, downloaded from the NCBI database and the BOLD database (Ratnasingham and Hebert 2007;NCBI 2020). See supplementary material S2 for further details of the metabarcoding and bioinformatics.

Classification of species in acoustic recordings
To classify bat calls from 2016, we followed a protocol developed by the Natural Sciences Museum of Granollers, as described in Tuneu-Corral et al. (2020). Bat recordings were split into 5-second files using Kaleidoscope v.4.3.2 (Wildlife Acoustics, USA) and were run through Tadarida, an automated acoustic classifier for European bat species (Bas et al. 2017). A manual post-validation of all files was completed using Avisoft-SASLab Pro 5.2.12 software (Avisoft Bioacoustics, Germany) to correct or validate the automatic classification. Each file was identified to species or conservative sonotypes (mixed-species groups) when unambiguous species determination was not possible. This method results in a standardised measure of bat activity as bat passes, defined as a minimum of two distinguishable echolocation calls per species within a 5-second sound file (e.g. Torrent et al. 2018). To classify Pipistrellus calls recorded during the trapping surveys in 2019, we automatically validated and used bat passes with a probability score of ≥ 0.9 (see supplementary material S3).
Additionally, to test the reliability of manual species classification by chiropterologists, we prepared an anonymous online questionnaire (see supplementary material S4).
Participants were asked to analyse 20 sound files from the Ebro Delta and identify the bat species. The specific purpose of the questionnaire (testing the reliability of classifying pipistrelle calls) was undisclosed, and other species were included to avoid biases in the answers.
Assessing echolocation plasticity in acoustic recordings of Pipistrellus sp.
For each confirmed Pipistrellus sp. bat pass, the following measurements of individual pulses were extracted with the Tadarida software: peak frequency, start frequency, end frequency, duration and modulation (or bandwidth). All extracted measurements were filtered to remove call pulses from other species and to retain only P. pipistrellus/P. pygmaeus call pulses using the conservative filters: peak frequency 42-65 kHz, end frequency 42-61 kHz, start frequency 45-110 kHz and duration 4.5-10 ms. These filters were selected based on Russ (2012) and Barataud (2015) as guidance, accounting for parameters suitable for the open habitats of the study area.
To examine if the peak frequency of echolocation calls was related to the density of bats flying, we assessed the density of bats in each bat pass (single bat vs multiple bats). To standardise classification, considering an inter-pulse interval (IPI) of~100 ms for P. pygmaeus (Russ 2012), we assumed that in one 5second file, a single bat could emit a maximum of 50 pulses. Files containing more pulses than this threshold were considered to have multiple bats. The IPI of 100 ms used accounts for an increase in IPI in uncluttered situations (Russ 2012). Subsequently, we verified this method by manually checking 100 files from each of 12 frequency categories, starting at a category of 43-44 kHz and incrementing by 1 kHz until 54-55 kHz. These categories cover Pipistrellus call peak frequencies recorded below the typical literature cut-off for P. pygmaeus calls, up to the reported mean peak frequency (~55 kHz).

Statistical analyses
To visualise the distribution of echolocation peak frequencies during emergence surveys at maternity roosts before and after the juveniles started flying, we used a kernel density estimation. We used data from six maternity roosts where the first emergence survey was conducted while juveniles were nonvolant (confirmed as newborn bats still unable to fly during visual inspections of the roosts), and the second survey was conducted when both mothers and juveniles were flying out. To examine the effect of the age and colony rearing stage (explanatory variable: before and after juveniles started flying) on the echolocation peak frequency (response variable), we used a linear mixed model. Roost was fitted as a random effect to account for roost variation in echolocation peak frequency. The model log(Peak frequency)~survey timing + (1| roost) was fitted using the "lme4" R package (Bates et al. 2015). A log transformation was used to normalise the frequency variable.
To examine the relationship between the echolocation pulse peak frequency (explanatory variable) and the modulation of the pulses (response variable), we performed a linear regression.
To assess the effect of habitat type, total bat activity and detector type (explanatory variables) on the number of echolocation pulses with peak frequencies of < 51 kHz (those below the typical literature cut-off for P. pygmaeus) per site per night (response variable), we used a generalised linear model. To deal with overdispersion of count data, we used a quasi-Poisson distribution. A post hoc Tukey's HSD test was used to examine differences between habitats using the "multcomp" R package (Hothorn et al. 2008). The "effects" R package (Fox 2003;Fox and Weisberg 2019) was used to visualise and plot the results of the models.
To estimate the proportion of files containing multiple versus single bats (binary variable), bootstrap sampling was used with the "boot" R package (Davison and H i n k l e y 1 9 9 7 ; C a n t y a n d R i p l e y 2 0 1 9 ) . T h e bootstrapping took random samples from 12 subsets of files that contained echolocation pulses with peak frequencies of 43 kHz to 55 kHz (each frequency integer value was a subset that contained calls of that particular frequency, with repetition between subsets where files contained call pulses of multiple peak frequencies) and calculated the proportion files in the subset that had multiple bats. This was done with 10,000 replicates for each subset. To model the effect of the peak frequency of the echolocation pulses in the files (explanatory variable; 12 categories) on the proportion of sound files containing multiple bats (response variable), a quasibinomial generalised linear model was used.
All data management and analyses were done using R version 3.6.1 (R Core Team 2019).

Trapping and box check surveys
Across 20 sites (Fig. 1), and after 160 h of trapping over 40 nights, a total of 233 bats of five species were captured, including 207 P. pygmaeus. All bats morphologically identified as P. pygmaeus were confirmed genetically via the DNA analysis of wing punch samples and buccal saliva swabs. Of the 390 boxes checked, 50.5% were occupied by P. pygmaeus. No P. pipistrellus were trapped at any site or recorded in any bat box or roost surveyed throughout the full summer survey period.

Roost emergence surveys
Emergence surveys at 12 P. pygmaeus roosts (species confirmed via DNA analysis of faecal samples) were acoustically recorded. Eight of the 12 roosts were bat boxes that were opened during the day, and the bats were also morphologically identified as P. pygmaeus. Pulses with peak frequencies below the typical lower end for P. pygmaeus (51 kHz) were recorded in all roosts and comprised 28.2% (n = 27,996) of the pulses. The lowest peak frequency recorded was 44 kHz, and the average was 54.9 kHz (SD = 4.8, n = 99,294). However, the mean peak frequency of emergence echolocation was significantly different before and after the juveniles started flying (GLMM, mean before = 51.0 kHz, mean after = 55.4 kHz, p < 0.001, n = 55868; full model results in supplementary material S5) with 3.6% of peak frequency variance explained by differences between roosts. Before the juveniles started flying (in other words only adults were recorded), 59.9% of the pulses had peak frequencies < 51 kHz, while after the juveniles started flying (both adults and juveniles recorded), only 17.1% of the pulses were below this value (Fig. 2).

Static passive acoustic monitoring
The 2016 acoustic data comprised a total of 257,031 bat passes, and 213,077 of them contained P. pygmaeus bat passes. For the 2019 acoustic data recorded during the trapping surveys, a total of 30,825 P. pygmaeus bat passes were confirmed. Individual echolocation pulses with peak frequencies of < 51 kHz comprised 36.3% (n = 2,391,007) of the extracted pulses from the bat passes (n = 6,572,334) from both years. The frequency bandwidth (or modulation) of the pulses was significantly correlated with the peak frequency of the pulse (lm, R 2 = 0.18, n = 6,572,334, p < 0.001), with lower peak frequency echolocation pulses generally less modulated, or flatter, than those at higher peak frequencies (see supplementary material S6).

Echolocation peak frequencies recorded during trapping surveys (2019)
Recordings of Pipistrellus echolocation calls collected alongside the trapping surveys over the 40 nights had an average peak frequency of 53.6 kHz and SD of 3.8 kHz (Fig. 3).

Effect of habitat type and night hour on echolocation peak frequency
We found that habitat had a significant effect on the number of echolocation pulses with peak frequencies < 51 kHz when modelled considering total bat activity and acoustic detector type (glm, chisq = 332,171, df = 3,162, p < 0.001, full model results in supplementary material S7). Echolocation pulses with peak frequencies below 51 kHz were recorded in all four habitats, and further post hoc tests showed that the number of pulses with peak frequencies < 51 kHz was significantly higher between lagoons and the other three habitats (rice paddies, salt marshes and urban areas), but no significant differences were found between the other habitats ( Fig. 4 and supplementary material S7).
In terms of acoustic activity across time of night, pulses with peak frequencies < 51 kHz were recorded throughout all hours of the night in similar proportions (Fig. 5). On average, during the night long recordings in 2016, 29.0% (range 24.0-33.2%) had peak frequencies below 51 kHz, but no hourly trends were detected.

Effect of flying bats density on echolocation peak frequency
Bat echolocation peak frequencies were affected by the density of flying bats (Fig. 6). We found a significant effect of the proportion of files with multiple/single bats (binary variable) on the peak frequencies of echolocation of the bats within the file (glm, p < 0.001, full model results in supplementary material S8). Lower peak frequency echolocation calls were found more commonly in files with multiple bats, whereas the higher peak frequency echolocation calls were more common in files with single bats. The analyses showed a clear pattern of decrease from 36% of files with pulses at peak frequencies of 43-44 kHz containing multiple bats to 19-20% of files with peak frequencies > 51 kHz with multiple bats (Fig. 6). Fig. 4 Predicted effect of habitat type on the amount of Pipistrellus echolocation pulses with peak frequencies < 51 kHz per site per night, balanced by the total number of pulses per site per night and the acoustic detector type used, modelled using a generalised linear model. Data collected in summer 2016 and 2019. Significant effects (p < 0.001) against all other habitats are indicated with an asterisk. Error bars indicate 95% confidence intervals Fig. 5 Pipistrellus bat echolocation pulses per hour split by peak frequency of echolocation for pulses below 51 kHz (those below the typical literature cut-off for P. pygmaeus) and grouped together for all pulses above 51 kHz. Data from summer 2016 across 20 sampling points in four different habitat types

Misidentification of automatic identification classifiers and manual classification
From the 2016 bat passes (n = 257,031), 31.3% (n = 80, 407) were assigned to P. pipistrellus by the Tadarida automatic identification classifier (false positives). From the 62 manual classification responses collected from the online questionnaire, all Pipistrellus sp. echolocation call files contained both P. pipistrellus and P. pygmaeus classifications. Therefore, under the hypothesis of echolocation plasticity of P. pygmaeus and no P. pipistrellus in the Ebro Delta, the answer P. pipistrellus resulted in false positives in 6.5% of answers for a call sequence with frequencies of 50-57 kHz, 67.8% of answers for a call sequence with frequencies 48-49 kHz and 90.3% of answers for a call sequence with frequencies 46-50 kHz (see supplementary material S9).

Discussion
This is the first study to document and examine the behavioural and ecological causes of extreme echolocation plasticity in P. pygmaeus, leading to a considerable overlap with the call frequencies of its cryptic sibling species P. pipistrellus. All evidence gathered supports the hypothesis that P. pygmaeus exhibit extreme echolocation plasticity in response to factors including the density of conspecifics flying in the area and the age of the bats. Our surveys did not confirm any P. pipistrellus, and the species identification of all captured and roosting P. pygmaeus bats was confirmed genetically, echoing the results of the last 20 years of monitoring.

Cryptic pipistrelles versus echolocation plasticity hypotheses
P. pipistrellus is readily trapped using traditional mist net and harp trap methods (Flaquer et al. 2007), with higher capture rates using acoustic lures (Lintott et al. 2013), and has been recorded in bat boxes (Park et al. 1996;Mcaney and Hanniffy 2015). Although we implemented a combination of these survey methods intensively during the whole summer period, no P. pipistrellus were captured or confirmed at any site, bat box or roost during this study. Nevertheless, as we predicted, echolocation pulses within the typical peak frequency range for P. pipistrellus were recorded across all sites and habitats and all hours of the night. Despite a third of the echolocation calls corresponding to call parameters of P. pipistrellus, no individuals of this species were found. Lintott et al. (2013) reported a significant positive relationship between bat activity recorded with acoustic detectors and the relative abundance of P. pipistrellus and P. pygmaeus from captures. Therefore, our acoustic records would suggest the presence of a large population of P. pipistrellus, spread evenly across all of the Ebro Delta, yet none have ever been found. Across all habitats and detectors, the first low frequency calls were recorded around sunset, indicative that bats come from nearby roosts, yet the availability of roosts outside of the urban areas (i.e. rice paddies, salt marshes and lagoons) is limited to bat boxes (Flaquer et al. 2006) that have been Fig. 6 Percentage of files with multiple Pipistrellus bats for files containing echolocation pulses of each frequency category. Boxplots were built by a bootstrap sampling analysis (replicates = 10,000) and show five summary statistics (the thick black lines represent the medians, the boxes encompass the interquartile ranges, and the whiskers extend to the data points within × 1.5 the interquartile range outside the box, and the circles show data points beyond the whiskers) monitored since 2000 and have never been occupied by P. pipistrellus during inspections. Moreover, as expected, the distribution of peak frequencies of calls recorded using passive acoustic detectors during trapping surveys (Fig. 3) follows a unimodal distribution ranging from 44 to 65 kHz, in contrast to the bimodal distribution that is typically noted in the presence of both cryptic pipistrelles (Jones and Van Parijs 1993;Park et al. 1996;Barratt et al. 1997;Russo and Jones 2000), or other cryptic species pairs (Thabah et al. 2006). P. pygmaeus were recorded echolocating as low as 44 kHz during emergence surveys, with 29% of the pulses within the typical frequency range for P. pipistrellus. This result, backed with genetics, proves that P. pygmaeus can lower its echolocation frequency and overlap with those of P. pipistrellus. Previous research of pipistrelles also found emergence call frequencies to be highly variable compared to other flight situations (i.e. foraging, commuting and returning to roosts) and attributed this to the dichotomy between calls before and after take-off (Berger-Tal et al. 2008). In Tadarida brasiliensis, a distinct echolocation call is used during emergence from roosts with high densities of individuals (Gillam et al. 2010), exemplifying echolocation plasticity during emergence.
Based on the results of this study, coupled with bat monitoring results from the last 20 years, we assume the likely absence of P. pipistrellus in the Ebro Delta and support the hypothesis that P. pygmaeus exhibits a higher degree of echolocation plasticity than ever reported.

Echolocation plasticity in P. pygmaeus in the Ebro Delta
According to our results, there are several factors that could drive this plasticity in echolocation calls: (a) habitat and environmental conditions, (b) absence of congeneric species and the availability of empty acoustic niches, (c) high flying bat densities and (d) age of the bats and, possibly, the energetic cost during pregnancy and lactating periods.
Narrowband, longer calls with lower peak frequencies are better suited for prey detection in open area habitats (Kalko 1995;Jones 1997;Schnitzler and Kalko 2001). It is well documented that higher frequency calls are more costly for bats as they are less energy efficient compared to low frequency calls that travel further and suffer less from attenuation (Lawrence and Simmons 1982). Furthermore, in areas of high humidity and temperature such as river delta wetlands or other aquatic habitats, attenuation of high frequency sounds is more pronounced, decreasing their detection distance compared to lower frequency calls (Kober and Schnitzler 1990). The habitats in the Ebro Delta, mainly open areas, markedly favour P. pygmaeus, which is associated with water habitats and riparian areas, compared to P. pipistrellus, which are more associated with woodlands (Vaughan et al. 1997a;Davidson-Watts et al. 2006). This combination of open habitats and certain environmental conditions might favour the use of lower frequency calls in P. pygmaeus. Although our study focuses primarily on the peak frequency parameter of echolocation calls because it is the most distinguishing variable between these two species (Russ 2012;Barataud 2015), we also found that the frequency bandwidth was narrower in calls characterised by a lower peak frequency. These flatter, lower frequency calls are typically emitted in open habitats (Kalko 1995;Schnitzler and Kalko 2001;Russ 2012), showcasing habitat type as an important factor in the frequency plasticity we see in P. pygmaeus echolocation calls.
In the absence of ecologically similar congeneric species, bats may selectively lower their peak frequencies to maximise their call efficiency. Buckley et al. (2011) suggest that this occurs with Nyctalus leisleri in Ireland in the absence of its congeneric species, Nyctalus noctula, that typically occupies the lower frequency niche. Studies comparing acoustic signals where species occur in sympatry and allopatry in various taxa report larger differences between signals in sympatric populations, for example, in African tinkerbirds (Kirschel et al. 2009), and green tree frogs (Höbel and Gerhardt 2003). In bats, Russo et al. (2007) demonstrated that in Sardinia, where three rhinolophid bat species occur in sympatry, the species echolocating at the highest frequency (Rhinolophus hipposideros) and at the lowest frequency (R. euryale) increased and decreased their frequency respectively to diverge their signals away from those of R. mehelyi which calls at an intermediate frequency between the two. In contrast, in southern Italy, in the absence of R. mehelyi, the gap in frequency between the other two rhinolophids is not as wide. This frequency partitioning enables occupation of separate acoustic niches (Jones and Holderied 2007) and facilitates intraspecific communication (Jones and Barlow 2004;Russo et al. 2007).
We suggest that in the absence of P. pipistrellus, P. pygmaeus can lower their echolocation peak frequency as a behavioural response as they do not need to diverge acoustic signals to facilitate communication or occupy separate foraging niches. The absence of its congeneric species may allow P. pygmaeus to occupy a broader acoustic niche, using calls with a broader range of peak frequencies, including those typically occupied by P. pipistrellus. Nevertheless, it is also possible that the frequency variation that we observe in P. pygmaeus echolocation also occurs where the two species occur in sympatry. In these situations, the extreme variance would be masked by the presence of heterospecifics, making it more difficult to detect. Further studies systematically comparing the patterns of echolocation frequencies in areas of sympatry versus allopatry of P. pipistrellus and P. pygmaeus are required to systematically test this. Comparative studies between sympatry and allopatry in bats are rare and are needed to improve our understanding of resource partitioning (e.g. of acoustic space) in relation to behaviour and echolocation traits (Salinas-Ramos et al. 2020).
Additionally, the peak frequencies of call pulses were significantly affected by the density of flying bats. As predicted, P. pygmaeus exhibited higher echolocation plasticity and echolocated at lower frequencies when flying in groups. Moreover, lagoon habitats had the highest bat densities and had significantly more echolocation pulses with lower (< 51 kHz) peak frequencies compared to the other three habitats. This suggests that the high abundance of P. pygmaeus bats in the Ebro Delta influences echolocation plasticity, causing conspecifics to modify their call frequency to avoid overlap with each other. P. pygmaeus, P. pipistrellus and P. nathusii have all been shown to vary their echolocation call peak frequencies in the presence of conspecifics, enhancing the differences in calls between individuals (Bartonička et al. 2007;Necknig and Zahn 2010), although not to the extent found in this study. Some studies argue that bats modify their echolocation frequency mainly as a response to clutter (Obrist 1995;Cvikel et al. 2015;Fawcett and Ratcliffe 2015;Götze et al. 2016), whereas others report a jamming avoidance response to the presence of other bats (Ulanovsky et al. 2004;Bartonička et al. 2007;Gillam et al. 2007). Altogether, most studies examining these echolocation behavioural responses report an upward shift in peak frequency in the presence of conspecifics, rather than downward as has been recorded in our study. Lowering peak frequency when flying with conspecifics is likely multifactorial and requires further study including the application of novel biological modelling (e.g. Mazar and Yovel 2020) to understand the response at different bat densities.
Finally, the peak frequencies of echolocation calls of female adult bats emerging from maternity roosts while newborn juveniles were non-volant were significantly lower than those emitted by adults and juveniles emerging together later in the season. We cannot conclusively determine why this occurs, but it suggests that echolocation plasticity could be linked to bat age and rearing stage, with echolocation calls from adults accounting for a higher proportion of the lower frequency calls we recorded in P. pygmaeus. This frequency shift could be related to social information contained in the echolocation calls for a communication purpose (see review Jones and Siemers 2011) or again, to higher energy demands during pregnancy/lactation which may favour the use of lower frequency calls. Echolocation calls during roost emergence once juveniles started flying were higher perhaps due to adults teaching juveniles the "typical" echolocation structure of the species (e.g. Liu et al. 2007).

Implications for acoustic bat surveys with automatic/manual classifications
As with all acoustic studies of free-flying bats, we can never be completely sure that none of the echolocation recordings in our study area correspond to P. pipistrellus, despite strong evidence that the species does not occur there. Nonetheless, numerous studies have highlighted the importance of a combination of survey methods to accurately assess species diversity in bats (e.g. Kuenzi and Morrison 1998;Flaquer et al. 2007) and to confirm the presence of species identified acoustically via other methods (e.g. Hill and Greenaway 2005). Current EUROBATS guidelines recommend bat surveys with detectors as the primary method for monitoring P. pygmaeus and P. pipistrellus (Battersby 2010). Our study shows that this may not always be suitable and that plasticity in echolocation of P. pygmaeus can be extreme and complex and is not yet accounted for in current guides of call parameters for acoustic identification. We showcase that species distribution maps based solely on acoustic records may be flawed for some cryptic species and ground truthing with complementary survey methods is required.
We showed that both the use of an automated classifier and manual identification resulted in classification of P. pygmaeus sonograms as P. pipistrellus, highlighting a limitation of both methods in accounting for echolocation frequency plasticity. Manual classifications were also incredibly varied, reflecting the subjectivity that is widely recognised as a limitation in manual classifications of acoustic recordings (Russo and Voigt 2016). We therefore emphasise the words of caution from other authors that have highlighted the variability of echolocation calls and the difficulty of species identification from acoustic recordings (Barclay 1999;Barclay and Brigham 2004;Russo et al. 2018). Echolocation plasticity also poses a challenge for the increasingly used automatic classifiers for bat acoustic identification. Extensive field testing is necessary not only to compile reference calls to train automatic classifiers (Vaughan et al. 1997b;Obrist et al. 2004) but also to test their performance in the presence and absence of congeneric species and their ability to account for differences in echolocation call structure and characteristics based on flight situation (Berger-Tal et al. 2008). Importantly, regional differences must be taken into consideration-particularly where there are changes in species composition that affect the frequency partitioning of bat echolocation as they occupy variable acoustic niches.
Further research is required to examine how behavioural echolocation plasticity relates to different factors individually (i.e. density of conspecifics, presence/absence of congeneric species, type of habitat, reproductive status and age) and in different bat cryptic species pairs with similar niches and echolocation calls. We raise an important limitation in the use of acoustics alone for conducting species inventories within open area habitats with high densities of only one species of a cryptic species pair. As the use of acoustic technologies continues to develop and is used more widely in ecological studies, understanding and accounting for echolocation plasticity is vital for researching bat ecology, monitoring species distributions and population trends and implementing effective conservation measures. invaluable guidance and knowledge of the study area. We would also like to thank all landowners and farmers who provided land access and without whom this work would not have been possible. Thank you to all those who contributed to the field surveys and laboratory analysis and to Huma Pearce for the loan of an acoustic lure. We would also like to thank Joseph Tobias, Josh Hodge and two anonymous reviewers for their valuable suggestions on improvements to this manuscript.
Funding This work was funded by the Catalan government and the Granollers Natural Sciences Museum (convention reference DB201804).
Code availability Not applicable.

Declarations
Ethics approval Data used in this study were collected under permits from the Catalan Government Ministry of the Environment (SF/0451/ 2019). Bat capture and handling were conducted following guidelines approved by the American Society of Mammologists (Sikes and Gannon 2011). All applicable international and national guidelines for the surveying, care and handling of bats were followed.

Consent for publication Not applicable.
Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.