Using Approximate Bayesian Computation to infer sex ratios from acoustic data.

Population sex ratios are of high ecological relevance, but are challenging to determine in species lacking conspicuous external cues indicating their sex. Acoustic sexing is an option if vocalizations differ between sexes, but is precluded by overlapping distributions of the values of male and female vocalizations in many species. A method allowing the inference of sex ratios despite such an overlap will therefore greatly increase the information extractable from acoustic data. To meet this demand, we developed a novel approach using Approximate Bayesian Computation (ABC) to infer the sex ratio of populations from acoustic data. Additionally, parameters characterizing the male and female distribution of acoustic values (mean and standard deviation) are inferred. This information is then used to probabilistically assign a sex to a single acoustic signal. We furthermore develop a simpler means of sex ratio estimation based on the exclusion of calls from the overlap zone. Applying our methods to simulated data demonstrates that sex ratio and acoustic parameter characteristics of males and females are reliably inferred by the ABC approach. Applying both the ABC and the exclusion method to empirical datasets (echolocation calls recorded in colonies of lesser horseshoe bats, Rhinolophus hipposideros) provides similar sex ratios as molecular sexing. Our methods aim to facilitate evidence-based conservation, and to benefit scientists investigating ecological or conservation questions related to sex- or group specific behaviour across a wide range of organisms emitting acoustic signals. The developed methodology is non-invasive, low-cost and time-efficient, thus allowing the study of many sites and individuals. We provide an R-script for the easy application of the method and discuss potential future extensions and fields of applications. The script can be easily adapted to account for numerous biological systems by adjusting the type and number of groups to be distinguished (e.g. age, social rank, cryptic species) and the acoustic parameters investigated.


