Our study focused on two wetlands used by wintering waterbirds in the UK; the Wildfowl & Wetland Trust (WWT) Centre reserves at Slimbridge (51°44ʹ29.3ʺ N, 2°24ʹ21.52ʺ W) and Caerlaverock (54°59ʹ2.4ʺ N, 3°30ʹ0ʺ W), in southwest England and southwest Scotland, respectively. Both wetland sites feature a mosaic of aquatic and terrestrial habitats, most notably small lakes used by waterbirds for feeding and roosting. In recent winters Slimbridge has supported both Bewick’s Swans and Mute Swans, whilst Caerlaverock has supported Mute and Whooper Swans (Black and Rees 1984). The mean of the peak winter counts of individual swans between 2014/2015 and 2018/2019 was 149 Bewick’s Swans (range 116–212) and 441 Mute Swans (range 374–500) in the Severn Estuary (which includes Slimbridge), and 337 Whooper (range 257–487) and 74 Mute Swans (range 70–83) in the Solway Estuary (which includes Caerlaverock) (Frost et al. 2020). Swans at both sites share feeding and roosting habitat with a range of smaller-bodied waterbird species, including geese such as Canada Geese (Branta canadensis) and Greylag Geese (Anser anser), dabbling ducks such as Northern Mallard (Anas platyrhynchos), Eurasian Teal (Anas crecca), and Northern Pintail (Anas acuta), diving ducks such as Tufted Duck (Aythya fuligula) and Common Pochard (Aythya ferina), Rallidae such as Eurasian Coot (Fulica atra) and Common Moorhen (Gallinula chloropus), as well as Common Shelduck (Tadorna tadorna), and gull species (Larus spp.). The non-migratory Mute Swans are resident at both Slimbridge and Caerlaverock throughout the year, and Bewick’s Swans are typically present at Slimbridge between November and February, whilst Whooper Swans use Caerlaverock between October and March (Black and Rees 1984; Rees 2006).
We used two ways of collecting data concurrently from observations made during periods of one hour in duration: one way was used to record aggression between any two individuals of any waterbird species to address our conspecific hypothesis, and the second way was used to record swan-specific aggression to address our density and winter decline hypotheses. First we used all occurrence sampling to record all incidents of aggression between individuals over the course of the hour-long observation period (Altmann 1974), based on an ethogram of the aggressive behaviours observed during preliminary observations at both sites in October 2018. These preliminary observations indicated a number of aggressive behaviours consistent with previous work (Johnsgard 1965), including strikes made with the bill, wings, or body, as well as chasing and lunging at another individual. For each aggressive behaviour that was observed, the species identities of both the aggressor and its opponent were recorded; all species were identified on the basis of size, body shape, and plumage characteristics, with the aid of an online photo-identification guide (RSPB 2018).
To address our density and winter decline hypotheses regarding the variations in aggression with swan numbers and between months, we used focal sampling (Altmann 1974) to quantify the total time that Bewick’s and Whooper Swans spent engaged in aggressive interactions with other swans. During each hour-long observation period, an observer selected a swan at random and used a stopwatch to record the duration of each aggressive interaction with another swan in a 10-min observation period; hence, six individual swans could be observed during each hour-long observation. A focal observation duration of 10 min was selected in order to make our study comparable with earlier time-activity budget of swans that used an observation duration of 10 min (e.g. O’Hare et al. 2007; Tatu et al. 2007; Wood et al. 2019b). Both immediately before and after each hour-long observation period, the number of swans that could be observed was counted, with the mean average taken as the number of swans present during that observation.
All behavioural data were collected remotely via live-streaming webcams (AXIS Q6035-E PTZ Dome Network Camera), which were fixed in place at both sites. Both webcams faced directly outwards over the study lake from the shore, and each webcam maintained the same zoom so that the field-of-view was standardised across all observation periods (Additional file 1: Figure S1). As our study was conducted entirely via the webcams, we do not know the precise numbers of birds outside of the cameras field-of-view, but we suspect it to be low as the cameras covered major proportions of the surface area of each lake (Additional file 1: Figure S1). Webcams have been shown previously to be useful tools in behavioural studies, which allow remote collection of data with limited disturbance to the focal birds (e.g. Anderson et al. 2011; Schulwitz et al. 2018). During winter months in 2018/2019 (November 2018–February 2019, inclusive) and 2019/2020 (November 2019–December 2019, inclusive), webcam footage from WWT Slimbridge was watched for 1 h on an average of 7.5 days per month at either 08:30 a.m. to 09:30 a.m. (hereafter “AM”), 11:30 a.m. to 12:30 p.m. (hereafter “MID”) or 14:30 p.m. to 15:30 p.m. (hereafter “PM”); these times were selected to achieve a balanced study design and to avoid coinciding with the periods when food (wheat grains) is provided as part of a public engagement programme. Winter storms at the start of January 2020 regrettably damaged the webcam and prevented further data collection. We aimed to use the same methodology for WWT Caerlaverock, however, the webcam suffered a technical failure after data for only November 2018 could be collected. Similar technical issues precluded data collection at WWT Caerlaverock during winter 2019/2020. Our sampling methodology of ten observations per hour, on an average of 7.5 days per month, over 6 months, yielded a total of 450 observations at Slimbridge; however, swans were only present during 282 of these. In contrast, swans were present at Caerlaverock during every observation period. In total, we therefore obtained 282 focal observations (of 10 min each) from WWT Slimbridge and 42 focal observations (of 10 min each) from WWT Caerlaverock. While we cannot discount the possibility that some individuals were observed on multiple occasions, the aforementioned large numbers of individual swans present at both sites mean that this is unlikely to have been a major issue.
All statistical analyses were carried out using R version 3.6.3 (R Core Team 2020). To address our conspecific hypothesis that intraspecific interactions would be observed more frequently than interspecific interactions, we used two-tailed binomial tests to assess the statistical significance of the deviation of the proportion of intraspecific aggressive interactions of each species from 0.5 (which would indicate equal numbers of intraspecific and interspecific interactions). For each of our focal species, separate tests were carried out for aggressive interactions (i) targeted at other individuals, and (ii) received from other individuals. Statistically significant differences between proportions were attributed where P < 0.05, after P values had been adjusted using Holm-Bonferroni corrections to account for multiple comparisons (Holm 1979). Our binomial tests also allowed 95% confidence intervals to be estimated for the proportions of intraspecific and interspecific interactions (Clopper and Pearson 1934). In addition, we used two-sample binomial tests for equality of proportions to assess the significance of the differences between species in their proportion of interspecific aggressive interactions recorded towards other birds; P values were again adjusted using Holm-Bonferroni corrections (Holm 1979).
To address our density and winter decline hypotheses regarding the variations in aggression with swan numbers and between months, we used zero-inflated generalized linear mixed effects models (ZIGLMMs), using the glmmTMB R package (Brooks et al. 2017). Zero-inflated models were required because of the relatively high proportions of zeros in the data sets (i.e. observations during which no aggression was recorded). In each model, our response variable was the number of seconds spent in aggressive interactions by each individual recorded during our focal observation. As the aforementioned issues with the webcam at WWT Caerlaverock prevented data collection after November 2018, the resulting datasets for Bewick’s and Whooper Swans were markedly different in terms of the sample size and the temporal replication. Therefore, we analysed the data for Whooper and Bewick’s Swans separately. For each species, we ran and compared candidate models that comprised all possible combinations of additive and two-way interactions between our explanatory variables, as well as the null model. For Whooper Swans, we considered the following explanatory variables: (i) the time of day of the observation, a categorical factor with three levels (AM, MID, and PM), and (ii) the mean number of swans present during the observation. For Bewick’s Swans we considered the following explanatory variables: (i) the time of day of the observation, a categorical factor with three levels (AM, MID, and PM), (ii) the mean number of swans present during the observation, (iii) the month of observation, a categorical factor with four levels (November, December, January, and February), and (iv) the winter of observation, a categorical factor with two levels (“A” = 2018/2019 and “B” = 2019/2020). In addition, a categorical variable unique to each observation block (termed the ‘observation identity’), was fitted as a random intercept in each model to account for the non-independence of individual swans observed within the same hour-long observation period. Preliminary comparisons of global models using second-order Akaike’s Information Criteria (AICc; Burnham et al. 2011) showed that for the zero-inflated generalized linear mixed effects models of the Whooper Swan data, a negative binomial distribution in which the variance increased linearly with the mean (Brooks et al. 2017) performed better than either a negative binomial distribution in which the variance increased quadratically with the mean (ΔAICc = 4.31), or a Poisson distribution (ΔAICc = 169.40). For Bewick’s Swans there was little difference between models with either of the two negative binomial distributions (ΔAICc = 0.24), whilst the Poisson distribution did not converge. Therefore, in all subsequent models we used the negative binomial distribution in which the variance increased linearly with the mean.
To ensure that collinearity did not confound our modelling (Dormann et al. 2013), we tested for covariance among our explanatory variables. To our knowledge, Variance Inflation Factors (VIFs), which are typically used to identify collinear variables in linear models (Dormann et al. 2013), are not currently available for ZIGLMMs, and so alternative methods were used. One-way Analysis of Variance (ANOVA) was used to test for covariance between our continuous variable (number of swans present) and each of our categorical variables: time of day, month, and winter. Significant covariance was inferred where statistically significant differences in the mean number of swans present per categorical variable level were detected. Values for the number of swans present were square-root transformed so that model residuals satisfied the assumptions of the ANOVA tests (Zuur et al. 2010). Associations between the frequencies of pairs of categorical variables were tested using χ2 tests where all frequencies were ≥ 5, or Fisher’s exact test where one or more frequencies were < 5 (Crawley 2013). Significant covariance between variables was inferred where statistically significant associations between the variables were found. Using these methods, we found for Whooper Swans a significant effect of time of day on the number of swans present (ANOVA: F2,39 = 47.27, P < 0.001). For Bewick’s Swans we detected covariance between the number of swans present and both the time of day (ANOVA: F2,279 = 12.24, P < 0.001) and winter (ANOVA: F1,280 = 22.81, P < 0.001), as well as between winter and month (Fisher’s exact test: P < 0.001); all other tests were non-significant (P > 0.05). Consequently, these collinear variables were not permitted within the same candidate models. Therefore, in total we ran 3 and 11 candidate models for the Whooper Swan and Bewick’s Swan datasets, respectively, accounting for all possible non-collinear combinations of additive and two-way interactions between our explanatory variables.
We compared the relative support of each candidate model using AICc, calculated using the MuMIn R package (Barton 2019). Typically the model with the lowest AICc value is considered to be the best-supported by the data, but we also considered models to be competitive where AICc values < 6.0, following the advice of Richards (2008) for dealing with overdispersion. Furthermore, to avoid selecting models with uninformative parameters, we considered that a model with one additional parameter was competitive only if the associated AICc value was lower than the more parsimonious model (Arnold 2010). Three further metrics for each model were used as indicators of the relative strength of support in the data, to facilitate more detailed comparisons among our candidate models: (i) the probability of a model being the best-fitting model compared with the best-supported model shown by AICc (Relative Likelihood RL), (ii) the ratio of ΔAICc values for each model relative to the whole set of candidate models (Akaike weight wi), and (iii) how many more times less likely a model is to be the best-fitting model compared with the best-supported model shown by AICc (Evidence Ratio ER) (Burnham et al. 2011). In addition, to quantify the explanatory power of each model (Mac Nally et al. 2018), the conditional and marginal R2 values, which represented the proportion of the between-swan variance in the time spent on aggression that was accounted for both the fixed and random effects combined and the fixed effects alone, respectively (Nakagawa et al. 2017), were calculated for each model using the sjstats R package (Lüdecke 2020). Finally, Tukey’s post hoc comparisons of the estimated marginal means of variables within our best-supported models were carried out using the emmeans R package (Lenth 2020).