Passive acoustic monitoring reveals seasonal patterns in European green toad calling activity but fails to accurately reflect population abundance

Study site

The study sites were located within the urban area of Poznań, a city in western Poland with a human population of over 500,000. The focal species exhibits a patchy distribution throughout the city, with a few relatively large populations thriving in city parks and smaller and/or ephemeral populations occurring in urban wastelands, construction sites, and similar habitats—a pattern also observed in other Central European cities such as Vienna and Košice^27,33. Overall, the number of European green toad breeding sites in Poznań has decreased substantially since the 1990s³⁹. Data were collected at three locations. Site 1 and Site 2 are both located within Cytadela Park—a large urban park established on the grounds of a former Prussian fortress. Site 1 – Rosarium [52°25’26.02″N, 16°55’56.32″E] is a shallow artificial pond with a concrete bottom (approx. 5000 m²), built in the 1970 s, which hosts a large European green toad population that has been extensively monitored over the past decade^{29,30,31,32,33,34,35,36,37,38,39}. Site 2 – Water Garden [52°25’15.48″N, 16°56’22.81″E] is a smaller artificial pond (approx. 400 m²) with a concrete bottom and a fountain, reconstructed in 2019 and subsequently colonized by European green toad. The two sites are located over 600 m apart and calls from one site are inaudible at the other. Site 3 – Hippodrome [52°25’26″N, 16°55’56″E] is a shallow artificial pond (approx. 900 m²) with a geomembrane bottom covered with gravel, created in early 2022 for the purpose of translocating European green toad population from a breeding site scheduled for demolition. The population was translocated in 2022, with mass breeding events occurring in both 2022 and 2023, the latter being the year when the acoustic data were collected.

Soundscape recordings aimed at evaluating the effectiveness of automated acoustic methods and analysing patterns of vocal activity were conducted at all three sites. However, field surveys to assess the relationship between calling intensity and the number of observed males were carried out only at Sites 1 and 2.

Field survey

Breeding activity of European green toad was monitored in 2023 (April 14th – July 2nd) at Site 1 and Site 2. Field surveys were conducted by a single observer. Observations began approximately one hour after dusk (around 21:30 CEST from April to mid-May, 22:00 CEST in late May, and 22:30 CEST in June and July). Each site was surveyed 21 times, with each visit lasting approximately 30 min. During each survey, the observer searched for individuals in the breeding pond using a 120-lumen flashlight, recording the number of males, females, and amplexed pairs. Additionally, vocal activity, presence of egg strings, and tadpoles were noted.

Soundscape recording

At all three locations, we recorded soundscapes using three autonomous sound recorders (Song Meter Mini, Wildlife Acoustics) between April 12 and July 10, 2023. This period should cover the entire breeding season of the European green toad. Each recorder (one per location) was mounted on a tree approximately 5 m above ground level and positioned 10–15 m from the edge of the water body. Recordings were made daily, beginning two hours before sunset and continuing until two hours after sunrise. We recorded 60-minute WAV files at a sampling rate of 22.05 kHz/16-bit, with an 18 dB gain applied. The microphone sensitivity was tested before the start of recording using VOLTCRAFT SLC-100 sound level calibrator.

Manual spectrogram scanning and listening to recordings

In the first step, we selected one-minute sound samples (one minute every 20 min), considering only the days when a field observer surveyed European green toads on site. We analysed the first minute of each 60-minute recording using manual spectrogram scanning and listening to recording (AS) in Raven Pro 1.6.5 software (window: Hann, FFT: 512 points, 50% overlap). Each vocalization of the European green toad was annotated, and we measured its peak frequency, along with the 5th (L5) and 95th (U95) percentiles of the energy distribution (i.e., the frequencies below which 5% and 95% of the acoustic signal’s energy are concentrated, respectively). These manual detections were then used to evaluate the effectiveness of automatic call detection and the performance of acoustic indices.

Automatic detection of European green toad calls

The remaining one-minute sound samples (from the 20th and 40th minute of each hourly recording, totalling 954 samples) were used to develop an algorithm for the automatic detection and classification of European green toad calls using Cluster Analysis in Kaleidoscope Pro 5.6.8 software. We calculated the average spectral parameters of the calls: L5 (1229 ± 114 Hz), peak frequency (1341 ± 110 Hz), and U95 (1409 ± 106 Hz) (Fig. 1). To constrain the recognizer to the frequency range characteristic of the studied species, we set the minimum and maximum frequency limits to 1000 Hz and 1700 Hz, respectively. Additionally, we specified the following parameters: a minimum call duration of 1 s, a maximum duration of 4 s, a maximum inter-syllable gap of 0.1 s, a maximum cluster distance of 1, and an FFT window size of 128. Using these settings, we scanned the recordings to detect and cluster vocalizations. In the next step, each detection was manually reviewed and classified as either a European green toad call or another type of sound. In this way, we removed false positive detections to obtain a more reliable recognizer. Subsequently, we re-scanned the recordings using the manually edited clusters to build the final recognizer.