Introduction
The proportion of males (POM) in animal populations is of great interest to ecologists and conservationists, because sex ratios influence mating systems and effective population size, the latter being crucial for the maintenance of genetic diversity [1]. In species whose sex cannot be PLOS  easily identified visually because they are elusive or nocturnal, or lack visible secondary sex dimorphism, other means of determining the POM are required. Molecular sexing can offer an alternative [2,3], but is usually time and cost-intensive, and often invasive due to the challenges arising from non-invasively collected samples, e.g. the need to process replicates of each sample and the associated increase in costs [3,4].
In species with vocal sexual dimorphism, these issues can be circumvented by acoustic sexing. In the last decades, passive acoustic monitoring has become a popular tool for collecting information about biodiversity [5], population densities [6,7], animal movement and behaviour [8], and how these factors are impacted by anthropogenic activities [9]. Both in terrestrial and marine environments, acoustic data have gained in importance [9,10] across taxonomic groups from insects [11], amphibians [12], reptiles [13], fish [14,15], and birds [16] to mammals [17][18][19][20][21][22]. Alongside species identity [21,23,24] acoustic data can encode information about sex [25][26][27][28], individual identity [25,29], body size, age, and social group [28,[30][31][32] or geographic origin [33][34][35]. Sex identification based on calls, for example, has successfully been applied to 25 of 69 investigated bird species lacking external sex dimorphism [36]. Acoustic sexing was however less reliable or even impossible in many of the studied species due to varying degrees of overlap in acoustic characteristics of males and females. Such overlap, which is commonly encountered across taxonomic groups [36,37], strongly limits the scope of the application of acoustic data. Therefore, novel approaches are required that allow acoustic sexing and inferences of sex ratios despite such an overlap.
To investigate the possibility to infer sex ratios and to provide sexing methods from acoustic data, we use Rhinolophus hipposideros (Rhip) as a biological model showing acoustic sexual dimorphism with overlap between sexes [28]. This bat species is of high conservation concern: it is classified as near threatened in the European Union Red List and listed in the annex II and IV of the EU Habitats and Species Directive. Hence population monitoring is legally required. In Rhip, the reproductive output of colonies is usually estimated by dividing the number of juveniles by the number of adult females. The latter however cannot be easily determined from visual counts of adults in the colony, because most bat species, including Rhip, show no visible secondary sex dimorphism. As in many mammalian species, sex identification in bats mostly relies on catching individuals and inspecting their sexual organs. For regular monitoring, this approach is unsuitable due to the large time effort required and handling stress for the animals, especially for vulnerable or threatened species. Current estimates of reproductive output therefore usually assume all adults present in the colony to be females. Males have however been reported to be present in Rhip maternity colonies in the past [38][39][40]. Zarzoso-Lacoste & Jan et al. [3] further demonstrated that the POM in maternity colonies can be substantial and importantly, differ between colonies. Traditional estimates are therefore prone to underestimate the number of offspring per female to an unpredictable degree and therefore call for alternative methods to be developed and tested.
The high calling rates of echolocating bats make them ideal for acoustic monitoring (e.g. [23,41]) and for developing methods of acoustic sexing. However, sexing based on simple call parameters such as the call frequency has not yet been possible, because these acoustic parameters overlap between sexes [28,42,43]. Hence, an approach that would allow a quick and simple determination of sex ratios in free-ranging populations despite such an overlap would greatly benefit monitoring programs and ecological studies. While different mixture model approaches exist to address this situation (cf. section "Screening for and testing available approaches"), we did not encounter an approach allowing reliable inference of the POM for overlapping ranges of acoustic parameters. As a consequence, we use Approximate Bayesian Computation (ABC) to infer sex ratios from passively recorded calls in a species with vocal sex dimorphism, where acoustic value ranges of males and females overlap. The developed novel approach moreover allows the investigation of spatial or temporal changes in sex composition, and the assignment of sex to calls. We additionally develop a simpler approach that only considers calls outside the overlap zone. The performance of all methods is tested with simulated data. We additionally validate the reliability of acoustic sex ratio determination with empirical data (recorded echolocation calls) from four molecularly sexed Rhip colonies.
All methods are implemented in the provided R-scripts (S1 and S2 Files). Users do not need to have programming experience if their study system meets the assumptions of the currently implemented approach, because the adjusted priors are the only input required. More experienced users can change additional settings to relax various assumptions and extend the approach to other study systems. Consequently, our approach has a wide range of applications in ecology and conservation biology.

Screening for and testing available approaches
We assessed the list of R-based cluster analyses and finite mixture modelling approaches available on CRAN (https://cran.r-project.org/web/views/Cluster.html) for their suitability concerning our study question. The maximum-likelihood (ML) based package mclust [44] was considered a promising candidate and its performance was tested (see S3 File for details) with the simulated data described in the section "Simulated data set". mclust performed relatively poorly in determining the POM, especially for low sample sizes, i.e. 100 or 200 calls (S1 Fig). The mclust package was therefore judged unsuitable for reliable estimation of the POM from acoustic data in Rhip colonies, where datasets of fewer than 200 calls are not uncommon.
Other ML approaches were not tested as they did not meet the requirements of ratio inference from univariate, normally distributed data.
A major advantage of Bayesian approaches compared to ML approaches is the specification of priors [45] by the user, which permits the algorithm to draw on existing information, hence potentially greatly increasing performance. In our study system, the value ranges of male and female mean peak frequencies are known, providing information that can be used as priors to improve performance in the estimation of POM. Existing Bayesian clustering approaches however do not allow the specification of uniform, or flat, priors where all possible values within the range are equally likely a priori. To overcome this issue, we developed an ABC approach that incorporates both the required feature of ratio inference from univariate, normally distributed data and the flexibility to define uniform priors.

Inferring sex ratios with ABC
The approach we developed is based on an ABC framework inferring the most likely parameters given some observed distribution by comparing summary statistics of the observed data to summary statistics of simulated datasets. In our case, the main parameter investigated was the POM in a given dataset (set of echolocation call recordings from many individuals). The characteristics of the simulated datasets are described in the section "Simulated data set" below. To generalize our approach and to allow for some uncertainty in prior knowledge in acoustic parameters for both sexes, we included three additional parameters: the mean peak frequency of each sex separately (two parameters) and the standard deviation (sd) of peak frequency of each sex. The latter was assumed to be equal for both sexes, and hence was represented by a single parameter, an assumption that can be relaxed. Based on our own observations (I.Karst, W. Schorcht, M. Biedermann) and published data on our focal species [28,[46][47][48], the following parameter priors (with uniform distribution) were used: POM: 0-1, mean peak frequency of males: 104-107 (kHz), mean peak frequency of females: 108-111 (kHz), sd of peak frequency: 0.8-1.2 kHz. We used two additional constraints stating that the mean peak frequency of males was at least 2.5 kHz, but not more than 5 kHz lower than the mean peak frequency of females [47]. The comparisons between observed and simulated datasets used the following summary statistics a) median of peak frequencies, b) mean of peak frequencies, c) standard deviation of peak frequencies (all irrespective of sex), d) Kolmogorov-Smirnov distance between the observed and simulated distributions of peak frequencies.
We used the adaptive population Monte-Carlo ABC algorithm developed by Lenormand et al. [49]. This algorithm was preferred over others as it was designed to minimise the number of models necessary to reach a given quality of the posterior distribution and was the only one providing reasonably good estimates when sex ratios were severely imbalanced (data not shown). The ABC algorithm was implemented via the EasyABC package [50] in R version 3.4.2 [51]. Unless otherwise stated, the input arguments used for the algorithm were 1000 for the initial number of simulations (nb_simul), 0.4 for the proportion α of best-fit simulations to update the tolerance level at each step (alpha), and 0.01 as the stopping criterion (p_acc_min). The 95% highest density intervals (HDI) of the estimates obtained via the ABC approach were computed using the 'hdi' function from the 'HDInterval' package [52].

