Coordination of care reduces conflict and predation risk in a cooperatively breeding bird

When two or more individuals cooperate to provision a shared brood, each carer may be able to maximize their payoffs by coordinating provisioning in relation to what others are doing. This investment ‘game’ is not simply a matter of how much to invest, but also of the relative timing of investment. Recent studies propose that temporal coordination of care in the forms of alternation (i.e. turn-taking) and synchrony (i.e. provisioning together) function to mitigate conflict between carers and reduce brood predation risk, respectively. Such coordination is widespread in biparental and cooperatively breeding birds, yet the fitness consequences have rarely been empirically tested. Here, we use a long-term study of long-tailed tits Aegithalos caudatus, a facultative cooperatively breeding bird with active coordination of care, to assess the support for these hypothesized functions for coordination of provisioning visits. First, we found evidence that turn-taking mitigates conflict between carers because, in cooperative groups, provisioning rates and offspring recruitment increased with the level of active alternation exhibited by carers, and with the associated increase in provisioning rate parity between carers. In contrast, offspring recruitment did not increase with alternation in biparental nests, although it was positively correlated with parity of provisioning between carers, which is predicted to result from conflict mitigation. Secondly, synchronous nest visits were associated with a reduced probability of nest predation and thus increased brood survival, especially when provisioning rates were high. We attribute this effect to synchrony reducing carer activity near the nest. We conclude that temporal coordination of provisioning visits in the forms of alternation and synchrony both confer fitness benefits on carers, and despite being intrinsically linked, these different kinds of coordination appear to serve different functions.

Methods
Methods

Study system and general field protocol

Data were collected during the breeding seasons (March–June) of 1994–2022 from an intensively studied wild population of long-tailed tits (Rivelin Valley, UK; c.3km2; 53°23′N, 1°34′W). Nests were located by following adults building nests, typically in low-lying (≤3m) shrubs (c.73% of nests; Higgot 2019) such as Rubus fruticosus, but sometimes in inaccessible tree forks. Nests were monitored at 1–3 day intervals but daily approaching incubation, hatching, and fledging. Incubation starts once the clutch is complete, lasts c.15 days, and all eggs that hatch do so within 24h of the first. After hatching (d0), both parents deliver prey items, such as flies, spiders, and caterpillars. Helpers may join a group at any point between hatching and fledging (d16–18), although helping is more common later in development. Long-tailed tits suffer only minimal chick starvation (0.2% daily per chick; Hatchwell et al. 2004) but lose 71.9% of clutches and broods to predation (Hatchwell et al. 2013), typically by corvids (e.g. Eurasian jay Garrulus glandarius) and mammals (e.g. stoat Mustela erminea). Therefore, most brood mortality can be attributed to complete depredation events, although broods may occasionally be partially depredated (3.9% of successful nests; 9/233).

Clutch size (median = 10; range 4–12; N = 293) in accessible nests was recorded during incubation. Brood size (median = 9; range 1–11; N = 275), and the mass (to 0.1g; mean = 7.4g ± 0.0142 SE; N = 1970) and tarsus length (to 0.1mm; mean = 18.3mm ± 0.0188 SE; N = 1970) of nestlings was recorded on d11. Each chick was ringed under British Trust for Ornithology (BTO) license with a unique combination of color rings. We took 5–20µl of blood by brachial venipuncture (under UK Home Office license) for genetic sex determination of nestlings using the P2–P8 sex-typing primers (Griffiths et al. 1998). Because clutch and/or brood size were important covariates, we limited our analysis to low nests where these metrics were sampled. The biometrics and sex of every chick was known in 77.7% (185/238) of broods. All applicable international, national, and institutional guidelines for the use of animals were followed and all regulated procedures were approved by the Animal Welfare and Ethical Review Body at the University of Sheffield.

Successful local recruitment was recorded when a fledgling attempted to breed in the field site in a subsequent year. Long-tailed tits exhibit female-biased dispersal, so 20–25% of male, but <10% of female fledglings recruit locally (Sharp et al. 2008, 2011). Therefore, we estimated recruitment success from resightings of male fledglings breeding in subsequent years. We used the number of fledged males that were (median = 1; range = 0–5; N = 170 broods), or were not resighted per brood (median = 3; range = 0–7; N = 170) to model the proportion that recruited locally (mean = 0.241 ± 0.0212 SE; N = 170). In our open population, c.40% of breeders were immigrants to our field site per annum; these were captured in mist-nets and ringed. Since the project started, 1531 individuals have been recorded breeding in the site, but many of these never raised a brood and only 576 (37.6%) birds were recorded provisioning broods in watches matching our criteria, which included 239 unique breeding females, 227 breeding males, and 171 helpers (61 of which were also breeders).

