Abstract In most species, males do not abandon offspring or reduce paternal care when they are cuckolded by other males. This apparent lack of adjustment of paternal investment with the likelihood of paternity presents a potential challenge to our understanding of what drives selection for paternal care. In a comparative analysis across birds, fish, mammals, and insects we identify key factors that explain why cuckolded males in many species do not reduce paternal care. Specifically, we show that cuckolded males only reduce paternal investment if both the costs of caring are relatively high and there is a high risk of cuckoldry. Under these circumstances, selection is expected to favour males that reduce paternal effort in response to cuckoldry. In many species, however, these conditions are not satisfied and tolerant males have outcompeted males that abandon young.

Author Summary In most species where it has been studied, males do not abandon or reduce paternal care when they are cuckolded by other males. These observations have presented a long-standing challenge to our understanding of what drives selection for paternal care. Our analysis of cuckolded fathers from 50 species of birds, fish, mammals, and insects, however, shows that sometimes it pays for males to stick around. In the case of humans and burying beetles, this is because females are relatively monogamous—by deserting, it is most likely the case that fathers will be deserting their own young. In species such as the chacma baboon, males face a significant risk of cuckoldry, and face potentially high penalties in terms of future breeding success by wasting precious resources on the young of other males. Unlike in humans, promiscuous females in these species will almost certainly lose the support of her mate in the effort to raise her young to adulthood.

Citation: Griffin AS, Alonzo SH, Cornwallis CK (2013) Why Do Cuckolded Males Provide Paternal Care? PLoS Biol 11(3): e1001520. https://doi.org/10.1371/journal.pbio.1001520 Academic Editor: Allen J. Moore, University of Georgia, United States of America Received: June 6, 2012; Accepted: February 12, 2013; Published: March 26, 2013 Copyright: © 2013 Griffin et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Funding: ASG is supported by the Royal Society (royalsociety.org); CKC by a Browne fellowship at The Queen's College, Oxford University (www.queens.ox.ac.uk), and the Swedish Research Council (Vetenskapsrådet: www.vr.se); SHA by a visiting fellowship from Corpus Christi College, Oxford University (www.ccc.ox.ac.uk), the National Science Foundation (www.nsf.gov), and Yale University (www.yale.edu). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist. Abbreviations: CI, credible interval; SE, standard error

Introduction Parental care is demanding: the effort it takes a typical garden bird to raise a clutch of chicks to adulthood is equivalent, in human terms, to cycling the Tour de France [1]. Intuition suggests that a male should only embark on this feat if he is the father of the chicks in his nest—if he has been cuckolded, he should avoid wasting resources on enhancing the reproductive success of his rivals and reduce paternal effort. Forty years of empirical research, however, have failed to provide consistent support for this prediction [2]. While there are exceptions, in the majority of species, it is reported that males do not significantly decrease care when cuckolded [2]–[5], challenging our understanding of how natural selection favours males that provide parental care. If information about paternity is unavailable or unreliable, selection may favour males that continue to provide care to avoid potential costs of abandoning their own offspring [6],[7]. Consequently, attempts to test explanations for the lack of paternal care adjustment have focussed on the ability of males to accurately assess paternity, but this has yielded abundant unexplained variation between species [2]. Another explanation is that variation in adjustment reflects differences in the strength of selection on males to reduce paternal care in response to loss of paternity. Firstly, theory predicts that the behaviour of males should optimise lifetime reproductive success rather than just paternity [8],[9]. This is formalised in Hamilton's Rule, rb−c>0, where r is the relatedness between the caring male and offspring in this context, b is the fitness benefit of care to offspring, and c is the cost of care to future reproductive success of the male [10]. Cuckolded males are, therefore, predicted to be relatively tolerant of cuckoldry if paternal care does not reduce future reproductive success (low c) [8],[11],[12], and/or paternal care has little effect on offspring fitness (low b) [10],[13]. When b is low, variation in r has relatively weak effect on variation in selection (if we substitute b = 0 into rb−c>0, rb is always 0; if we substitute b = 1 into rb−c>0, rb varies from −1 to 1 depending on the value of r). Secondly, when there is little variation in the risk of cuckoldry between breeding attempts, or cuckoldry is rare, there will be weak selection for adjustment [11]. This is because rare or low variation in cuckoldry within males reduces the opportunity for selection to favour individuals that adjust paternal care [14]. Despite well-developed predictions about when cuckolded males should reduce care, our understanding of this problem remains limited because empirical studies have focused on reporting the presence or absence of adjustment without further formal analyses of causation.