Simulated data set
Simulations used in the testing phase of the ABC algorithm and to perform the ABC analyses were implemented via a custom R script (S1 and S2 Files). For a given total sample size (here, number of calls) and POM, we calculated the number of data points (calls) from females and males, respectively, and simulated normally distributed peak frequency data for each sex separately. We assumed the ratio of male to female calls to be equal to the ratio between the number of males and the number of females. The samples sizes used were 100, 500, 1000, 2500, 5000, and 10000; the POM ranged from 0 to 1 in incremental steps of 0.05, with 100 replicates per parameter combination. The mean peak frequency was either set by the user to create test datasets or was automatically chosen by the ABC algorithm. For the test datasets, the following parameter values were used: male mean peak frequency = 106 kHz, female mean peak frequency = 109 kHz, and sd of peak frequency = 1 for each sex. When running the ABC algorithm, priors were used as defined in the above section. The normal distribution was used as it fitted well to the distribution of peak frequencies observed in Rhip [47].

Inferring sex ratios by excluding the overlap zone
Mirroring approaches from species without overlapping acoustic value ranges, we developed a second method determining sex ratios that solely considered calls that can be classified with high confidence as male or female. The POM was then simply the proportion of male calls in the dataset. Calls falling within the overlap range were therefore excluded in this method (exclusion method). To investigate trade-offs between confidence levels and reductions in sample size, we tested the performance of this method using different thresholds for the exclusion zone. For simplicity, we always considered males to call lower than females, so that the exclusion zone encompassed the right tail of the distribution of male calls and the left tail of the distribution of female calls. In the most stringent case (widest exclusion zone), the lower boundary of the exclusion zone was set to the frequency that was lower than the frequency of 99.9% of all female calls, and the higher boundary to the frequency that was higher than the frequency of 99.9% of male calls (Fig 1). The less stringent approaches used thresholds of 99 and 95%.
This exclusion method used normally distributed data with the mean of the acoustic parameter for each sex and the sd given by the ABC estimation. It can also be applied without using the ABC method but then the value of these parameters must be provided by the user based on prior knowledge of the system. We therefore examined how the accuracy of the assumed parameter values affected the estimated POM. For this purpose, we explored the effect of assuming mean peak frequency values that were up to 3 kHz higher or lower than the true (simulated) value. For simplicity, we did not investigate errors in the assumed sd of peak frequencies and the difference in mean peak frequency between sexes, hence our calculations are most likely underestimating the true error. The 95, 99 and 99.9% thresholds were calculated as explained above based on the user provided values of mean peak frequencies of males and females. Calculations were carried out for male proportions between 0 and 1 in incremental steps of 0.05.