Calculating coordination

Provisioning watches (hereafter ‘watches’) were typically performed from d2 every other day until nest failure or fledging. Coordination is not possible when females are brooding young chicks (≤d5), so analysis was restricted to watches when both parents provisioned full-time (median = 3 watches per nest; median brood age = d10; range d6–18; N = 894). Each season c.95% of carers were identifiable by their unique combination of color rings. If an unringed carer provisioned during a watch (4.36% of watches; 39/894) we assumed all feeds were by the same unringed individual and omitted watches including >1 unringed birds.

The watch protocol was consistent throughout the study. Following a c.10 min habituation period, watches usually lasted for one hour between 04:00–19:30, unless curtailed by inclement weather (minimum duration: 30 mins). When a carer provisioned the brood its identity and time was recorded to the nearest minute, either by direct observation through binoculars (15–25m away) or video camera (1–5m away). Watch duration was the time between the first and last feed (mean = 54 mins 1 ± 25 seconds SE; range = 30–117 mins; N = 894). We excluded any watch where identities of feeds or hatch date were unknown, and from nests subjected to experimental manipulation. Long-tailed tits provision their nestlings frequently, with a mean rate of 23.8 feeds/h ± 0.320 SE (range 4.53–69.2; N = 894) per group, in the sample used in this study. The total number of alternated and synchronized feeds was calculated per watch. An alternated feed was any that occurred following a feed by another carer (median = 15 per watch; range 1–68; N = 894), meaning that alternation did not require a consistent pattern of feeds (e.g. A-B-C-A-B-C), just non-consecutive feeds (e.g. A-B-A-C-B-A). A synchronized feed was any non-consecutive feed that occurred within 2 minutes of the previous feed (median = 8 per watch; range 0–55; N = 894). This 2-minute window was chosen to facilitate comparison with prior studies of coordination in this species, which found that rates of synchrony using different window lengths were highly correlated (rp ≥ 0.94; Halliwell et al. 2022) and key results were qualitatively the same (Bebbington & Hatchwell 2016; Halliwell et al. 2022). In watches where one carer performed >50% of feeds some cannot be alternated or synchronized, so we calculated the ‘Maximum possible coordination’ (mean = 87.4% ± 0.414 SE; N = 894) for each watch, which functions as a proxy for provisioning rate parity between carers.

Because some alternation and synchrony are expected by chance (Schlicht et al. 2016, Ihle et al. 2019a, Santema et al. 2019) we used a null model to estimate levels of expected coordination (e.g. Halliwell et al. 2022; Ihle et al. 2019a, b; Johnstone et al. 2014). We randomized the order of each carer’s intervisit intervals within a watch 1000 times, which were then recombined to produce 1000 randomized sequences, and the median number of alternated and synchronized feeds were the expected levels for that watch. Observed and expected alternation and synchrony were used to generate measures of how much each watch deviated from expected, termed ‘Active alternation score’ (mean = 0.0663 ± 0.00499 SE; N = 894) and ‘Active synchrony score’ (mean = 0.206 ± 0.00884 SE; N = 894), respectively. Alternation scores were calculated from log(observed alternated feeds + 0.5) – log(expected alternated feeds + 0.5), and synchrony scores likewise. Therefore, a positive score means that carers coordinated more than expected by random chance (i.e. if they provisioned independently), while a negative score denotes less coordination than expected, which may occur if carers provision in a manner that actively avoids alternation (i.e. in bouts of successive uninterrupted feeds by the same carer). We added 0.5 to each value to avoid taking the log of zero.

Statistical analysis

Statistical analysis was performed on R version 4.2.3 (R Core Team 2023). Models were built using lme4 (Bates et al. 2015), coxme (Therneau 2022), and analyzed with lmerTest (Kuznetsova et al. 2017). Figures were produced using ggplot2 (Wickham 2016), survminer (Kassambara et al. 2021) and cowplot (Wilke 2020). When investigating the influence of alternation and synchrony on reproductive success we used ‘Active alternation score’ and ‘Active synchrony score’, respectively. These were analyzed in separate models because they are intrinsically correlated as synchronized feeds are, by definition, alternated. Additionally, because previous studies found that active alternation decreased with helper presence, and active synchrony with provisioning rate (Halliwell et al. 2022), we included these terms and their interactions with each coordination score as explanatory terms in our alternation and synchrony analyses, respectively. A significant interaction term indicates that the importance of coordination varied with the presence of helpers (‘Helped during watch?’ – binary factor denoting whether a pair was assisted by helpers in each watch, and ‘Nest helped’ – a binary factor denoting whether a pair was helped within the range of watches analyzed) or ‘(aggregate) provisioning rate’ – a continuous numerical variable denoting the number of feeds performed per hour by all carers per watch (or across several watches).

