Bird population declines and species turnover are changing the acoustic properties of spring soundscapes

Bird data

North America: we used annual bird count data collated under the North American Breeding Bird Survey (NA-BBS: https://www.pwrc.usgs.gov/bbs/) from 1996 to 2017. NA-BBS survey routes, consisting of 50 survey points (hereafter sites) evenly distributed over ~24.5 miles, are distributed across the United States and Canada and are usually surveyed in June. At each site, skilled volunteers conduct a three-minute point count, recording all birds seen or heard within a 400-m radius⁵⁹.

Europe: we used annual bird count data from 23 survey schemes across 22 countries collated under the Pan-European Common Bird Monitoring Scheme (PECBMS: https://pecbms.info) from 1998 to 2018. In each scheme, skilled volunteers carry out either line transects, point counts or territory mapping at survey sites during the breeding season and record all birds encountered⁶⁰ (Supplementary Table 5); while methods vary between survey schemes, they are consistent within schemes across the time period included here.

Where count data were reported for subspecies, these were aggregated to species level and any records of hybrid species or specifying genus only were removed. The longitude and latitude of each survey site (just the first site of each NA-BBS survey route) were also provided by NA-BBS and PECBMS. Not all sites were surveyed in every year and only sites surveyed at least three times during the defined time period were included in analyses. Note that similar results were found when restricting data to sites surveyed in at least 10 years during the defined period.

Sound recordings

Sound files for all species detected on NA-BBS and PECBMS surveys were downloaded from Xeno Canto, an online database of sound recordings of wild birds from around the world (http://www.xeno-canto.org). Specifically, we identified all files longer than 30 s, with associated metadata categorising them as high quality (category “A”) and as either “song”, “call” or “drumming” types; sound files whose type category including the term “wingbeat”, “flap”, “begging”, “alarm” or “night” types were excluded. Sound files downloaded for NA-BBS species were restricted to those recorded in North America and those from PECBMS to recordings made in Europe. If no sound files met these requirements for a given species, we downloaded all files of shorter duration for that species that met the quality and type criteria and stitched repeats of these together to produce files longer than 30 s. Where more than 50 sound files for a given species met our criteria for inclusion, a random selection of 50 was taken for use in subsequent analyses. We used multiple sound files for each species to capture, where possible, between-individual variation in song and call structure, with the sound file(s) for inclusion in specific soundscapes randomly subsampled from this set. If no sound files for a species were available, the sites where that species was detected were removed from subsequent analyses; this represented <1.5% NA-BBS sites and <3.5% PECBMS sites. Each downloaded sound file was then standardised to ensure consistent sampling rate, duration and volume. Each file was clipped to the first 27.5 s, with the first 2.5 s of this then removed to produce a 25 s recording. These sound files varied in the quantity of vocalisation they contained according to the song and call characteristics of the focal species. Thus, some included 25 s of continuous song while others included just a single, short burst of sound. The sampling rate was set to 44.1 kHz, and each file normalised with a −6 dB gain before being saved as a mono mp3 output.

It is important to recognise that the sound recordings used here are taken in the wild and thus inevitably contain some background noise in addition to vocalisations of the target species, and that this may influence the acoustic properties of the constructed soundscapes to some extent. To minimise this, we selected only Quality “A” recordings and clipped out 25 s from the beginning of each of these for use in soundscape construction, on the assumption that the named focal species will be more dominant in these recordings and that it is most likely to be vocalising towards the beginning of a submitted recording. Furthermore, any background noise is expected to be both random in acoustic structure and randomly distributed across the sound files of species considered here; we see no plausible reason why, for example, the field recordings of increasing or declining species would be more or less likely to contain background noise. Our systematic approach to soundscape construction and our analyses of trends in standardised site-level acoustic metrics also limits the potential of background noise to cause directional bias in the results reported and, if anything, it is expected to have reduced our ability to detect changes in soundscape characteristics.

In total, count data were available for 202,737 sites and 620 species in North America, with a mean ± SE of 15.62 ± 0.6 sound files available per species. For Europe, count data were available for 16,524 sites and 447 species, with 21.05 ± 0.9 sound files per species.

Soundscape reconstruction

This is described in detail in the main text.

Soundscape characteristics

Four acoustic indices were used to explore changes in the acoustic properties of reconstructed soundscapes. The Acoustic Diversity Index (ADI) uses the Shannon–Wiener index to estimate acoustic diversity, dividing spectrograms into frequency bands and calculating the proportion of each band occupied by sounds above a set amplitude threshold³⁰. Higher values represent a more even distribution of sound across frequencies and are associated with increased species richness. The Acoustic Evenness Index (AEI) uses a similar approach, dividing spectrograms into frequency bands but using the Gini coefficient to measure the evenness of sound distribution across them³⁰. It is therefore negatively related to ADI, with higher values representing a greater unevenness between frequency bands, suggesting dominance by fewer species. Increases in abundance are expected to have less impact on ADI and AEI than increases in species richness as the songs of individuals from the same species will broadly occupy the same frequency space. The Bioacoustic Index (BI) measures variation in amplitude across a range of frequencies by calculating the dB spectrum across frequencies and quantifying the area under the curve³¹. BI is expected to increase with both increases in abundance and species richness. Total Acoustic Entropy (H) is defined as the product of spatial and temporal entropies and quantifies variation in amplitude across frequency bands and time using Shannon–Wiener index³². It increases with both species richness and abundance following a logarithmic model^28,32. As soundscapes become saturated, the influence of additional species and/or individuals on BI and H is expected to decrease. Default settings were used for each acoustic index except BI, where the maximum frequency was set to 22,050 Hz.

We initially generated soundscapes for a series of simulated communities to confirm that the acoustic indices respond as expected when calculated from artificial soundscapes. Firstly, we calculated ADI, AEI, BI and H for soundscapes derived from communities comprising 1 to 10, 20, 30, 40 or 50 individuals of each species in turn. Given the randomised selection of sound files, insertion point and playback volume, we iterated this process 1000 times for each species-abundance combination. Next, we constructed communities containing 2, 3, 4, 5, 10, 20 or 50 species, with 1–10 individuals of each species present, i.e. 70 communities in total. We iterated this process 100 times for each species richness-abundance combination, randomly selecting species for inclusion from the NA-BBS species pool, and a further 100 times, randomly selecting species from the PECBMS species pool. Again, the four acoustic indices were calculated for each soundscape produced.

Annual soundscapes for each NA-BBS and PECBMS site were constructed from each site-year count file and the four acoustic indices were calculated for each. Given the randomised selection of the specific sound file, insertion point, and playback volume used to represent each individual during the construction of each soundscape, this process was iterated five times, with each acoustic index averaged across these five site-year iterations for use in subsequent analyses. For all PECBMS sites and for the first site of each NA-BBS route, the soundscape generated from the fifth iteration was saved as an .mp3 file. All sound file processing and soundscape construction was undertaken using Sound eXchange programme (SoX: http://sox.sourceforge.net/) and acoustic indices were calculated using R packages ‘seewave’³², ‘soundecology’⁶¹ and ‘tuneR’⁶² in R v3.5.1⁶³.

Finally, we tested the sensitivity of soundscape characteristics to key parameters imposed during construction. While predominately driven by community composition, the acoustic properties of constructed soundscapes could also be influenced by rules that influence the degree of the overlap between individual sound files and their amplitude. First, we generated a community of 10 randomly selected European bird species and specified declines in each species from 10 to five individuals over a 6-year period. For each year, we then constructed four soundscapes and extracted the associated acoustic indices for each. The first soundscape type was built using the methods described above. The second was built by inserting sound files into a 3-min soundscape, to increase the degree of overlap, while the third was built by inserting sound files into a 10-min soundscape to decrease the degree of overlap. Finally, we reverted to a 5-min soundscape but randomly sampled playback volume for each sound file from a half-normal distribution. This increased the relative proportion of distant vocalisations and may be more representative of point count data, where the area surveyed increases with increasing distance; though note this is likely to be offset by reduced detectability at greater distances. This process was iterated for 1000 randomly sampled communities of 10 species.

Statistical analyses

Response of acoustic indices to changes in community structure

To confirm that acoustic indices respond to changes in species richness and abundance, we fitted General Linear Models (GLMs) to outputs for the simulated single and multi-species communities. In each model, the mean acoustic index across all iterations was fitted as the response variable. For the single-species communities, the log number of individuals was fitted as the explanatory variable and for the multi-species communities, the log number of individuals, log number of species and their interaction were fitted as explanatory variables. Separate models were fitted to the North American and European data and for each acoustic index in turn.

Site-level changes in acoustic indices

We standardised each acoustic index within each site (by subtracting the mean site-level measure from the annual value and dividing by the site-level standard deviation⁶⁴) prior to analysis to account for any potential differences in detectability or observer effects between sites, differing sampling protocols across survey schemes, and for initial community structure. In all analyses, separate models were constructed for North American (204,813 sites on 4197 routes spanning 22 years) and European data (16,524 sites spanning 21 years), and for each acoustic index in turn. To explore large-scale temporal trends while accounting for any geographic differences in acoustic characteristics, we fitted Gaussian General Linear Mixed Models (GLMMs) via the R package ‘lme4’⁶⁵. Standardised annual site-level values for each acoustic index were fitted as the response variable, with latitude, longitude and year (continuous) as fixed effects. To account for non-independence of soundscapes from the same site, random effects of site and year were included in all models, along with route and state (North America models, Eq. (1a)) or country (Europe models, Eq. (1b)). To assess the importance of fixed effects, we performed a likelihood ratio test by comparing models with and without a particular term, reporting the χ² value and associated significance. Spatial autocorrelation of modelled residuals was examined by Moran’s I, separately for each year, using the package ‘ape’⁶⁶. While significant spatial autocorrelation was found, the sizes of the estimates were negligible (Supplementary Table 6) and therefore this is subsequently ignored. To explicitly explore how temporal trends in the acoustic properties of reconstructed soundscapes varied geographically, we refitted the models described above, including latitude*year and longitude*year interaction terms. To visualise the large-scale annual variation in acoustic properties we refitted these models with year included as a categorial rather than a continuous variable, with predictions from these models providing continent-level annual estimates for each acoustic index (Fig. 3).

To explore the relationships between site-level trends in each acoustic index, we fitted GLMs with the standardised annual values for each index as the response variable and year (continuous) as the explanatory variable (Eq. (2)). This resulted in an independent estimate of the rate of change in each acoustic index at each site. For all six possible pairwise comparisons between acoustic indices, we used Pearson’s correlation coefficients to estimate the magnitude of the association between their site-level trends. All statistical analyses were carried out in R v3.5.1⁶³.

where i = site, t = year, j = state, k = route

where i = site, t = year, j = country

where t = year

To explore large-scale temporal trends in the total number of individuals and species recorded on NA-BBS and PECBMS surveys, we fitted two additional GLMMs. Standardised annual site-level values of the total number of (a) individuals or (b) species were fitted as response variables, with latitude, longitude and year (continuous) as fixed effects. To account for non-independence in community structure from the same site, random effects of site and year were included in all models, along with route and state (North America models) or country (Europe models). Model structures were therefore equivalent to those set out in Eqs. (1a) and (1b), albeit with different dependent variables. We then refitted these models including year as a categorial rather than a continuous variable to visualise the large-scale annual variation, and used predictions from these models to provide continent-level annual estimates for total abundance and species richness (Supplementary Fig. 5).

To explore the site-level relationships between trends in total number of individuals, total number of species and acoustic indices, we first fitted GLMs with either the standardised total number of (a) individuals or (b) species as response variables and year (continuous) as the explanatory variable at each site. These models were therefore equivalent in structure to that described in Eq. (2) and resulted in independent estimates of the rates of change in the total number of individuals and species at each site. We then fitted separate GLMMs for each acoustic index, in each continent, in turn with site-level trend in acoustic index as the response variable and site-level trends in the total number of individuals and the total number of species and their interaction as fixed effects. State was included as a random effect in the North American models and country as a random effect in the European models. To incorporate the error associated with site-level trend estimates we used a bootstrapping procedure in our assessment of the significance of the modelled effects. We generated 1000 new estimates for each variable (site-level trend in: acoustic index, total number of individuals and total number of species) by randomly sampling from a normal distribution with a mean equal to the site-level trend and standard deviation equal to the standard error of the site-level trend. The GLMMs were then fitted over each of the 1000 datasets separately. We present the results of a final model carried out on the original site-level estimates, as well as the proportion of times each fixed effect included in the final model was significant across the 1000 bootstrapped datasets (p < 0.05). Non-significant interaction terms were removed from the models.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com