Methods' performance
The performance of the different methods in estimating the POM was evaluated by computing the root-mean-squared error (RMSE): where p i is the estimate of the simulated (≙true) proportion of males P for the ith data set (i = 1, 2, . . ., R).

Two-step approach for small data sets
As the precision of estimates obtained with the ABC method diminished with decreasing sample size (see Results), precise sex ratio estimation might be precluded for small datasets (i.e. low numbers of calls). This applies to studies that are interested in temporal or spatial changes in the POM, e.g. over a month or between different habitat patches, thus requiring estimates on a daily basis or for restricted spatial areas. Using simulated data, we explored the use of a two-step ABC procedure whereby we first estimated the acoustic parameters for each sex (mean and sd) across the whole dataset (e.g. for a whole month). A second step used the 95% highest posterior density interval of these estimates as new priors to run the ABC procedure for a subset of these data (e.g. each day separately). This two-step approach assumes that the sex-specific mean of acoustic parameters does not change in time/space within the dataset considered, i.e. the mean peak frequencies and sd do not change for example from day to day (time) or between two foraging patches (space). We investigated the performance of the twostep procedure using simulated datasets mimicking a POM increase over time in steps of 0.05 from 0 to 1 (≙21 steps, each representing a certain time period). The number of calls per period was kept constant. Simulations were carried out for 25, 50 and 100 calls per period so that the total number of calls per dataset was 525, 1050 and 2100, respectively. A total of 100 simulated datasets was used per parameter combination.

Inferring the sex per recording
The ABC framework does not provide information about the sex of each data point (here, call). Nevertheless, the information on acoustic parameters and the POM it provides can be subsequently used to probabilistically assign a sex to each data point. To do so, we calculated the likelihood ratio of a call being from a male versus a female as the conditional probability of a call being from a male given its frequency divided by the conditional probability of a call being from a female given its frequency. Following Bayes' theorem, the conditional probability of observing a male given the peak frequency can be written as: The conditional probability of observing a female given the peak frequency can be written as: The likelihood ratio (LR) can thus be reduced to: Calls with a likelihood ratio >1 or <1 were considered as being from males or females, respectively.
Using simulated datasets, we investigated the proportion of calls assigned to the correct sex using this probabilistic approach. Simulated datasets and simulation parameters were as detailed in 'Simulated data set' and the POM ranged from 0.1 to 0.9 in steps of 0.1. To investigate the performance of the approach with varying levels of overlap between male and female call distributions, differences of 1, 2, 3, 4, and 5 kHz between male and female mean peak frequencies were considered. Sample sizes of 25, 50, 100, 500 and 1000 calls were investigated with 1000 replicates per parameter combination (225,000 simulated datasets).

Empirical dataset
The empirical data set consisted of Rhip echolocation calls recorded in four maternity colonies in Thuringia, Central Germany, in the summers of 2015 and 2016 (see S1 and S2 Tables for details). Accessing the roosts was approved by local nature conservation authorities (permit Jena AV09_AGO7_17). One automatic recording device (Anabat SD2 bat detector, Titley Scientific) was positioned inside each of the four studied roosts and directed towards the entrance (ca. 3 m away) to record bats entering between 1:00 am and 6:00 a.m. using a zero-crossing division ratio of eight. We used the Anabat CFCread program (Titley Scientific) to split the continuous data into recordings, i.e. sequences of calls considered to be emitted by one animal passing the detector (for details see S4 File: What is a recording in our study?). Recordings were then filtered via the software AnaLook version 4.2.n to include only Rhip calls (S5 File). Subsequently, calls with an average frequency (Fmean) below 100 kHz were removed as they were outside the value range of the call parameter for adult Rhip and to exclude social calls, which have a lower frequency [53]. Fmean and sd were calculated for each recording. The mean of Fmean calculated over all calls for each recording was considered representative of the acoustic parameter (i.e. peak frequency) of the emitter.

