Evidence of garden bird feeding industry change

Garden bird feeding industry data were derived from advertisements featured in Birds magazine; a free publication, widely distributed to members of the Royal Society for the Protection of Birds (RSPB; recent membership figures totalled over 1.2 million49). Our assessment focused on all advertisements promoting garden bird food and/or feeder products in the Autumn editions of Birds from 1973 to 2005. After 2005, advertising was biased toward RSPB-branded products. Over the 33-year period considered, the number of pages featuring all forms of advertising in Birds increased linearly (GLS \(\chi _1^2\) = 26.19, p < 0.001), correlated with a general increase in magazine length (r = 0.901, n = 33, p < 0.001). For every advertisement (n = 179), we extracted data for its size (proportion of page covered), and the names and descriptions of individual food and feeder products. Foods were also allocated to one of 21 different food type categories (see Supplementary Table 1) using product descriptions, images and online information.

To test for temporal changes in supplementary food quantites, the food and feeder product ranges from all advertisements were summed (respectively) per magazine (n = 33). To test for temporal changes in supplementary food diversity, we quantified the number and diversity (using Simpson’s diversity index) of food types represented in each magazine. We ruled out the possibility that observed patterns were confounded by increased advertising space by controlling for the total number of pages featuring bird feeding industry advertisements (advertising pages). Number of bird food products and number of feeder products were modelled seperately as a function of a year smooth using generalised linear models (GAMs), with a Poisson error distribution. Advertising pages was also log-transformed and included as an offset, since the numbers of products advertised was expected to be proporitional the the total amount fo advertising space available. The number of different food types and food diversity over time were modelled as a function of year using generalised least squares (GLS) with advertising pages included as a covariate. Food type number was log-transformed to reach normality of the residuals. In a further test for overall growth in the bird food industry, we modelled the proportion of total advertising space dedicated to feeding products as a function of a year smooth using a GAM. Full details about model specifications are given below.

Garden Bird Feeding Survey

Data from the Garden Bird Feeding Survey (GBFS) were used to investigate changes in bird communities at feeders in mainland Britain between winter 1973/4 and winter 2012/13 (40 years). GBFS is an annual survey monitoring the number of birds visiting feeders in domestic gardens over a 26-week period from October to March (www.bto.org/gbfs). The survey comprises an average of 217.8 (s.e.m. 7.1) gardens each year, covering a representative range of garden types (suburban/urban and rural) and a consistent geographic distribution. Participants leaving the scheme are replaced with new volunteers from the same region, and with gardens of a similar type and size.

Survey participants record the maximum number of each bird species observed simultaneously using feeders (i.e. at/on feeders and in their vicinity) each week or, in the case of sparrowhawk, feeding on birds using feeders. Data for all species known to occur in Britain according to the British Bird List22, except vagrants (n = 6 species), were used to estimate community indices (n = 133; Supplementary Table 2). Scarce migrants, summer migrants (recorded at the beginning or end of the winter period) and species not traditionally considered as garden visitors (e.g. wetland birds) were retained to avoid removing evidence of community change. Records for domestic and aviary species were excluded (n = 21). Species with ambiguous identification, such as marsh tit and willow tit, were combined. Changes in community composition are unlikely to have been biased by long-term changes in the arrival and departure of British breeding migrants, since estimated phenological shifts are not large enough to have noticibly increased the probablity of these migrants being recorded by GBFS50. Similarly, although volunteer field-skills are not formally controlled, there is no reason to suspect spatial biases or temporal change in the accuracy of species’ identification and counts.

Since diversity estimates are influenced by sampling effort51, the data were restricted to gardens with at least 20 weekly submissions for a given winter (mean = 25.13 weeks), as this number of replicates was seen to produce reliable species richness and species-specific abundance measures. More specifically, species accumulation curves, averaged across gardens each winter, reached an asymptote with 20 weeks of surveying. Further, we used general linear mixed models (GLMMs) controlling for garden and year random effects to verify that increasing sampling effort over 20 weeks had no effect on winter abundance (defined as the maximum count observed in a garden across all weeks surveyed). Winter abundance was independent of sampling effort in 94% of species (n = 125/133), significantly more than expected by chance (z-test, χ2 = 101.2, p < 0.001, 95% CI = 89.2 – 100.0%). We note that maximum counts provide an accurate (asymptotic) means of comparing changes in relative abundances across species, but probably under-represent true abundance.

The filtered data set used to conduct the analyses comprised 1,001 GBFS gardens in mainland Britain (mean ± s.e.m. per year = 185.8 ± 6.1) that contributed a total of 186,825 weekly submissions over the 40-year survey period (Supplementary Fig. 1). We do not expect the data set to contain any biases that might influence the observed findings, given that the data were collected across a consistent set of gardens using structured sampling to a defined protocol, and that the species list has been carefully validated and sampling effort controlled. Any erroneous data, for example due to species misidentification, should be a source of noise rather than bias.

Potential drivers of bird community change

