In species with biparental and cooperative brood care, multiple carers cooperate by contributing costly investment to raise a shared brood. However, shared benefits and individual costs also give rise to conflict among carers over investment. Coordination of provisioning visits has been hypothesized to facilitate the resolution of this conflict, preventing exploitation, and ensuring collective investment in the shared brood. We used a 26-year study of long-tailed tits, Aegithalos caudatus, a facultative cooperative breeder, to investigate whether care by parents and helpers is coordinated, whether there are consistent differences in coordination between individuals and reproductive roles, and whether coordination varies with helper relatedness to breeders. Coordination takes the form of turn-taking (alternation) or feeding within a short time interval of another carer (synchrony), and both behaviors were observed to occur more than expected by chance, i.e. ‘active’ coordination. First, we found that active alternation decreased with group size while active synchrony occurred at all group sizes. Secondly, we show that alternation was repeatable between observations at the same nest, while synchrony was repeatable between observations of the same individual. Active synchrony varied with reproductive status, with helpers synchronizing visits more than breeders, although active alternation did not vary with reproductive status. Finally, we found no significant effect of relatedness on either alternation or synchrony exhibited by helpers. In conclusion, we demonstrate active coordination of provisioning by carers and conclude that coordination is a socially plastic behavior depending on reproductive status and the number of carers raising the brood.

Study system and data collection

We used data from a long-term study of a population of long-tailed tits in the Rivelin Valley, Sheffield, UK (53°23′N, 1°34′W) from 1994 to 2019. The field site is ~3km² with a population of 25-72 breeding pairs (Hatchwell 2016). Each year ~95% of adult birds were marked (under British Trust for Ornithology license) with a unique combination of two color rings on one leg and a BTO metal ring on the other. The adult annual mortality rate is ~50% (Meade and Hatchwell 2010), and ~20% of new recruits into the adult population were ringed as nestlings in the study site, while the remaining ~80% of new recruits were unringed adult immigrants that dispersed into the population. Unringed birds were captured in mist-nests during the nest-building period and DNA samples were collected (under Home Office license) for genotyping and social pedigree reconstruction. Nests were found by following adults and once located, were monitored every 2-3 days, with daily visits around the expected hatch date. Median clutch size is 10 eggs (range: 4-12), which are incubated for ~15 days (Hatchwell 2016). Hatching is extremely synchronous within clutches, with all chicks typically hatching within 24 hours of the first. Initial hatch date was recorded as day 0, and chicks were ringed and counted on day 11. Protocols for provisioning watches (hereafter ‘watches’) were broadly consistent throughout the study. In most cases, watches of duration ~60 minutes were carried out every other day, starting on day 2, either by direct field observation or by video camera, for later review (69% of watches were between 45 and 65 minutes). Watches were carried out between 04:00 and 18:00, with 89% starting between 06:00 and 14:00. Watches were performed until a nest was predated, abandoned or chicks fledged, typically on day 16-18.

For ~5 days post-hatching nestlings are brooded regularly by their mothers, who provision offspring only occasionally, while fathers either feed the offspring directly or give food to the mother, who then feeds the chicks. We restricted our analysis, therefore, to watches at day 6 and older, when both parents provision offspring directly. Long-tailed tits exhibit facultative cooperative breeding (Lack and Lack 1958, Hatchwell 2016), meaning nests may be uniparental (1 carer, in the rare event of a parent dying), biparental (2 carers) or cooperative (>2 carers). For this study we restricted analysis to watches of biparental and cooperative nests with up to 5 carers (i.e. social parents and up to 3 helpers). Our dataset contained 65% (516) of watches from biparental nests and 21% (171), 11% (88) and 3% (20) from nests with 3, 4 and 5 carers, respectively. Before starting a watch, ~10 minutes was usually allowed for birds to recover from observer disturbance and we restricted analysis to watches of total duration ≥ 30.0 minutes and ≤ 180.0 minutes, with duration defined as the time between first and last observed feeds. Mean watch duration (± SD) was 54.8 ± 14.4 minutes (range 30-118 minutes, N = 795 watches). We omitted watches where the identity of any provisioning visit was unknown, and from nests that were manipulated for other behavioral studies (e.g. Meade et al. 2011). Watches were used from 24 years between 1994 and 2019, with 2007 and 2009 excluded because experiments conducted in those years meant that they contained no watches matching our criteria. In total, our dataset included 795 watches performed at 250 unique nests, involving 192 different breeding males, 203 breeding females and 144 helpers.