Molecular sexing
Bat faeces were non-invasively collected in the studied colonies during the summers of 2015 and 2016 (see Table 1). Newspapers were spread on the ground beneath the roosting sites and faeces were collected approximately 10-20 days later (Table 1). Upon collection, faeces were stored in plastic boxes (1 per colony) pre-filled to 1/3 with absorbing silica-gel beads to prevent DNA degradation until analysis [54,55]. Animals were not touched during the sampling, and accessing the roosts was approved by local nature conservation authorities (permit Jena AV09_AGO7_17). DNA extraction and amplification, as well as multilocus genotyping were carried out as described by Zarzoso-Lacoste & Jan et al. [3], but employing centrifugation instead of a vacuum pump during DNA extraction (see S6 File for a detailed protocol). We also used the same bioinformatics pipeline, with slight modifications: In the current study, all relative fluorescence unit (RFU) peaks corresponding to sex specific alleles were visually reinspected and validated by plotting the RFUs for the appropriate marker with the aid of the Rpackage Fragman [56]. Furthermore, we applied a different rule for the peak of the Y-linked allele to be accepted: it was counted only if the corresponding peak was present in all three replicates and higher than the peak of the X-linked allele in at least one replicate. An exception was made in the rare case where a multilocus genotype (MLG) was found to be identical to Table 1. Proportion of males and peak frequencies of the studied colonies. The proportion of males was estimated with the genetic (Gen.), acoustic ABC (ABC), and acoustic 99.9% exclusion methods. The 95% highest density interval (HDI) is presented for parameters estimated via the ABC approach. Peak frequencies were estimated with the acoustic ABC approach.

Samples
Proportion of males ABC estimate of the mean peak frequency in kHz another MLG in all loci but the sex specific one. In that case, the Y-linked allele was also counted if it was missing in one replicate and if the corresponding peak was not higher than the peak of the X-linked allele. In general, only MLGs with missing data at no more than two loci were kept, provided that the peaks at the sex specific locus were unambiguous.

Estimation of acoustic parameters and sex ratios with ABC
The estimated male and female mean peak frequencies matched the simulated values. Inaccurate estimates were obtained however if the proportion of the focus sex was zero (Figs 2 and 3). Both precision and accuracy increased with sample size (Fig 3). The estimated sd of call frequencies was slightly underestimated when only one sex was present (Figs 2 and 3). Generally, there was a slight upward bias in the estimated sd that disappeared with increasing sample size (Fig 3). POM values estimated via ABC showed a high degree of concordance with the real values in the simulated data set for all sample sizes (Fig 4). Both precision and accuracy increased with sample size (Figs 4 and 5). The RMSE did not exceed 0.05 for the lowest sample size (100 calls), and dropped to less than 0.02 when increasing the number of calls to 500 (Fig 5). Further increases in the number of calls resulted in even lower RMSEs for all sample sizes tested, but beyond 5000 calls the improvement became marginal (Fig 5).
For the four sampled Rhip colonies (empirical datasets), mean peak frequencies were inferred to be between 105.0 and 105.9 for males, and between 108.3 and 109.1 for females ( Table 1). The POM estimates ranged between 0.21 and 0.63 depending on the colony, differing from the genetically determined one by 0 to 0.11.

Two-step ABC approach
Applying a two-step ABC procedure improved the performance of POM estimates compared to an approach where acoustic parameters (male and female peak frequencies and sd) were estimated separately for each subset (Fig 6). RMSEs below 0.04 were obtained for datasets of 50 calls, which matched the performance of the one-step approach but with a sample size reduced by half (Fig 6).
The proportion of calls assigned to the correct sex using the probabilistic approach was influenced by the difference in mean peak frequency between sexes, the number of calls as well as the POM (Table 2; Fig 7). For minor overlaps between sexes (i.e. 4-5 kHz difference in mean), the method performed extremely well (> 98% calls correctly assigned; Table 2). When the overlap increased, the method still performed relatively well, with e.g. 78% calls correctly assigned when the overlap zone comprises 95% of calls in the dataset (1 kHz difference, N = 1000 calls; Table 2). The variance across datasets, but not the mean correct assignment, decreased with increasing number of calls. The POM affected the results in two ways; first, the more common sex had a larger proportion of correct assignment (Fig 7); second, the overall proportion of correct assignments was higher for more unbalanced sex ratios (Fig 8).

