Turn-Taking Between Provisioning Parents: Partitioning Alternation

Introduction

In species with bi-parental care, there is an evolutionary conflict of interest between the two parents (Trivers, 1972). This occurs because parents caring for common offspring share the benefit of their joint investment but only pay the cost of their own care, with the consequence that each parent is selected to exploit its mate by providing a smaller share of the care (Trivers, 1972; Lessells, 2012). A number of theoretical models have been developed to investigate how sexual conflict affects the evolutionarily stable amount of care that parents devote to offspring. Despite the diversity of patterns of parental investment that these models represent, including both “sealed bids” (Houston and Davies, 1985) and “negotiation” (McNamara et al., 1999; Johnstone and Hinde, 2006; Lessells and McNamara, 2012), they all predict that the evolutionarily stable outcome of sexual conflict is a decrease in parental care and reduction in parent and offspring fitness compared with completely cooperating parents. More recently, however, Johnstone et al. (2014) have developed a model in which the ESS involves a form of conditional cooperation between the provisioning parents. The negotiation mechanism is specified by the rates at which each of the parents makes provisioning visits to the brood, depending on whether it is, or is not, the last to visit the nest. The evolutionary stable outcome of this strategy is a “turn-taking” rule in which each parent does not provision when it is the last to feed, but only after a visit by the mate, leading to strict alternation of the nest visits. The ESS is remarkable in that it results in completely cooperative behavior, with each parent having maximum fitness.

Johnstone et al.'s (2014) model has triggered studies on several avian species which have shown that alternation of nest visits does indeed occur at above the rate expected if the two parents provision independently at a constant rate (Johnstone et al., 2014; Bebbington and Hatchwell, 2016; Koenig and Walters, 2016; Iserbyt et al., 2017; Savage et al., 2017; Baldan et al., 2019). However, although this alternation of nest visits is consistent with a turn-taking strategy by the parents, other processes may give rise to alternation. First, alternation can arise simply because parents cannot make successive visits to the nest within a short interval. This latency between consecutive visits by the same parent, referred to as a “refractory period” by Johnstone et al. (2014), can per se produce some degree of alternation in nest visits because it increases the likelihood that the next visit is made by the other parent. Although Johnstone et al. (2014) focused on the presence of a refractory period, any reduction in the variance in the length of inter-visit intervals (IVIs, the time intervals between two consecutive provisioning events by the same parent) will increase the amount of alternation. In order to appraise the effects of a refractory period on the amount of alternation, Johnstone et al. (2014) suggested a randomization procedure that shuffles the sequence of IVIs while maintaining the original frequency distribution of IVIs. Empirical studies using this randomization technique have revealed that the IVI frequency distribution (including the refractory period) accounts for only a small amount of the observed alternation (Johnstone et al., 2014; Bebbington and Hatchwell, 2016; Savage et al., 2017).

A second process other than turn-taking that can produce alternation of nest visits is temporal environmental variation that affects the provisioning rates of both parents, so that long IVIs of one parent tend to be matched with long IVIs of the other (Schlicht et al., 2016). We refer here to this correlation in the provisioning rates of the two parents as provisioning CTH (correlated temporal heterogeneity). Schlicht et al. (2016) used an example in which the provisioning rates of both parents decreased over the observation period to demonstrate that alternation could be generated in this way. However, changes in provisioning rate can be sudden or gradual, occur multiple times, be of different durations and occur in either direction (Ihle et al., 2019; Santema et al., 2019) and all of these kinds of changes can potentially produce alternation—what is critical is whether these changes are correlated within pairs (Ihle et al., 2019). Temporal environmental variation (in the broadest sense) capable of generating provisioning CTH encompasses a range of behavioral and ecological factors, including offspring begging (Ottosson et al., 1997; Hinde, 2006; Hinde and Kilner, 2007), predation risk (Fontaine and Martin, 2006; Ghalambor et al., 2013), food availability (Naef-Daenzer and Keller, 1999; Tremblay et al., 2005), weather conditions (Radford et al., 2001; Wiley and Ridley, 2016), and simple diurnal variation. The diversity of possible factors means that it is difficult or impossible to be sure that all relevant sources of temporal environmental variation have been identified, and hence that it may be impossible to exclude these as an explanation for alternation by provisioning parents (Ihle et al., 2019; Santema et al., 2019).

The aim of this paper is to add to the previous studies of alternation by provisioning parents which have attempted to understand the kinds of conclusions that can be reached from empirical measurements of alternation rates. First, we use a simulation model to investigate whether turn-taking per se can produce statistical patterns (a refractory period and provisioning CTH) that have previously been implicitly assumed to be the product of other ecological and behavioral processes (foraging constraints and temporal environmental variation, respectively). Second, we incorporate the insights gained from this into a conceptual framework, and combine this with Johnstone et al.'s (2014) randomization analysis and a novel analytical approach in which “pseudo-pairs” are created by analytically pairing the provisioning sequence of a parent at one nest with the contemporaneous provisioning sequence of a parent of the opposite sex at a nearby nest, allowing us to partition the observed amount of alternation into different components. Third, we apply these techniques to data from parental great tits Parus major provisioning their broods in the wild. To facilitate these analyses we used a recently developed alternation score (Baldan et al., 2019) and a novel measure of provisioning CTH.