Calculating coordination

We analyzed alternation and synchrony as the absolute number of alternated and synchronized feeding visits in a provisioning watch, respectively. We defined an alternated visit as any non-consecutive provisioning visit (i.e. a visit occurring after the provisioning visit of any carer other than itself) and a synchronized visit as an alternated visit occurring within 2-minutes of the previous feed (Figure 1). We chose an interval of 2-minutes in accordance with previous studies (Mariette and Griffith 2015, Bebbington and Hatchwell 2016, Ihle et al. 2019a), and further analyses revealed that number of synchronized visits was highly correlated for 1, 2 and 3-minute intervals (Pearson correlations: 1 v. 2 min, r = 0.97, df = 793, P < 0.001; 2 v. 3 min, r = 0.97, df = 793, P < 0.001; 1 v. 3 min, r = 0.94, df = 793, P < 0.001), and analyses of synchrony with different intervals produced qualitatively the same results.

We calculated observed alternation and synchrony directly from visit sequences and times recorded through field observation, generating coordination measures per watch and for each individual carer present in each watch (Figure 1). We generated expected data by null model randomization of observed data, with the binary factor ‘Data type’ specifying whether data were observed or expected. In accordance with the most conservative method of calculating expected alternation and synchrony recommended by Ihle et al. (2019b), our null models used a within-watch, within-individual randomization procedure in which the order of provisioning visits within a watch was randomized in a manner that preserved the length and identity of each period between feeding visits (inter-visit intervals) (Figure S1; supplementary material). We calculated expected numbers of alternated and synchronized visits, both for group total and for individual carers, from the median of 1000 iterations of the null model applied to each provisioning watch. We used median values to preserve integer values for subsequent analysis in Poisson-distributed linear models; mean and median values were highly positively correlated (Pearson correlations: alternated visits, r = 0.99, df = 793, P < 0.001; synchronized visits, r = 0.99, df = 793, P < 0.001).

Calculating kinship

To calculate pairwise values of pedigree relatedness of helpers to parents we constructed an additive relationship matrix using the R package NADIV (Wolak 2012), partially reconstructed using molecular genetic data from up to 17 microsatellite loci to perform offspring-parent reconstruction on CERVUS v. 3.0.7 (Kalinowski et al. 2007) and sibling-sibling reconstruction on KINGROUP v.2 (Konovalov et al. 2004). Building on the social pedigree and protocol used in Leedale et al. (2018, 2020) we expanded the pedigree to include 2018 and 2019 data. Our study population is open, so even after reconstruction the social pedigree remained incomplete; therefore, where necessary we omitted data with incomplete pairwise relatedness metrics to either social parent.

Statistical analysis

All statistical analysis was performed in R version 4.0.2 (R core development team, 2020). All models were built and analyzed using the lme4 package (Bates et al. 2015) and lmerTest (Kuznetsova et al. 2017), except for our repeatability models which were built and analyzed using the rptR package (Stoffel et al. 2017).

Collective coordination models (Alt-C and Sync-C)

To investigate collective alternation and synchrony performed by all carers at a nest we defined two Poisson-distributed generalized linear mixed effects models (GLMM) named ‘Alt-C’ and ‘Sync-C’, respectively. The response variables to these models were the number of alternated visits (collective) and synchronized visits (collective) by all carers at each watch, respectively. To control for observation and population structure, these models were built with the following random effects: ‘Year’, ‘Nest ID’, ‘Watch ID’, ‘Male ID’, ‘Female ID’, ‘Helper1 ID’, ‘Helper2 ID’, ‘Helper3 ID’ and ‘Row reference’ (see Table 1 for explanation). The fixed effects tested were as follows: ‘Data type’ (observed vs. expected values of alternation and synchrony), ‘Provisioning rate (collective)’, ‘Carer number’, ‘Watch duration’, ‘Brood size’, ‘Time of day’, ‘Brood age’, ‘Hatch date’ and ‘AMax (or SMax)’ (Table 1). We focused our analysis on ‘Data type’ and 2-way interactions with other fixed effect terms, as a disparity between observed and expected data represents the level of active coordination performed.