Estimation of sex ratios via the exclusion method
Sex ratios estimated via the exclusion method varied in their performance in relation to the thresholds considered for the overlap zone. The 99% and particularly 95% approach performed as well as the ABC method for rather balanced sex ratios only, but resulted in higher RMSEs for unbalanced ones (Fig 5). A threshold of 99.9% yielded similar RMSEs as the ABC approach (Fig 5). However, the exclusion method was highly sensitive to the accuracy of the prior information, i.e. assumed mean peak frequencies of males and females (Fig 9). The sensitivity was dependant on the difference between the true and assumed mean peak frequency, the sign of this difference and the true POM.
The POM estimated via the 99.9% exclusion method in the four sampled Rhip colonies was between 0.19 and 0.67 depending on the colony. These estimates differed from the corresponding ABC results by 0.01 to 0.04, and from the genetically determined sex ratios by 0.02 to 0.12 (Table 1).

Discussion
We have developed and tested two methods to infer the sex ratio of groups of vocalizing animals where acoustic parameters of the signals overlap between males and females. The newly developed Bayesian (ABC) approach infers male and female mean peak frequencies and the according sd as well as the proportion of males (POM). These inferences can then be used to separately analyse subsets of the data that correspond e.g. to certain time periods of interest and/or to probabilistically assign a sex to single data points (e.g. calls). Finally, the inferred acoustic parameters can be used to define the exclusion zone for a second, simpler approach based on the exclusion of calls from the overlap zone.   could also be investigated if they differ (at least on average) for some acoustic parameter, which is often the case [61][62][63][64][65][66][67]. Our ABC approach can be used to determine the values which are characteristic for vocalizations of subgroups and to estimate their ratio in abundance. Even though the method has been developed for obtaining ratios for two groups, it can be extended to any number suitable for the study question.

Sex ratio estimation
We confirm the robustness of the ABC and to some extent the exclusion method for sex ratio estimation with simulated data sets of varying sample size. Obtaining larger numbers of calls clearly improves precision and accuracy, but even for as few as 100 calls, errors are minor for the Bayesian approach (Fig 5). For our simulated datasets, the 99.9% exclusion method is slightly less accurate than the Bayesian approach for small sample sizes, but performs better compared to the 95% and 99% exclusion approaches (Fig 5). The exclusion method is very straightforward, providing an advantage for practitioners (easily applicable and very quick). However, the mean peak frequencies of males and females (and their sd) used in this approach must be known with very high accuracy to obtain reliable POMs (Fig 9). If this is not the case, we strongly recommend using the ABC method instead. We additionally tested the ABC and exclusion methods with empirical data. For three of the four Rhip colonies, the POM determined from recorded calls is nearly equal to the ratio obtained via molecular sexing (Thu22, Thu26, Thu47; Table 1). These results provide clear evidence that the reliable estimation of the POM from acoustic recordings demonstrated via  Fig 7) and sexes are combined. simulation also applies to empirical datasets. Our method thus can be used to obtain sex ratios and hence, the number of females, providing more reliable information to calculate the number of offspring per female in Rhip colonies from passive acoustic monitoring data. An overestimation of the POM via ABC compared to molecular sexing was observed in the remaining colony (Thu35; Table 1). Theoretically, the molecular sexing approach could suffer from a bias introduced by different sampling probabilities of males and females, e.g. because they prefer different hanging sites or due to differences in roost fidelity [3]. Due to the very good match between molecular and acoustic sexing for the other three colonies we however rather suspect an underrepresentation of some individuals in the recorded calls. The roost of this colony has two entrances, but only one was acoustically monitored. Hence, if more females than males preferred the unmonitored entrance over the one where recordings were made, their calls might be missing or less frequent in the data set, resulting in an overestimation of males. Therefore, we strongly recommend considering sex-specific or individual behavioural differences in the study design of future applications. Importantly, our method assumes a similar acoustic detection probability for individuals belonging to different groups of interest. If for example males call louder or more often than females, they will be overrepresented in the dataset and their proportion will be overestimated. If such differences in detection probability are well quantified for the study species/context, they can be corrected for when simulating datasets within the ABC method. For our target species, the calls are used for orientation and hence likely to be emitted equally often and at similar intensity by individuals of both sexes. Nevertheless, for a successful application of our method to other contexts or species, similarities in signal detection probability should be considered and discussed.  (99.9% threshold). The x-axis shows the difference between the assumed and true mean peak frequency of male calls. Green lines represent an overestimation of the proportion of males, blue lines an underestimation. The dashed, broken, and dotted lines depict sex ratios estimated for different assumed male mean peak frequencies when the simulated POM is 0, 0.5, or 1, respectively.