Measures of Alternation and Provisioning CTH

Alternation Score

We used an alternation score which measures the deviation of the observed amount of alternation from that expected given the proportion of visits by the two parents (Baldan et al., 2019):

All the visits within a sample period are given an alternation status based on the sex of the parent making the previous visit (same sex = “non-alternated”; different sex = “alternated”). The expected numbers of alternated and non-alternated visits are calculated from a 2 × 2 contingency table with the sex of the parent at the focal visit cross-tabulated against the sex of the parent at the previous visit.

An alternation score of zero represents the amount of alternation expected by chance, a value of <0 means that the observed alternation of the visits is lower than expected by chance, and a value of >0 means that the observed alternation of the visits is greater than expected by chance. Because the score is based on log odds, alternation score is expected to be additive on a linear scale, allowing it to be partitioned between different factors and processes (Baldan et al., 2019). We did not attempt to calculate an alternation score separately for the two parents because the 2 × 2 contingency table on which the score is based has only one degree of freedom.

Calculating Provisioning CTH

We evaluated the amount of provisioning CTH in the provisioning sequence of a pair by calculating the Pearson correlation coefficient between the provisioning rates of the male and female parent sampled at uniformly spaced points in time within an observation period. We used the inverse of the length of the IVI of each parent in which each sampling point fell as the provisioning rate of the parent at that point in time. In the analyses below, we selected 20 points which, because the observation periods were 4 h long, were spaced at 11.43 min intervals. When the sampling interval is less than the longest IVI it is possible for successive sampling points to fall within the same IVI of one or both parents. These values are included to avoid biasing the correlation coefficient, and their inclusion does not inflate sample size in analyses because the variable that is analyzed is the Pearson correlation coefficient, not the individual pairs of values. We chose 20 sampling points because this gave a reasonable sample size on which to calculate provisioning CTH values, while limiting the number of IVIs that were sampled more than once. A brief analysis using other sampling intervals produced broadly similar results, but confirmed that using longer sampling intervals (and hence fewer sampling points) produces estimates that lack sufficient precision, while shorter sampling intervals may create problems related to sampling close to the beginning and end of the sampling period (cf Baldan et al., 2019, Appendix 2). In calculating provisioning CTH on our data for randomized pseudo-pairs (see Provisioning CTH), we obtained a value that was significantly less than zero from randomized data (mean CTH = −0.003 ± 0.001 (SE), t = −2.566, p = 0.021), whereas the expectation for randomized data is zero. We do not know how this bias arose, but one possibility is that it was generated by sampling effects related to the IVIs that are “cut” by the start and end of the sampling period (cf Baldan et al., 2019, Appendix 2). To put the tiny bias in context, the mean of −0.003 is based on the values from the individual randomizations which range from −0.80 to 0.93—a range more than 500 times larger than the value of the overall bias. In addition, provisioning CTH is a correlation coefficient. When r = 0.003, the resultant r² value implies that <0.001% of the heterogeneity in provisioning rate of one member of the pair is explained by heterogeneity in provisioning rate of the other. For these reasons we have ignored the bias in our analyses below.

Provisioning CTH can potentially occur at a range of time scales (Ihle et al., 2019; Santema et al., 2019). Our measure assesses provisioning CTH at the shortest possible time scale. This is the time scale relevant to alternation (the length of individual IVIs), and it also captures provisioning CTH at longer time scales, whereas the reverse is not the case: the correlation between provisioning rates measured over longer periods will not necessarily detect provisioning CTH occurring on a shorter time-scale.

Simulation Model: Does Turn-Taking Affect the Frequency Distribution of IVIs (Including Producing a Refractory Period) and Provisioning CTH?

We investigated whether turn-taking in nest visits affects the frequency distribution of IVIs, including the presence of a refractory period, and provisioning CTH by simulating and analyzing sequences of provisioning visits by pairs. The sequences simulated either (i) independent provisioning by the male and female at a constant rate, λ, or (ii) turn-taking as in Johnstone et al.'s (2014) game theory model, in which each parent provisions at a constant rate, λ', when its mate is the last to provision, and does not provision when it itself is the last to provision. By simulating constant rates of provisioning, we excluded constraints due to foraging behavior and temporal environmental variation as sources of any refractory period and provisioning CTH that we found in the simulated sequences. Thus we can conclude that any difference between the presence of a refractory period or provisioning CTH between our two types of simulation must be generated by turn-taking.

Simulation Methods

Independent provisioning at a constant rate by the two parents was simulated by drawing a series of IVIs for each parent randomly from a negative exponential distribution with rate parameter, λ, equal to 1 min⁻¹ (Choosing a different value of λ is equivalent to changing the units in which time is expressed, and does not affect the calculated value of alternation or provisioning CTH). We used a similar procedure for the turn-taking simulations. However, turn-taking as in Johnstone et al.'s (2014) model assumes that each parent does not feed at the nest when it is its partner's turn to feed, so a single provisioning sequence for the pair combined was simulated by drawing a series of intervals from a negative exponential distribution with rate parameter λ', so that alternate intervals represent the periods when one of the sexes is provisioning and the other is not. We refer to these shorter intervals between one parent's nest visit and the other parent's nest visit as provisioning intervals (PIs) to distinguish them from IVIs, which are the intervals between successive visits by an individual parent. This means that the lengths of the IVIs for each parent consist of two consecutively drawn PIs summed together. For instance, the first IVI for one sex is equal to the sum of the first and second randomly drawn PIs, whereas the first IVI for the other sex is equal to the sum of the second and third randomly drawn PIs. Because each parent provisions for only half the time, we used a rate parameter, λ', of 2 min⁻¹, so that the overall rate of provisioning in each type of simulation was the same.