Spectrogram of a male European green toad call. The call consists of a series of pulses repeated multiple times. In this study, the frequency range for analysis was limited to 1.0–1.7 kHz (grey band), which is characteristic of the target species’ vocalizations. Spectrogram generated using tuner and seewave packages in R.

The final recognizer was then applied to the one-minute sound samples in which all calls had been manually identified. The recognizer was trained on one dataset (20th and 40th minute of each hourly recording) and tested on a separate dataset (1th minute of each hourly recording) to objectively assess its performance on unseen data and to prevent overfitting. To evaluate its performance, we calculated the Pearson correlation coefficient between the number of calls detected automatically and manually in each one-minute sample, along with recall (TP/[TP + FN]) and precision (TP/[TP + FP]). The abbreviations are defined as follows: TP – true positives; FN – false negatives; FP – false positives. Recall represents the proportion of correctly detected vocalizations relative to all vocalizations identified manually. Precision indicates the proportion of correctly detected vocalizations relative to all detections made by the algorithm. A high recall value suggests that the algorithm detected most vocalizations, while a high precision value indicates a low rate of false detections.

Acoustic indices calculations

We calculated two commonly used acoustic indices: the Acoustic Complexity Index (ACI) and the Bioacoustic Index (BI). The ACI quantifies the absolute difference in amplitude between two adjacent time samples within a frequency bin, relative to the total amplitude⁴⁰. Higher ACI values are generally indicative of a more complex soundscape, which is typically associated with greater species diversity and abundance⁴¹. The BI measures the area under the log-amplitude spectrum curve, with the minimum decibel level set to zero. It was originally developed to estimate avian abundance and community composition in Hawaiian ecosystems⁴². More recent studies have shown a positive correlation between BI and both species richness and abundance in other regions^41,43. The selection of these two indices from among the dozens used in ecoacoustics was motivated by the desire to provide a comprehensive description of the acoustic diversity of the soundscape, capturing both the complexity and temporal dynamics of sounds (ACI) as well as their overall intensity (BI).

In our study, we restricted the frequency range for both indices to match that of European green toad calls (1000–1700 Hz) and used a 256-point FFT window. For the ACI, we additionally set the j parameter to 5 and applied an amplitude threshold of −50 dB (Fig. 1).

Seasonal and daily pattern of calling intensity

We calculated the BI and ACI for one-minute sound samples and averaged the values for each full hour, starting two hours before sunrise. Recordings from the end of each session that were shorter than one hour were excluded from the analysis. Simultaneously, we used the number of calls automatically detected for each hour of recording. Additionally, since the sound recorders logged temperature data every minute (built-in temperature logger inside the recorder), we were able to calculate hourly average temperatures. This analysis included soundscape data from all three locations.

Statistics

To examine the relationship between the number of males detected by a human observer at the two study locations and the outputs of an automatic acoustic approach, we conducted two sets of generalized linear mixed models (GLMMs) in glmmTMB package⁴⁴. In the first set, we considered the number of calls detected automatically, along with the average ACI and BI calculated during the second hour after sunset. In the second set, we used the total number of calls and the average ACI and BI values for the entire night. In each model, the dependent variable was the number of males detected by the observer. In this way, we focused both on the specific time when the field surveys were conducted by the observer (approximately the second hour after sunset) and on the entire night during which the field surveys took place. Predictors included the number of calls, ACI, BI, time within the season (in days), and temperature (either the daily minimum or the average during the second hour after sunset). Location ID was specified as a random effect. All models were fitted using a negative binomial distribution with a log link function. The negative binomial distribution was used because the response variable (number of males detected) is count data that exhibited overdispersion. Variance inflation factors in all models were below 5.0.

To analyse daily and seasonal patterns of vocal activity in relation to temperature, we conducted three separate GLMMs. The dependent variables were the average number of calls detected automatically, the BI, and the ACI—each calculated as an hourly average. Predictors included time within the season (as a continuous variable), time relative to sunset (a continuous variable ranging from two hours before sunset to the last full hour of recording, up to two hours after sunrise), and the average hourly temperature. Recording point ID was included as a random effect. We used a negative binomial distribution for the number of calls and a gamma distribution with a log link function for ACI and BI. All statistical analyses were performed in R 4.3.3⁴⁵, and model diagnostics were conducted using the DHARMa package⁴⁶.