Two-step sex ratio determination for data subsets
We have developed and tested a method to track changes in the POM via a two-step ABC approach. With this method, sex ratios can be estimated quite reliably even for relatively low numbers of acoustic signals, e.g. short time periods such as days, or small spatial scales such as single foraging/commuting sites. This allows the detection of temporal or spatial patterns in the POM, which could provide interesting insights into the flexibility of social organisation and behavioural differences of males and females in habitat use [68]. The method could be used to compare sex ratios obtained for groups of migrating animals passing a recorder station at different time periods (temporal) or to compare data simultaneously collected from different stations (spatial), for example. Migratory birds are a particularly promising target for such an approach, as acoustic signals even of flying flocks can be used for sex determination in certain species [69]. Furthermore, existing acoustic data from monitoring programs (e.g. [70]) could be used to test for within-species (or between-species, in the case of cryptic species) temporal/spatial segregation or habitat preferences. Similarly to the estimation of acoustic parameters, this temporal/spatial resolution method can also be applied to biological groups other than sex. To ensure accurate estimates, the choice of the temporal/spatial scale will depend on (1) the number of calls that can be obtained and, (2) the confidence in the dataset to meet the assumptions of the method (i.e. the mean of acoustic parameters for each group remaining constant through time/space).

Inferring the sex per recording
We have additionally established a method to probabilistically assign a sex to single acoustic signals based on the acoustic parameters and the POM inferred from the ABC approach. While not error-free at the individual level, our method provides a good strategy for detecting overall patterns of sex-specific behaviour whose temporal/spatial resolution exceeds by far even that of the two-step Bayesian method. As the sex assignment of a single acoustic signal is based on the overall probability of signals being emitted by a male or female within the studied group, which in turn depends on the sex ratio, the two-step ABC method should be used prior to assigning a sex to single signals to confirm temporal/spatial homogeneity in sex composition, or to identify appropriate time periods/spatial scales over which the POM is stable in the population.