To estimate the potential for species’ spatial range shifts to impact measures of community temporal change, distribution data from the 1981/82–1983/84 Winter Atlas25 were compared with equivalent data from Bird Atlas 2007–201124 to identify areas where a species had apparently colonised or apparently disappeared at a 10-km resolution. Data were available for 98% (n = 130) of all species represented in the GBFS data set (Supplementary Table 2); three species were not recorded during either winter atlas period. GBFS gardens were assigned data from the 10-km squares in which they were located, allowing the net change in the number of GBFS gardens located within a species’ range to be calculated and averaged across all species (using the median due to data skew).

In addition to examining bird community changes through time, three potential drivers of garden bird feeder use were also considered: number of feeders, winter temperature18 and local habitat28. The numbers of hanging feeders, bird tables and ground feeding stations (collectively termed feeders) were recorded each week in surveyed gardens. Weekly feeder numbers were averaged annually by feeder type and in total, to provide an indicator of supplementary food availability throughout winter in each garden. To examine changes in the provision of garden bird feeders over time, numbers of feeders, in total and by feeder type, were log-transformed and fitted separately against year using linear mixed models (LMM) with garden identity included as a random effect. The log-linear functional relationship, which specifies the percentage change in numbers of feeders over time, was applied to achieve normality of the residuals.

Gardens were classified as either suburban/urban or rural according to their surrounding habitat. Suburban/urban includes gardens in areas with a mix of built cover and green space, or in dense urban areas with little vegetation, such as town centres. Rural includes gardens in areas away from towns, with just a few scattered houses, farms or other isolated buildings. The difference in garden habitat types was verified using Land Cover Map 199052, which showed that rural gardens were located in 1-km squares with significantly less urban cover on average than suburban/urban gardens (rural = 12.04% ± 0.70 s.e.m.; suburban/urban = 42.66% ± 1.06; χ2 = 40265, p < 0.001). Gardens were re-classified from rural to suburban/urban if urban encroachment occurred (n = 6 gardens), but garden identity remained unchanged.

We expected that weather conditions throughout the whole bird data collection period would have a greater influence feeder use than extreme weather events. Therefore, annual measures of average winter temperature were estimated using mean monthly temperature for October – March and used to test for climatic effects on feeder use. Mean monthly temperature (°C) data were extracted from the UK Meteorological Office Climate Projections (UKCP09) 5 × 5-km resolution gridded data set53. Gardens were assigned averaged winter temperature data for the 5-km square in which they were located.

Evidence of bird community change

We used annual measures of species richness, Simpson’s diversity index and k-dominance to examine bird community patterns. Species richness was the total number of species observed throughout the winter. Simpson’s diversity index was used to provide a robust and meaningful measure of community diversity per winter51. Since Simpson’s diversity incorporates both the number of species present and their relative abundances, its comparison with species richness was also used to infer changes in bird community evenness. k-dominance curves—which plot the cumulative abundances of all species in a community (as percentages) against their species rank (logged)—were used to study changes in community evenness over time51. We used species abundance and rank, averaged annually across gardens, to compare k-dominance curves from each year of the time-series. The higher the curve, the less diverse and more uneven the community it represented.

To estimate national indices of species richness and Simpson’s diversity, we compiled data from all gardens into a single time-series, then applied sample-based rarefaction to standardise sampling effort through time54. Specifically, 115 gardens (equivalent to the minimum number surveyed in a single year) were randomly resampled without replacement from the total pool of gardens surveyed per year to achieve a consistent sample size over time. For each year, data from all resampled gardens were pooled and species richness and Simpson’s diversity calculated. Resampling was repeated for 1000 iterations and diversity measures averaged. Confidence intervals were not generated, since estimates derived from rarefaction are dependent on the size of the subsample and are therefore not informative about sample variability. The rarefied measures of Simpson’s diversity and species richness were modelled separately to quantify national-scale bird community trends. Simpson’s diversity was fitted against a year smooth using a GAM, and species richness was fitted against year using GLS regression.

To assess bird community trends at the garden scale, annual measures of species richness and Simpson’s diversity per garden were fitted using generalised additive mixed models (GAMMs) with garden identity included as a random term (n = 7433 garden-years). To evaluate the influence of other garden use drivers on bird community change, we also included number of feeders, winter temperature and habitat as fixed effects. These terms were standardised to a mean of 0 and s.d. of 0.5 to enable effect sizes to be compared directly55.

Linking feeder use to national population change

Since feeder use would need to be reasonably prevalent within a population to incur national-scale impacts, we focused on species that regularly used feeders when testing for associations between changing feeder use and changes to population size. Data were combined across all gardens per year to derive a single, intuitive index of overall feeder use per species in Britain, defined as the proportion of GBFS gardens in which a species was observed using feeders. Using site occupancy to derive feeder use, as opposed to species abundance, produces an easily interpreted measure of the scale of feeder use nationally, while also minimising the influence of stochastic variation in species counts. We conservatively defined species that regularly used feeders (feeder-users) as those with a mean feeder use of ≥ 0.1 (i.e., observed in an average of 10% of surveyed gardens per year across the study period; n = 39 species; Supplementary Table 3).