Discussion In this study, we address problems arising from the lack of empirical tests of theoretical predictions about the evolution of paternal care adjustment. In particular, measuring variation in paternity and the costs and benefits of paternal care within species is a major undertaking and it remains a challenge to think of ways how these factors can be experimentally manipulated to test existing theoretical predictions. By adopting a comparative approach, we have been able to exploit variation between species and our results suggest several important considerations for future studies of paternal care adjustment. Firstly, our study provides some guidance about the characteristics of species best suited to studies of paternal care adjustment. Specifically, the expectation of adjustment should be lowered in species with either low costs of care or low promiscuity. Secondly, we suggest that studies more closely address theoretical predictions by linking measures of adjustment with measures of costs of care (none of the measures in our adjustment dataset were directly linked with the measures of costs [Tables S1 and S2]). Thirdly, the interaction between promiscuity and cost reported here could be tested empirically within a single species by characterising the relationship between the strength of adjustment and the residual reproductive value of males. This could be achieved by exploiting individual variation in factors such as male age, timing of breeding or quality, which are expected to correlate with residual reproductive value. When we see males caring for the offspring of another male, it is possible to assume that they are doing so simply because of a failure to accurately determine paternity success. Although our analysis shows that males are not as constrained by lack of reliable cues to paternity as often thought, it remains the case that females seem to get away with promiscuous mating in some species, with cuckolded males continuing to provide care. Are cuckolded males that maintain care blissfully ignorant or selfless dupes? Our analyses suggest they are neither; instead, whether or not cuckolded males reduce paternal care is readily explained by the cost of paternal care and the risks of cuckoldry. More generally, cuckolded males provide a good example of how selection favouring tolerance can lead to the appearance of losing in an evolutionary arms race with cheats [16]. The finding that males show some degree of response to cuckoldry in the majority of species studied has implications for the evolution of paternal care more generally: we estimate that the extent to which females reduce paternal investment by engaging in promiscuous copulations is 12% on average within-species. By changing relatedness between nest mates, it has been shown that female promiscuity can drive transitions to and from cooperative breeding in birds [17]. By lowering relatedness between males and offspring, female promiscuity also has the potential to drive selection for reduced levels of paternal investment [18]–[20] and may ultimately cause the breakdown of biparental breeding systems [21],[22].