Technical considerations for future applications
We have developed a method to infer the relative ratio in abundance of two (or potentially more) groups of animals with partially overlapping value ranges of acoustic parameters. The choice of parameters and distributions used in our investigation was driven by the biological model used to empirically test the method. Our work therefore only covers a small range of parameter combinations compared to what can be encountered in empirical datasets across taxa. For example, normally distributed data are ubiquitous in biology and we expect many datasets to conform to this distribution. Nevertheless, the use of other distributions is very straightforward within the proposed framework which can therefore accommodate datasets with a wide spectrum of distributions. The distributions can be any probability distribution with a known mathematical function (e.g. normal) but also any empirical distribution obtained by fitting an empirical dataset, offering within a single tool greater flexibility than with conventional maximum likelihood or conventional Bayesian models. The range of estimated parameters used in the ABC approach can be customised depending on the dataset. In our case, we used the same sd for male and female peak frequencies but these can be decoupled if needed. Similarly, we selected four summary statistics to compare the simulated and the observed dataset in the ABC procedure but different statistics could be used. We chose the mean, median and sd of peak frequencies (irrespective of sex) and the Kolmogorov-Smirnov statistic between the observed and simulated dataset as these are summarising well the variations observed in different distributions. However, the summary statistics can be easily and quickly customised to improve performance if needed.
When developing the algorithm, Lenormand et al. [49] found that smaller alpha and p_acc_min improved the quality of the posterior approximation but also increased the number of runs and hence time for completing the analysis. The authors recommended to use alpha = 0.5 and p_acc_min between 0.01 and 0.05 depending on the desired level of convergence. Our simulations show a limited influence of the values chosen for those parameters, probably because of the simplicity of our datasets and quick convergence. Nevertheless, it remains important to test the influence of these parameters and select those providing the most accurate estimates based on simulated datasets mimicking the empirical datasets of interest. More generally, we recommend that prior to studying sex ratio differences in empirical datasets, simulations mimicking those empirical datasets should be carried out to investigate the reliability and limitations of the method.
Based on previous investigations of lesser horseshoe bats returning to their roost from foraging trips (data not shown), we chose the spatial arrangement of the recording device and the temporal resolution of the recordings in a manner that limited the probability of two or more individuals being recorded simultaneously. While the simultaneous recording of multiple individuals cannot be ruled out, the good agreement between the ABC estimates of POM and the independent non-invasive genetic estimates suggests that those occurrences are not problematic in our empirical dataset. This suggests that when two animals are unlikely to be recorded simultaneously, one could simply use the mean peak frequency per recording (as done in our study). An alternative would be to filter out those recordings or eventually, one could use the peak frequency of each call instead of the mean peak frequency per recording. We did not investigate this strategy in our empirical dataset because the prevalence of recordings with high sd (suggesting simultaneous recording of more than one individual) was low (data not shown). However, this issue should be investigated in other organisms/situations where such recordings might be common, and especially when the acoustic parameters of calls of the emitter are altered in the presence of other individuals [71,72].
Bats were recorded when entering the roost. Exiting bats were not recorded to avoid multiple recordings per individual introduced by individuals going in and out or circling around the entrance. Small entrances maximise the probability that individuals approach the recording device from similar angles and at similar speed, leading to comparable Doppler shift on all recordings. Setting up the recording device in a more open environment (e.g. foraging grounds or commuting routes) might result in individuals approaching from different angles and recordings being obtained from animals either approaching or going away from the microphone, possibly leading to more ambiguous datasets. Although not tested here, choosing a microphone with high directionality might limit the problem though at the cost of less recordings being collected. Other non-mutually exclusive strategies might be to use the minimum value of the peak frequency per recording (instead of the mean), which is likely to represent the call emitted when the individual is going away from the microphone. Datasets collected in different situations might benefit from different data pre-processing steps prior to sex ratio estimation and we encourage users to explore different options (cf. also the above paragraph). Filtering should also be performed when different types of acoustic signals (e.g. echolocation versus social calls) are present in a dataset. Furthermore, in its current form, the developed ABC approach is based on the assumption that the detection probability is equal for both groups whose ratio shall be inferred. This deserves special attention because noncompliance with this assumption could lead to biases.

Conclusions
The methods developed here will allow scientists and applied conservationists to investigate ecological questions dealing with the specific behaviour of groups of individuals in greater detail, provided that acoustic di-/polymorphism exists between the focal groups. These can comprise individuals that differ in an externally cryptic trait of interest, e.g. sex, species affiliation, age cohort, size or weight class, or social rank. We present a toolbox that combines the inference of acoustic parameters and sex ratio via ABC with downstream applications. The latter can be used to detect changes in group composition over time/space, or to assign single acoustic signals to one of the groups. To allow a broad field of application, we describe the conditions that must be fulfilled in order to apply our approach to other study systems and provide suggestions on how to overcome some challenges that may arise. The approach is based on acoustic data, which in many species can be acquired more easily than close-up morphological or genetic data, thus providing a non-invasive, time-efficient, and relatively low-cost approach to explore ecological traits that differ between groups of animals.