Individual coordination models (Alt-I and Sync-I)

To investigate the effect of carer status on alternation and synchrony performed by a given carer, we built two Poisson-distributed GLMMs named ‘Alt-I’ and ‘Sync-I’, respectively. The response variables to these models were the number of alternated visits (individual) and synchronized visits (individual), respectively. These models were built with the following random effects: ‘Year’, ‘Nest ID’, ‘Watch ID’, ‘Carer ID’ and ‘Row reference’ (Table 1). The fixed effects tested were as follows: ‘Data type’, ‘Carer status’, ‘Provisioning rate (individual)’, ‘Carer number’, ‘Watch duration’, ‘Brood size’, ‘Time of day’, ‘Brood age’, ‘Hatch date’ and ‘Amax’ or ‘SMax’ (Table 1). In this analysis, the focus was on the interaction of ‘Data type’ with ‘Carer status’ because this term represents the disparity in active coordination between carers of different breeding status.

Repeatability models (Alt-R and Sync-R)

To investigate the repeatability of active alternation and synchrony within nests and within individuals we constructed two Gaussian-distributed GLMMs named ‘Alt-R’ and ‘Sync-R’, respectively. In these models, response variables were the number of actively alternated (individual) and actively synchronized (individual) visits by an individual during a watch, respectively (active alternation range: -3 to 6; active synchrony range: -7 to 9). We used these metrics because repeatability analyses required active coordination to be the response variable, rather than using interaction terms with ‘Data type’ as in our other models. To control for the effect of confounding factors on active coordination we included all fixed effects previously found to significantly influence either individual alternation or synchrony (Alt-I, Sync-I) and, using the rptR function, ran models with 1000 bootstrapped simulations and 1000 permutations. We investigated both within-nest repeatability (‘Nest ID’) and within-individual repeatability (‘Carer ID’) in the same models. Additionally, we included ‘Year’ as a random effect to account for between-year variation. As active coordination was the response variable and a Gaussian error distribution was used, ‘Watch ID’ and the ‘Row reference’ random effects were not required for these models. We present our repeatability results as values of R and extracted 2.5% and 97.5% confidence intervals (CIs) in addition to P-values.

In our dataset many individuals were observed provisioning at only one nest, potentially confounding repeatability of an individual’s behavior with the potential effect of common nest factors. Therefore, we ran the repeatability analysis on a subset of data, restricted to carers observed provisioning at two or more nests (Table S2, supplementary material). Results from these models were qualitatively the same as those for the full dataset for both within-nest and within-individual repeatability for both alternation and synchrony models.

Kinship models (Alt-K and Sync-K)

To investigate the effect of kinship to the breeding pair on alternation and synchrony performed by helpers we constructed two Poisson-distributed GLMMs named ‘Alt-K’ and ‘Sync-K’, respectively. Just as with ‘Alt-I’ and ‘Sync-I’, the response variables to these models were the number of alternated visits (individual) and synchronized visits (individual) performed by an individual during a watch, respectively, however analysis was restricted to helpers whose pedigree kinship with breeders was known. These models were built with the same random and fixed effects as ‘Alt-I’ and ‘Sync-I’ but with the addition of three fixed effects: ‘Sex’, ‘Kinship with father’ and ‘Kinship with mother’ (Table 1). We focused our analysis on the interactions of ‘Data type’ with our kinship terms as these represent the relationship between the level of active coordination and relatedness.

Raw provisioning watch data are available as it was recorded directly from field observation prior to processing.

Collective (and individual) level processed data are available. These are the values per watch (and per individual) for terms of interest e.g. number of alternated visits for observed (from field observation) and expected (from null model randomization). These also contain nest-specific information e.g. brood size, hatch date.

Relatedness data are available in the form of each individual's known dam and sire, either from known brood associations or via reconstruction from genetic loci analysis.