Materials and Methods Data Collection Adjustment of paternal care. We quantified the adjustment of paternal care by calculating the effect size of the relationship between paternal care and paternity (r Adjust ) for all species with available data using the effect size calculator in metawin on test statistics reported in papers. We searched the Web of Science for published papers that included the keywords “care” and “paternity” and then conducted forward and backward searches through cited references of these papers (search performed on papers published up to and including 1 March 2011). We also contacted researchers within the field to check for availability of unpublished data. In total, we obtained 192 effect sizes from 62 papers representing 48 species, from 29 families and five different classes. For each effect size we recorded: (a) whether the amount or the probability (i.e., desertion) of paternal care was measured (two-level factor), (b) whether the data were observational or experimental (two-level factor), (c) whether molecular genetic techniques had to been used to assess parentage (yes/no: two-level factor), and (d) the potential cues available to males the study had measured, which included access to females during the fertile period (yes/no: two-level factor) and the presence of potential competitors (yes/no: two-level factor). For the species (n = 48) with data on the adjustment of paternal care we then collected information on the benefits of male care for offspring, the costs of care for male future reproductive success, and the risk of that females are promiscuous. Measuring the benefit of male care for offspring. To assess the benefits of paternal care for offspring we calculated effect sizes from studies that measured the relationship between male care and offspring fitness (r Benefit ; Figure S1). We located data by checking the references of the studies on the adjustment of paternal care, contacting authors working on target species, and searching the Web of Science. In the Web of Science we searched using the latin name of the species as a keyword and if the number of results were over 200 we also added the term “care OR paternal”. For each effect size we also recorded (a) if the probability or amount of paternal care had been measured (two-level factor), (b) whether the study was observational or experimental (two-level factor), and (c) how offspring fitness was measured (three-level factor: condition, survival, or recruitment). Measuring the costs of care to male future reproductive success. To quantify the cost of paternal care to male future reproductive success (opportunity costs) we calculated the effect size of the relationship between male care in a current breeding attempt with measures of success in the future (r Cost ; Figure S2). We located studies using exactly the same methods as those used to find studies on r Benefit . Once again we recorded (a) if the probability or amount of paternal care had been measured (two-level factor), (c) whether the study was observational or experimental (two-level factor), and (c) how male future success was measured (three-level factor: attracting mates, reproductive success once paired [number of offspring] or survival). We also recorded the type of care that males provided in studies measuring r Benefit , r Cost , and r Adjust , but it was not possible to analyze if this influenced effect sizes because it was confounded with taxonomic grouping. For example, only male fish care for offspring by fanning eggs. Measuring the risk of female promiscuity. We quantified the probability that males care for unrelated offspring using data from molecular genetic studies that have measured levels of multiple paternity. Previously we have collected data on birds [17], but for other species we located data using the same methods as those used to find studies on the benefits and costs of paternal care. We defined multiple paternity as the percentage of families in the population that had offspring sired by more than one father (range = 6.1%–97.5%). We also collected data on the percent offspring sired by males other than the caring male (range = 5.3%–71.2%). However, we used the probability of males being cuckolded as our measure of relatedness to offspring because it captures variation in multiple paternity at the level of the population, more data were available and it is highly correlated to the percentage of offspring fathered by other males (Pearson's correlation coefficient = 0.84). General statistical techniques. We analysed variation in r Benefit , r Cost , and r Adjust with taxonomic random effect meta-analyses performed using Bayesian linear mixed models (BMM) with Markov chain Monte Carlo estimation in MCMCglmm, R version 2.13.1 [23],[24]. Effect sizes were transformed to Fisher's Z (Zr) before analysis and weighted by the inverse variance to account for variation in sample size between studies. The variance associated with effect sizes was calculated as the reciprocal of the sum of the conditional variance, where n is the sample size of the study [25]. Prior to analyses covariate fixed effects were Z-transformed (mean = 0, standard deviation = 1) and two-level fixed effects were converted to binary coding −1, 1 so that the magnitude of parameter estimates could be compared and main effects could be interpreted in the presence of higher order interactions [26],[27]. The data in all analyses were from a taxonomically diverse range of species and for some species there were multiple studies and multiple effect sizes per study. We took into account the non-independence of data arising from multiple studies on the same species, and from the phylogenetic relationships between species by defining a nested random effects structure whereby study was nested within species, species was nested within family, and family was nested within class [23]. Only family and class were entered into models because in our dataset there was no replication at the level of genus and order. We ran each analysis for 3,000,000 iterations with a burn-in of 2,500,000 and a thinning interval of 100. This generated 5,000 samples from each chain from which model statistics (deviance information criteria [DIC], posterior mean ± standard deviation [SD], posterior mode and 95% CIs [lower CI–upper CI]) were calculated. Terms were considered statistically significant when 95% CIs did not span 0 and pMCMC values calculated in MCMCglmm (number of simulated cases that are >0 or less than 0 corrected for finite number of MCMC samples) were less than 0.05 [24]. Initially we tried using two different priors. First, we used an inverse gamma prior that is commonly used for random effects (V = 1, nu = 0.002). Second, we ran models with parameter expanded priors (half-Cauchy priors following [28]: V = 1, nu = 1, alpha.mu = 0, alpha.V = 25∧2) due to some variance components being close to 0. The inverse gamma prior led to better convergence as measured by the Gelman-Rubin statistic (see below) and produced almost identical results to equivalent frequentist models run with ASReml-R version 3 [29]. We therefore used priors with V = 1, nu = 0.002 for all models. Missing values in explanatory variables were imputed using the mean of the missing variable (Z-transformed scale = 0) so that it was possible to compare different models using DIC [30]. We checked the convergence of each analysis using two diagnostic tests in the R package “coda” [31]. First, we ran each analysis three times and used the Gelman-Rubin statistic (potential scale reduction factor [PSR]) to compare within- and between-chain variance [32]. When convergence is met PSR<1.1 and in all our analyses PSR was less than 1.01. Second, we used Geweke's convergence diagnostic, which calculates Z scores from mean parameter estimates ± standard errors (SEs) generated from the first 10% and the last 50% of the chain [33]. If Z scores follow an asymptotically standard normal distribution the samples are considered to be drawn from a stationary distribution. Tests for publication bias. We tested whether there was evidence for publication bias/small-study effects (smaller studies show greater effects than larger studies) in estimates of r Benefit , r Cost , and r Adjust using trim and fill analysis and funnel plot asymmetry tests. Trim and fill analysis. We conducted trim and fill analyses on r Benefit , r Cost , and r Adjust . It is currently not possible to implement trim and fill methods in the Bayesian mixed models that were used to conduct the meta-analyses. We therefore performed the trim and fill analysis using a random-effects meta-analysis with restricted maximum likelihood estimation (REML) in the R package “metaphor” [34]. It must be noted that this does not allow for the non-independence of multiple effect sizes per study and taxonomy to be taken into account, but still provides a test of the influence of small-sample effects across all data points on the mean effect size. For r Benefit and r Adjust the trim and fill analyses estimated that no points were missing due to small-study effects (Figure S3a and S3b). For r Cost 6 points were estimated to be missing resulting in an estimated effect size of 0.12 (95% confidence limits: 0.03–0.22) rather than 0.19 (95% confidence limits: 0.11–0.26) (Figure S3c). Together these results indicate that at the level of individual studies r Benefit and r Adjust show little evidence of being influenced by small-study effects and that significance of the mean effect size of r Cost is not changed by publication bias. Funnel plot asymmetry. We carried out regression tests of funnel plot asymmetry using Egger's regression method on the full datasets used in the meta-analyses [35]–[37]. Egger's regression tests the relationship between the effect size and its associated SE, weighted by the inverse of the estimated effect size variance. In order to control for the non-independence of effect sizes from the same study and taxonomy we used Bayesian mixed models with the same random effects structure and settings (priors, iterations, burn-ins, thinning, and convergence testing) that were used for the meta-analyses. Effect sizes were entered as the response variable weighted by 1/se2 and se was entered as a covariate. For r Benefit , r Cost , and r Adjust the confidence intervals of the relationship between effect size and SE all spanned 0 and pMCMC>0.80. This further indicates that small-sample effects did not have a strong influence on the estimation of r Benefit , r Cost , or r Adjust . Specific analyses. We conducted five sets of meta-analyses: (1) First, we tested if variation in r Benefit was explained by methodological (fixed effects: amount versus probability, observation versus experiment, measure of offspring fitness) and taxonomic/study effects (random effects) (Tables S4a–S4f); (2) Second, we tested if variation in r Cost was explained by methodological effects (fixed effects: amount versus probability, observation versus experiment, measure of cost), the benefits males provide to offspring (r Benefit and proportion of male care: covariate fixed effects) and taxonomic/study effects (random effects) (Tables S5a–S5h); (3) Third, we tested if r Adjust was explained by methodological effects (fixed effects: amount versus probability, observation versus experiment, use of genetic methods, male access to females, presence of competing males) (Tables S6a–S6h); (4) Fourth, we tested if r Adjust was explained by biological variables (fixed effects: r Benefit [covariate], r Cost [covariate], multiple paternity and multiple paternity [covariates] [23]), and taxonomic/study effects (random effects) after controlling for any methodological differences found in analyses 1–3 (Tables S7a–S7h); (5) Fifth, we tested if r Adjust and the effects found in analyses 1–4 were influenced by phylogenetic history not accounted for by the nested taxonomic random effects structure. To do this we had to restrict the data to only birds (32 species) for which there is a phylogeny. We used a bird supertree [17] that was pruned to include only species for which there were data on r Adjust using the “ape” package in R [38]. We used a Phylogenetic Bayesian mixed model (BPMM) that was exactly the same as in analysis 4, but instead of fitting a nested taxonomic random effects structure we fitted the phylogeny and study as random effects. For analyses 1–4 we used the following procedure. First, we tested fixed effects individually to estimate the mean effect of each variable. Second, we ran full models to estimate the marginal effects of each variable controlling for the effects of all other variables. Third, we ran full models sequentially adding in all biologically relevant interactions (r Cost * r Benefit to test if adjustment is increased by the multiplicative effects of costs and benefits: r Cost * multiple paternity to test if adjustment is greater when costs are high and there is a high chance that males are cuckolded: r Benefit * multiple paternity to test if adjustment is higher in species where males confer greater benefits to offspring and the risk of cuckoldry is high) and then selected the best model using DIC values. For each set of analyses we present a model summary table that details each model fitted, its DIC value, and the percentage of variation in effect sizes explained by each random term (Tables S4a, S5a, S6a, and S7a). The percentage of variation explained by each random effect, in conjunction with the estimates of variance components, gives an indication of whether different taxonomic levels and different studies share a common effect size (see Horvathova et al. [39] for similar approach to quantifying heterogeneity). We then present tables for each analysis with estimates of all fixed and random effects (Tables S4b–S4f, S5b–S5h, S6b–S6h, S7a–S7h). In the main text we transform effect sizes from Zr back to r and present posterior means, credible intervals and pMCMC values.

Acknowledgments This work would not have been possible without the generosity and expertise of the many researchers who supplied us with data or advised us on the interpretation of their data: Jim Bednarz, Daniela Canestrari, Mike Cherry, Javier Cuervo, Peter Dunn, John Faaborg, Jo Frommen, Peter Grant, Ian Hartley, Sjouke Kingma, Kia Lindström, Ken Levenstein, Matt Low, Rob Magrath, Andrea Manica, Imane Oubou, Nichola Raihani, Andy Radford, Michael Richardson, Tim Schmoll, Kirk Stodola, Ola Svensson, Rose Thorogood, Ben Sheldon, Per Smiseth, Juan Soler-Cruz, Alberto Velando, Geoff While, Dean Williams, and Derek Yalden. Mike Wilson helped us extract data from papers in written in German and Russian that we were able to access thanks to the collections of the Alexander Library. Ben Sheldon, Stuart West, Geoff While, and Hanna Kokko greatly improved the clarity of our manuscript.

Author Contributions The author(s) have made the following declarations about their contributions: Conceived and designed the experiments: ASG SHA CKC. Performed the experiments: ASG SHA CKC. Analyzed the data: ASG SHA CKC. Contributed reagents/materials/analysis tools: ASG SHA CKC. Wrote the paper: ASG SHA CKC.