Provisioning rate (Prediction 1a):

We used a normally distributed linear mixed effects model (LMM) to investigate the relationship between alternation and provisioning rate in the Full sample of watches (N = 871 at 275 nests; Figure S1). The response variable was the log-transformed provisioning rate per watch and the explanatory terms of interest were ‘Active alternation score’ and its interaction with ‘Helped during watch?’. Covariates and random effects used here and throughout are described below.

Chick mass (Predictions 1b):

We fitted normally distributed LMMs to investigate the relationship between alternation and chick mass using a subset of watches taken prior to biometric assessment from nests where all chicks’ biometrics and sexes were known (Chick mass sample sex known; N = 360 at 185 nests containing 1533 chicks). The response variable was each chick’s mass on d11. Terms of interest were ‘Active alternation score’ (aggregated across each appropriate watch) and its interactions with ‘Nest helped?’ and ‘Aggregate provisioning rate’.

Recruitment rate (Predictions 1c):

We fitted a binomially distributed generalized LMM to investigate the relationship between alternation and recruitment rate, using a subset of watches from successful nests where each chick’s sex was known (Recruit sample; N = 574 watches at 170 nests containing 719 male chicks). The response variable was a two-column variable (number of males recruited, number of males not recruited), which functions as a measure of proportion recruited. Terms of interest were ‘Aggregate active alternation score’ and its interaction with ‘Nest helped?’. Because this interaction term was significant, we also repeated this analysis on sub-samples of biparental (2 carers) watches (N = 331 at 101 nests containing 420 male chicks) and cooperative (>2 carers) watches (N = 243 at 69 nests containing 299 male chicks).

Predation and survival (Prediction 2):

To investigate survival time, we used a Cox proportional hazard mixed model (Therneau 2022) with the two-column response variable (days until event, fledged or failed) applied to the Predation sample, which included watches from nests depredated prior to ringing, for which brood size was assumed equal to clutch size (N = 894 watches at 293 nests). Data were right censored with all fledged broods defined as age 18. Terms of interest were ‘Aggregate active synchrony score’ and its interaction with ‘Aggregate provisioning rate’.

Covariates and random effects:

We controlled for biologically important covariates that could influence provisioning behavior. ‘Provisioning rate variation’ – a continuous numerical variable denoting variation in the provisioning rate during each watch (mean = 0.582 ± 0.00733 SE; N = 894); included because as carers’ intervisit intervals become more consistent, the null model’s ability to disrupt patterns of coordination fundamentally diminishes (supplementary material Tables S1 and S2; Figure S2). ‘Brood size’ (linear and quadratic) – an integer numerical variable denoting the number of live chicks on d11; included because it affects demand on carers. ‘Watch duration’ – a continuous numerical variable denoting the time (minutes) between the first and last recorded feeds per watch; included because the total number of feeds increases with watch duration. ‘Watch start time’ – a continuous numerical variable denoting the time each watch started (mean = 9:54am ± 5 mins SE; N = 894); included because long-tailed tits have higher provisioning activity soon after sunrise (Hatchwell et al. 2004; MacColl & Hatchwell 2002). ‘Maximum possible coordination’ – a continuous numerical variable denoting the highest theoretical percentage of feeds that could be alternated (or synchronized) during a given watch. ‘Hatch date’ (linear and quadratic) – an integer numerical variable denoting the number of days between March 1 and hatching (Median = May 3; range April 15–June 6; N = 293); included to account for within-season environmental variation. ‘Brood age’ – an integer numerical variable denoting the number of days between hatching and a watch; included because provisioning rate increases with age. ‘Chick sex’ – a binary factor denoting whether a chick was male or female; included because male chicks are typically heavier (Nam et al. 2011). ‘Brood sex ratio’ – a continuous numerical variable denoting the ratio of female:male chicks within a given brood (mean = 0.511 ± 0.0124 SE; N = 185); included to account for differences in chick mass between broods with different sex ratios. ‘Tarsus length’ (linear and quadratic) – a continuous numerical variable denoting chick tarsal length (mm). ‘Mean carer number’ – a continuous numerical variable denoting the mean number of carers observed provisioning during all watches of a given nest (mean = 3.21 ± 0.0798 SE; N = 69, cooperative nests only). Random effects were as follows. ‘Year’ – factor denoting the year a watch was performed. ‘Nest ID’ – factor denoting the identity of a nest. ‘Pair ID’ – factor denoting the unique combination of parents. ‘Female ID’ and ‘Male ID’ – factors denoting the unique identity of each mother and father, respectively. ‘Rowref’ – observation level random effect used to account for overdispersion in Poisson-distributed models. Full model details are available in the supplementary material (Table S3).