For all feeder-users, a binomial generalised linear model, testing the difference in feeder use between the first and last three years in the time-series (1973/4–1975/6 vs. 2010/11–2012/13), was used to estimate the value and significance of net change in feeder use. This approach was used as it could produce a measure of change that was analogous with estimates of national breeding population change over the same timeframe, while also minimising any influence of inter-annual stochasticity and avoiding assumptions about the shape of the temporal trend. More specifically, breeding population changes for feeder-users were similarly calculated as the difference between smoothed annual indices for 1974 and 2013 (i.e. the breeding seasons immediately following the beginning and end of the time-series), derived from the joint Common Bird Census/Breeding Bird Survey (CBC/BBS) trends for England30. CBC and BBS use structured, stratified protocols to monitor national bird populations and inform the UK Biodiversity Indicators. Trends were not available for eight feeder-user species, which included winter migrants and species not well covered by the CBC.

Phylogenetic generalised least squares regression (PGLS) was used to test the relationship between changes in feeder use and national population size (n = 31 species) while accounting for phylogenetic structure. Bird phylogeny was based on a pruned consensus tree produced by majority rules using 100 phylogenetic trees randomly extracted from the avian phylogenies developed by Jetz et al.56 (Supplementary Fig. 5). Within the pruned tree, Eurasian nuthatch (Sitta europaea) was represented by phylogenetically similar white-tailed nuthatch (S. himalyensis)57 and lesser redpoll (Carduelis cabaret) by common redpoll (C. flammea)58 since these species were absent from the global avian tree. Maximum likelihood was used to estimate the PGLS model’s Pagel’s lambda, giving a measure of the phylogenetic covariation between the predictor and response. Pagel’s lambda values of zero indicate that the predictor-response relationship is unrelated to phylogeny, whereas high lambda values indicate a strong similarity in the relationship between closely related species. To ensure that the error associated with each annual index value was accounted for within the final model outcome, we used a bootstrap procedure to produce 95% confidence limits around the PGLS regression line. For each bootstrap sample (n = 1000), new values of feeder use and population index for the beginning and end of the time series were drawn at random from the confidence limits around their original estimates, and then used to recalculate estimates of change. The PGLS model was fitted to each bootstrap sample with lambda set at the value estimated for the original model, then 95% confidence limits were calculated from the set of regression coefficients produced.

PGLS, using the pruned consensus tree and maximum likelihood to estimate Pagel’s lambda, was also used to test for differences in population trends between feeder-users and non-feeder users. Here, we used the 1994-2012 habitat-specific trends from Sullivan et al.31 for all breeding birds that are associated with urban areas of Britian and therefore have frequent access to garden bird feeders. These trends were available for 72 species, 33 (46%) of which had been defined as feeder-users. Feeder-users without trends were either winter migrants or did not have suitable data for population estimation. We aggregated the trends for 12 individual habitats to derive two broader trends of interest, urban and non-urban. More specifically, urban trends were estimated using a weighted average of the trends for suburban/urban settlements and rural settlements, accounting for their habitat availability. Non-urban trends were estimated using an weighted average of the trends from all other habitat types (deciduous woodland, mixed woodland, coniferous woodland, upland semi-natural open habitats, lowland semi-natural open habitats, arable farmland, pasture, mixed farming, wetlands and flowing water), accounting for their availability31.

Statistical modelling

To account for temporal auto-correlation, all trend analyses (described above) included an AR(1) correlation structure. The AR(1) correlation structure was found to be optimal for time series modelling across the different response variables, based on the comparison of models with and without different autocorrelation structures (AR1 or AR2) using AIC and the examination of auto-correlation plots for the model residuals. There was no evidence of spatial autocorrelation in bird community indices across gardens within years, according to spline correlograms fitted to the raw data using the ncf R package59. When using GAM(M)s to investigate non-linear temporal trends, year was always fitted in the form of a thin-plate regression spine with a maximum of five degrees of freedom and the gamma parameter was fixed at 1.4 to reduce over-fitting. Generalised least squares (GLS, national-scale data) and linear mixed models (LMM, garden-scale data) were used to determine the significance of linear temporal trends when GAM(M)s fitted to the same data did not indicate non-linearity (e.g. the smoothed trend had one degree of freedom, was not significant, or did not deviate enough from the linear trend to be deemed ecologically meaningful). Significance was determined using maximum likelihood, Wald statistics, χ2 and F-tests as appropriate with alpha set at 0.0560. To identify periods of significant change within non-linear trends, where the rate of change (the slope) was distinguishable from zero given the uncertainty of the model, we estimated the first derivatives of the GAM temporal smooth61,62. A significant change was assumed where the 95% confidence intervals of the first derivatives excluded zero61. All analyses were performed using R version 3.4.363. Trend analyses used the packages mgcv64 and nlme65, and phylogenetic comparative analyses used APE66, phytools67 and caper68.