For each type of simulation we simulated 1,000 four-hour sequences in the R environment (version 3.2.3; R Development Core Team, 2017) using the R function rexp to generate intervals randomly drawn from a negative exponential distribution (see ‘Baldan_et_al_2019_R_script' in the Supplementary Material for the R script), and calculated alternation scores and provisioning CTH as described in Measures of alternation and provisioning CTH above.

Simulation Results and Conclusions

As expected, our simulation produced alternation scores close to zero for the parents provisioning independently (alternation score = 0.003 ± 0.090 (SD); Figure 1), and perfect alternation in all cases for the turn-taking parents. Unsurprisingly (because they were randomly drawn from a negative exponential distribution), the IVIs of the independently provisioning parents follow a negative exponential distribution [mean IVI = 0.995 min ± 0.995 (SD); Figure 2A; because the two parents are exactly equivalent in their provisioning behavior, the distributions of IVIs for males and females have been combined in the calculation of the mean and SD, and in the figures, for both the independently provisioning and turn-taking parents]. In contrast, the distribution of IVIs of the turn-taking parents has the same mean, but lower variance, and clearly has what Johnstone et al. (2014) would refer to as a refractory period (mean IVI = 0.996 min ± 0.705 (SD); Figure 2B). The refractory period arises because each IVI is the sum of two PIs drawn from a negative exponential distribution. Although intervals of zero length are the most frequent in a negative exponential distribution, an IVI of zero length can only be created by combining two PIs of zero length, whereas an IVI of closer to the average IVI length can be created by many different combinations of PI length. Thus, these intermediate length IVIs are more common than the shortest IVIs. Lastly, the provisioning CTH is close to zero for the independently provisioning parents [mean provisioning CTH = −0.003 ± 0.226 (SD); Figure 3A], while there is substantial provisioning CTH for the turn-taking pairs [mean provisioning CTH = 0.545 ± 0.226 (SD); Figure 3B]. This occurs because any specific PI forms part of the IVIs of each of the two sexes: for one sex the IVI is the sum of the specific PI and the previous PI, and for the other sex the sum of the specific PI and the following PI. The specific PI is chosen at random from a negative exponential distribution, so if it is short by chance, the overlapping IVIs will both be short by chance. The opposite applies if the specific PI is long by chance, so the lengths of the IVIs are correlated.

Figure 1. The frequency distribution of alternation scores from simulated provisioning sequences (n = 1,000) when parents provision independently.

Figure 2. The frequency distributions of inter-visit intervals from simulated provisioning sequences (n = 1,000) when parents provision independently (A), and take turns (B).

Figure 3. The frequency distributions of provisioning CTH (correlated temporal heterogeneity) from simulated provisioning sequences (n = 1,000) when parents provision independently (A), and take turns (B).

In conclusion, turn-taking by the parents may generate differences in the frequency distribution of IVIs, particularly by creating or enhancing the presence of a refractory period, compared with those when the parents feed independently of each other and only their foraging behavior determines the shape of the distribution. Moreover, turn-taking is a source of provisioning CTH at the time-scale of IVI length.

A Conceptual Framework: Partitioning the Sources of Alternation

One of the main questions faced by empirical studies of turn-taking is the extent to which evidence for turn-taking can be garnered from the sequences of provisioning visits by individual parents and pairs. In this section, we provide a conceptual framework (Figure 4) which clarifies the relationships between these variables and incorporates the results of our simulation model above, and use it to infer what conclusions can be drawn from provisioning data about the occurrence of turn-taking.

Figure 4. A conceptual framework for understanding the behavioral and ecological processes and factors affecting alternation. Green boxes are behavioral or ecological processes or factors; mauve boxes are statistical properties of the sequences of provisioning visits by individual parents or pairs. IVIs are inter-visit intervals, the intervals between successive provisioning visits to the nest by a given parent. The sequence of IVIs box contains a dashed box containing a question mark, because it is not clear whether turn-taking produces statistical properties in addition to provisioning CTH that depend on the sequence of IVIs and affect the level of alternation.

In creating a framework, it is important to make a distinction between behavioral and ecological processes on the one hand, and the statistical properties of the provisioning patterns that they affect or create on the other (Figure 4). Here we use the terms foraging behavior (meaning foraging constraints and decisions), parental coordination (including turn-taking), and temporal environmental variation for the behavioral and ecological processes, and the frequency distribution of the IVIs (which includes any so-called refractory period), sequence of the IVIs (which includes provisioning CTH) and alternation (the statistical property that we are trying to understand) for statistical properties of provisioning sequences. We use the broader term parental coordination in Figure 4 because there are other forms of parental coordination than turn-taking which can cause alternation in nest visits. In particular, if parents visit the nest site synchronously (Raihani et al., 2010; Mariette and Griffith, 2012; Bebbington and Hatchwell, 2016), and then enter the nest one at a time in a random order before departing together, their nest visits will alternate more than expected by chance on the basis of their provisioning rates. This is because the visit of the first parent to enter the nest will be an alternated visit by chance 50% of the time, whereas the visit of the second parent to enter the nest will always be alternated, giving an overall alternation rate of 75% as against the random expectation with non-synchronous visits at equal rates of 50%. Previous studies have recognized that besides turn-taking and other forms of parental coordination, foraging constraints producing a refractory period (and, more generally, any processes affecting the frequency distribution of IVIs; Johnstone et al., 2014), and temporal environmental variation producing provisioning CTH (Schlicht et al., 2016), can also create alternation by provisioning parents. However, by explicitly separating the statistical properties of provisioning sequences from the behavioral and ecological processes generating them, it becomes clear that other causal links may exist. In particular, we showed above using a simulation model that turn-taking can affect the frequency distribution of IVIs (including creating a refractory period; Figure 2) and provisioning CTH (Figure 3). In addition, temporal environmental variation is expected to increase the variance in provisioning rate, and hence to affect the frequency distribution of IVIs. Lastly, turn-taking produces alternation through its effects on the frequency distribution of IVIs and its effects on the sequence of IVIs. However, within the latter, it is not clear whether provisioning CTH is the only relevant statistical property that turn-taking produces, or whether there are additional statistical properties produced by turn-taking alone that depend on the sequence of IVIs and produce alternation. For this reason, Figure 4 contains a dashed box containing a question mark representing these possible additional statistical properties of the sequence of IVIs.

Partitioning the Alternation Score

In this section we use the conceptual framework to discuss how the alternation score can be partitioned, and what the resulting components of the alternation score represent. This enables us to elucidate what can—or cannot—be concluded from observations of parental provisioning sequences. We first consider Johnstone et al.'s (2014) randomization analysis, and then introduce a new analytical approach–the use of ‘pseudo-pairs'—which allows one of the components created by the randomization analysis to be further partitioned, giving an insight into the role of temporal environmental variation.

Randomization Analysis

Johnstone et al. (2014) quantified the amount of alternation between the members of a pair in search of evidence for the turn-taking behavior that their model predicted but were aware that alternation could be the result of other processes, specifically the presence of a refractory period. Their randomization analysis involves randomizing the sequence of IVIs for each of the two parents separately and then re-pairing the sequences. The amount of alternation remaining in the randomized provisioning data can be attributed only to the frequency distribution of the IVIs, and the reduction in the amount of alternation only to the sequence of IVIs. The randomization analysis therefore partitions alternation into components due to these two statistical properties. This partitioning is represented in the left-hand column of Figure 4.

Johnstone et al. (2014) implicitly assumed that these two components represent constraints on foraging behavior and turn-taking, respectively. Subsequently, it has been realized that the second component also includes the effects of provisioning CTH caused by temporal environmental variation (Schlicht et al., 2016; Santema et al., 2019). However, the conceptual framework in Figure 4 indicates that the behavioral and ecological processes involved in each of the two components are more complicated than this (see upper half of Table 1).

Table 1. The behavioral and ecological processes contributing to the components of the alternation score when it is partitioned using a randomization or pseudo-pair analysis.

First, the component due to the frequency distribution of IVIs can be affected not only by foraging behavior, but also by parental coordination including turn-taking and temporal environmental variation. Moreover, turn-taking and temporal environmental variation are expected to have opposite effects on the component of alternation due to the frequency distribution of IVIs. This is because a decrease in the variance of the IVI frequency distribution is expected to increase the amount of alternation. For example, when provisioning visits are made entirely regularly (i.e., the variance in the IVIs equals 0) the amount of alternation is the maximum possible. Turn-taking (at least in our simulation model above) decreases the variance of the IVIs, so increases the amount of alternation via its effect on the frequency distribution of IVIs (in addition to the effect via the sequence of IVIs). In contrast, temporal environmental variation will increase the variance of the IVIs, and decrease the amount of alternation (Table 1). This complexity in the processes and factors that may be involved makes the component due to the frequency distribution of IVIs hard to interpret. However, empirically this component tends to be small, although positive (Johnstone et al., 2014; Bebbington and Hatchwell, 2016; Savage et al., 2017; Table 2).

Table 2. Mean and standard error of the alternation score and its IVI sequence component in real pairs, randomized real pairs, pseudo-pairs, and randomized pseudo-pairs of great tits.

The second component of alternation discriminated by Johnstone et al.'s (2014) randomization test is due to the sequence of IVIs. This component can be generated by parental coordination (including turn-taking) and/or temporal environmental variation (Figure 4; Table 1), both of which have a positive effect on this component. As noted by others (Schlicht et al., 2016; Ihle et al., 2019; Santema et al., 2019), the size of this component therefore sets only an upper limit on the extent of turn-taking or temporal environmental variation, with each having a lower limit of zero. This component therefore does not provide evidence for the occurrence of turn-taking (or of temporal environmental variation). However, this component can be further partitioned by the analytical use of ‘pseudo-pairs' as described in the following section, and this further analysis does lead to the isolation of a component representing just one factor.

Pseudo-Pairs Analysis

The new analytical approach that we propose here partitions the component of alternation due to the sequence of IVIs into a component due to parental coordination and/or nest-specific temporal environmental variation, and a second component due to general environmental variation (lower half of Table 1; Figure 5). “Nest-specific temporal environmental variation” is environmental variation limited to individual nests. If it affects the provisioning rates of both members of the pair it will produce provisioning CTH and hence alternation. Nest-specific temporal environmental variation might include, for example, the begging behavior of the chicks or some kind of disturbance at or near the nest that is limited to that nest. “General temporal environmental variation” affects the provisioning rate of parents over a wider area (for example, the passage of rain showers) resulting in synchronous variation in provisioning rates at nests that are sufficiently close together to have experienced the same conditions. Because general temporal environmental variation causes synchronous variation, “pseudo-pairs” created by analytically pairing the provisioning sequence of a parent at one nest with the contemporaneous provisioning sequence of the opposite sex at a nearby nest will exhibit alternation. In pseudo-pairs, this is the only source of alternation in the component due to the sequence of IVIs. In real pairs, the component due to the sequence of IVIs also contains the effects of parental coordination and nest-specific environmental temporal variation (Table 1).

Figure 5. Partitioning the alternation score. To partition the alternation score, the alternation scores of real pairs (entire left-hand bar; purple) and pseudo-pairs (entire right-hand bar; red) are first calculated. Then the alternation scores of randomized real pairs and randomized pseudo-pairs are calculated. These estimate the component due to the IVI frequency distribution (lower section of left- and right-hand bars; darker shade) in real pairs and pseudo-pairs, respectively. The component due to the IVI sequence is then calculated by subtraction for real pairs and pseudo-pairs (upper section of left- and right-hand bars; lighter shade). Finally, the component due to parental coordination (including turn-taking) and/or nest-specific environmental variation is calculated by subtracting the component due to IVI sequence in pseudo-pairs (i.e., the component due to general environmental variation) from the component due to IVI sequence in pairs (blue, because it is calculated by subtraction using components in both pairs and pseudo-pairs). The middle bar shows the final partitioning of the alternation score in real pairs.

Calculating the Components of Alternation

The alternation score that we used above in our simulation models is additive on a linear scale (Baldan et al., 2019) and thus can be partitioned numerically. This partitioning is carried out as follows (Figure 5): (a) For real pairs and pseudo-pairs: (i) calculate the mean alternation score; (ii) calculate the mean alternation score for their respective randomized sequences (“randomized real pairs” and “randomized pseudo-pairs”). These are the components of alternation due to the frequency distribution of IVIs for real pairs and pseudo-pairs respectively; (iii) calculate the component of alternation due to the sequence of IVIs by subtracting the respective component due the frequency distribution of IVIs from the respective alternation score. (b) The component due to general temporal environmental variation in real pairs is equal to the component due to the sequence of IVIs in pseudo-pairs. (c) The component due to turn-taking and/or nest-specific temporal environmental variation in real pairs is equal to the component due to the sequence of IVIs in real pairs minus the component due to the sequence of IVIs in pseudo-pairs. In this way, the alternation score of real pairs is partitioned into three components due to the frequency distribution of IVIs (this component is difficult to interpret), parental coordination and/or nest-specific temporal environmental variation, and general temporal environmental variation.

Provisioning CTH (based on the Pearson rank correlation coefficient) is not additive on a linear scale, so cannot be partitioned numerically in the same way. Moreover, randomized sequences of IVIs should not exhibit provisioning CTH. However, provisioning CTH can be generated in real pairs by turn-taking and the two forms (nest-specific and general) of temporal environmental variation, while in pseudo-pairs it can only be generated by general temporal environmental variation. Thus, provisioning CTH in pseudo-pairs provides evidence for the occurrence of general temporal environmental variation, and a difference between provisioning CTH in real pairs and pseudo-pairs for the occurrence of parental coordination and/or nest-specific temporal environmental variation.

Great Tit Provisioning Sequences

We collected provisioning sequences from wild great tit Parus major pairs in order to carry out randomization and pseudo-pair analyses as described above, and thereby investigate the behavioral and ecological processes producing alternation and provisioning CTH in free-living provisioning parents.

Data Collection

We collected provisioning data from a great tit population at the Hoge Veluwe National Park, The Netherlands (52°23′N, 5°51′E) in 2014. This area contains around 400 nest-boxes that were checked weekly from the beginning of April to determine the onset of egg laying and incubation. Active nests were then checked daily from the day before predicted hatching to determine the exact hatch date (day 0). Parental nest visit data were collected from 17 broods with day 10 chicks on 10–13 May 2014. The recorded nests were 183–1,467 m apart (n = 29, median = 790 m; see Table S1 for individual distances). Nest visits were recorded using a small video camera with infra-red illumination mounted in the roof of the nest-box and connected to an external video recorder at the foot of the tree. The camera was placed on day 9, and the recordings (720 × 576 pixels of resolution) started before 09:30 of the following day. Recorders were synchronized daily (to the nearest second) with a digital watch before each recording began. Four hours of video (10:00–14:00) were scored for each nest, and the time the bird entered the nest-box (to the nearest second) and the sex of the provisioning parent (determined from the blackness of the crown feathers, which is more sexually distinct under infra-red than daylight illumination) were ascertained from the video-recordings. Provisioning videos were scored blind with respect to our study questions by MC (see Acknowledgments), who was not otherwise involved in this study. Data from radio-tagged pairs in a nearby great tit population scored by a different observer confirmed the potentially complete reliability of sexing based on crown feather coloration (n > 2,600 nest visits; DB, personal observations).

Provisioning Sequence Analysis

From the data for “real pairs” that we collected in the field, we created three additional kinds of provisioning sequence: “pseudo-pairs” were created by matching the provisioning sequence of one individual with that of another individual of the opposite sex recorded at a different nest at the same time on the same day. Our 17 real pair sequences were collected over 4 days (3, 4, 5, 5 pairs on 10–13 May, respectively) and generated 58 pseudo-pairs (pseudo-pairs per day: 6, 12, 20, 20). “randomized real pairs” were created by randomly rearranging the order of the IVIs of each parental sequence and then pairing the two newly created parental sequences, and “randomized pseudo-pairs” were created in the same way from pseudo-pairs sequences. In the analysis, we generated 10,000 replicate randomizations for each pair and pseudo-pair.

We calculated the alternation score and provisioning CTH as described above (see Measures of alternation and provisioning CTH). When calculating provisioning CTH, we included IVIs that fell within the same IVI as the preceding sampling point (see above; in the visit sequences for each of the parents in the 17 real pairs, 14% (n = 680) of the sampling points fell within the same IVI as the preceding sampling point). For the purposes of calculating the mean and SE of the alternation score or provisioning CTH for randomized real pairs and randomized pseudo-pairs we used the mean of the 10,000 randomization replicates for each of the real pairs or pseudo-pairs. Data for individual nest visits are given in ‘Baldan_et_al_2019_Dataset_1.xlsx', and alternation scores and provisioning CTH for each pair in each of the four pair types in ‘Baldan_et_al_2019_Dataset_2.xlsx', in the Supplementary Material.

Partitioning of the alternation score and provisioning CTH was carried out as described above (see Calculating the components of alternation and Figure 5).

Statistical Analysis

There are two general issues related to making statistical comparisons between alternation scores or provisioning CTH for different pair types. The first issue is that the way in which pseudo-pairs are created potentially gives rise to pseudo-replication: 58 pseudo-pairs were created from 17 real pairs. We addressed this problem by using mixed models in which we fitted male ID and female ID as random effects. The second issue is that for each randomized pair and randomized pseudo-pair, we generated 10,000 replicate randomizations. In statistical comparisons we used the mean value of these replicates for each randomized pair or randomized pseudo-pair. We calculated the IVI sequence component of alternation scores for each real pair or pseudo-pair by subtracting the mean of the matching randomized real or randomized pseudo-pair from the overall alternation score of the pair. Statistical tests on the components of the alternation score and provisioning CTH were carried out further as specified in Tables 2–4.

Table 3. Mean and standard error of provisioning CTH rate in real pairs, randomized real pairs, pseudo-pairs, and randomized pseudo-pairs of great tits.

Table 4. Statistical tests for comparisons between the alternation scores and provisioning CTH for different pair types.

To explore the scale at which general environmental temporal variation influenced the alternation score and provisioning CTH, we investigated how these variables varied with the distance between the nests of the members of a pseudo-pair. For each of the pairs of nests there are two “reciprocal” pseudo-pairs consisting of the male from each nest analytically paired with the female of the other nest. We therefore used linear mixed models with “reciprocal pair ID” as a random factor. If general environmental temporal variation is acting at a scale smaller than the maximum distance between pseudo-pairs, alternation score and provisioning CTH should decrease with distance between the members of a pseudo-pair. In addition, we investigated whether provisioning CTH varied between days. This could occur if the amount of temporal environmental variation varied between days, for example some days having short showers of rain, and others having less variable weather conditions. We investigated whether provisioning CTH between days occurred using a linear mixed model with reciprocal pair ID as a random effect and day (as a factor) as a fixed effect.

All the statistical analyses were performed in R environment (version 3.2.3). All mixed models were performed with the lme4 package (Bates et al., 2015) and the significance of the main effect was calculated with the Kenward-Roger approximation implemented in the pbkrtest package (Halekoh and Hojsgaard, 2014). All the statistical analyses were two-tailed, and significance was taken at α = 0.05. Unless otherwise stated, we reported means and standard errors of the estimates as mean ± SE.

Results

Alternation Scores

The mean alternation score of real pairs was 0.341, and significantly greater than zero (Table 2, Figure 6). The mean alternation score of randomized real pairs (the component due to the IVI frequency distribution) was much lower (0.028) and not significantly greater than zero (Table 2, Figure 6). The IVI sequence component of the alternation score in real pairs (the pairwise difference between the previous two values) was 0.313, and significantly greater than zero (Table 2, Figure 6). The alternation scores for pseudo-pairs exhibit the same general pattern, except that the alternation score for pseudo-pairs was not significantly greater than zero, but the IVI sequence component (generated from the pairwise difference between pseudo-pairs and randomized pseudo-pairs) was (Table 2, Figure 6). This last result indicates a significant effect of general temporal environmental variation on the alternation score. Lastly, the alternation score, and its IVI sequence component, for real pairs were significantly higher than the corresponding values for pseudo-pairs (Tables 2, 4, Figure 6), indicating a significant effect of parental coordination (including turn-taking) and/or nest-specific temporal environmental variation on the alternation score.

Figure 6. Alternation score (A) and provisioning CTH (B) in real, randomized real, pseudo- and randomized pseudo-pairs. Points designate the mean ± SE. For pairs and pseudo-pairs, these statistics are based on the value for each pair (n = 17) or pseudo-pair (n = 58). For randomized pseudo-pairs and randomized real pairs, mean ± SE are based on the mean values of the 10,000 randomized distributions for each of the pairs and pseudo-pairs. An alternation score of zero (dashed line) represents the amount of alternation expected by chance, assuming that the probability of a provisioning visit by each parent is constant with respect to time.

We used the values of the alternation scores in Table 2 to calculate the components due to the IVI frequency distribution (0.043), nest-specific IVI sequence effects (0.251), and general environmental IVI sequence effects (0.062) (see Calculating the components of alternation and provisioning CTH and Figure 5). Based on these values, approximately 8% of the alternation score is due to the IVI frequency distribution, 74% to nest-specific IVI sequence effects, and 18% to general environmental IVI sequence effects.

Provisioning CTH

The values of provisioning CTH (Table 3) follow broadly the same pattern as the alternation scores in that provisioning CTH was significantly higher in real pairs than pseudo-pairs (Tables 3, 4). As expected, the values of provisioning CTH for randomized real pairs and randomized pseudo-pairs were close to zero, although the latter was significantly negative (see Calculating provisioning CTH above for a discussion of this negative bias).

Variation in Alternation Score and Provisioning CTH in Pseudo-Pairs With Distance Between the Members of Pseudo-Pairs and Between Days

Both alternation score and provisioning CTH of pseudo-pairs tended to increase with distance between the members of pseudo-pairs (Figure 7). This increase was significant for provisioning CTH [F_{(1, 27)} = 4.68, p = 0.040], but not for alternation score [F_{(1, 27)} = 1.10, p = 0.302].

Figure 7. Alternation score (A) and provisioning CTH (B) in pseudo-pairs in relation to the distance between the members of pseudo-pairs. Vertical bars join data for the two reciprocal pseudo-pairs for each pair of nests (the male of each nest pseudo-paired to the female of the other nest). Red, green, cyan and mauve points and lines indicate data collected each day from 10 to 13 May 2014, respectively. The shaded gray areas are the standard errors for the relationship with distance from a linear mixed model.

Alternation score, but not provisioning CTH, of pseudo-pairs varied significantly between the 4 days on which we recorded provisioning behavior [alternation score: F_{(3, 25)} = 3.49, p = 0.030. provisioning CTH: F_{(3, 25)} = 1.37, p = 0.275].

We expected a negative relationship between alternation score or provisioning CTH and distance between the members of pseudo-pairs, but the observed relationships are positive (significantly so for provisioning CTH). However, in our data, the separation between nests involved in pseudo-pairs varied between days (Kruskal-Wallis H test, χ² = 8.326, p = 0.040, df = 3) and decreased over the 4 days of data collection (Kendall rank correlation coefficient, τ = −0.431, p = 0.003, n = 29). These relationships resulted from the limited availability of day 10 broods. We do not have an adequate sample size to separate an effect of distance from the confounded effects of day or date, so have not attempted to carry out such an analysis. However, our uncontrolled results give little reason to expect a negative relationship for either variable over the range of distances that we studied if we had been able to do so.

Discussion

In this study, we investigated the processes creating alternation of nest visits. First, we found using a simulation model that turn-taking can produce a refractory period and provisioning CTH. Second, we created a conceptual framework that allows us to partition alternation into three components: the effects of the IVI frequency distribution; nest-specific effects of the sequence of IVIs; and general effects of the sequence of IVIs. Third we applied this analytical method to provisioning data from wild great tits and estimated that nest-specific sequence effects are the major contributors to alternation explaining 74% of the observed amount of alternation. The frequency distribution of IVIs and general sequence effects are responsible for the remaining amount of alternation (8 and 18%, respectively).

The first of the three components—the effect of the frequency distribution of IVIs—is quantified by a randomization test suggested by Johnstone et al. (2014). They implicitly assumed that this component was due to what they refer to as a refractory period caused by constraints on foraging behavior, and other authors have also equated a refractory period with foraging constraints (Schlicht et al., 2016; Savage et al., 2017; Ihle et al., 2019; Santema et al., 2019). However, our simulation and conceptual model indicate that this notion is ill-founded, and that this component may also include the effects of parental coordination (including turn-taking) and temporal environmental variation. This makes this component difficult to interpret, but fortunately the actual size of this component is small in empirical studies, including our present study of great tits (Johnstone et al., 2014; Bebbington and Hatchwell, 2016; Savage et al., 2017).

Our second and third components of alternation are derived from the other component produced by Johnstone et al.'s (2014) randomization test—the IVI sequence effect—using our pseudo-pair analysis to separate nest-specific and general sequence effects. Our second component—due to nest-specific IVI sequence effects—resembles the IVI sequence effect discriminated by Johnstone et al.'s (2014) randomization test in being generated by parental coordination (including turn-taking) or (in this case, nest-specific) temporal environmental variation. Other authors (Ihle et al., 2019; Santema et al., 2019) have argued for the unpartitioned IVI sequence effect that a role for turn-taking can only be demonstrated when all relevant environmental variation has been taken into account, and that this is a difficult-to-impossible task given the ubiquity of such variation. This argument also applies to our nest-specific IVI sequence component.

Our third component was isolated using analytically created pseudo-pairs. Because the members of such pairs are separated by a minimum of 183 m in our study of great tits, it is implausible that these individuals are reacting directly to each other, so that the significant IVI sequence effect in pseudo-pairs can only be attributed to a common response by the parents to some general environmental variation and cannot be due to turn-taking. Moreover, this general temporal environmental variation accounts for about a fifth of the alternation in real pairs. Because nest-specific sources of temporal environmental variation creating alternation may be as, or more, frequent than general sources, temporal environmental variation may account for a substantial proportion of the observed amount of alternation, leaving little room for turn-taking. Experimental manipulations of parental behavior may be more effective approaches than analyzing observational provisioning data to investigate the existence of turn-taking. Manipulations could involve brood size (Baldan et al., 2019; Griffioen et al., 2019a), begging playback at the nest (Santema et al., 2017) or handicapping of one of the parents (Griffioen et al., 2019b). In particular, manipulations directed at only one parent (e.g., selective playback or handicapping) can be used to investigate whether parents do, indeed, react to the provisioning behavior of their mate.

The significant IVI sequence effect in the alternation of pseudo-pairs raises some intriguing questions beyond the implications for the possible occurrence of turn-taking. The first of these is simply whether our result can be replicated, and further questions concern the scale over which this effect occurs. Pseudo-pairs whose members are closer together are more likely to share more temporal environmental variation (e.g., due to local rain showers or shared food patches between breeding pairs), so at some spatial scale we expect provisioning CTH and alternation in pseudo-pairs to decrease as the separation between members of pseudo-pairs increases, and the scale at which this happens may give some indication of the nature of the temporal variation. This is also the case for variation between days in the level of alternation in pseudo-pairs because higher levels of an IVI sequence effect in pseudo-pairs are expected to occur on days with higher levels of the relevant temporally varying environmental variable. We were unable to detect any decrease in the alternation scores of pseudo-pairs with separation between the members of the pairs of up to 1,467 m. However, our data did not allow us to control adequately for any confounded effects of date or days, and we encourage others carrying out similar studies to avoid the distance between pseudo-pairs being confounded with date and to obtain larger sample sizes both within days and in the number of days sampled. We did find significant variation between days in the amount of alternation in pseudo-pairs. Possibilities for further work on alternation and provisioning CTH between pseudo-pairs include investigating at what distance these can still be detected, and whether there is any relationship between alternation or provisioning CTH of pseudo-pairs on individual days and meteorological conditions.

Schlicht et al. (2016) and Santema et al. (2019) have emphasized the possible role of provisioning CTH generated by a common response to temporal environmental variation in producing alternation. However, as far as we are aware, we are the first to attempt to quantify provisioning CTH by calculating the correlation between the provisioning rates of the parents across time. We found an appreciable and significant level of provisioning CTH in the provisioning sequences of great tit parents which at first sight appears to argue that temporal environmental variation does indeed play a major role in generating alternation as the above authors have suggested. However, we also showed that turn-taking can, by itself, generate considerable provisioning CTH, so its presence does not provide evidence against turn-taking as the process behind alternation. It may simply be that provisioning CTH and alternation are inextricably statistically linked, and that any process generating alternation automatically generates CTH, and vice versa. Similarly, it is also currently unclear whether there are any additional sequence effects than provisioning CTH that are linked to, or generate, alternation.

In conclusion, our study contributes in multiple ways to understanding the kinds of conclusions that can be reached from empirical measurements of alternation rates. First, we have shown using simulations that turn-taking can, by itself, generate a refractory period and provisioning CTH. Second, we created a conceptual framework and combined it with a novel analytical approach using pseudo-pairs to provide a means of partitioning alternation into three components. Finally, we applied this new framework and methodology to data on parental great tits provisioning in the wild. This analysis demonstrated for the first time that alternation by provisioning parents can be generated by temporal environmental variation to which both parents respond, accounting, in the great tits, for about one fifth of the observed amount of alternation.

Data Availability Statement

All datasets generated and analyzed for this study are included in the article/Supplementary Material.

Ethics Statement

The animal study was reviewed and approved by KNAW Dier Experimenten Commissie (DEC).

Author Contributions

CL conceived the study and created the conceptual framework. DB and CH suggested additional analyses. DB carried out the simulation and collected the field data and performed the analysis with input from CH and CL. DB and CL drafted the manuscript. All authors contributed to revisions and approved the final manuscript.

Funding

This research was supported by a grant (823.01.005) from the Netherlands Organization for Scientific Research (NWO) to CL.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the board of the National Park De Hoge Veluwe for the permission to conduct our research there. We are grateful to Henri Bouwmeester for his help with the field work, Jorina Boer for her help collecting data, Mathias Cox for scoring the video recordings, and James Savage for providing the script for the randomization test. We thank Jenny Ouyang, James Savage and Marcel Visser for comments on various versions of this manuscript, and Antica Culina and Arie van Noordwijk for discussion. CL thanks Erna Harthoorn for discussion and support.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2019.00448/full#supplementary-material

References