{'fraction_non_alphanumeric': 0.0322282264738815, 'fraction_numerical': 0.019862713597195852, 'mean_word_length': 5.247414841849149, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Predicting evolutionary responses to selection on polyandry in the wild: additive genetic covariances with female extra-pair reproduction

The evolutionary forces that underlie polyandry, including extra-pair reproduction (EPR) by socially monogamous females, remain unclear. Selection on EPR and resulting evolution have rarely been explicitly estimated or predicted in wild populations, and evolutionary predictions are vulnerable to bias due to environmental covariances and correlated selection through unmeasured traits. However, evolutionary responses to (correlated) selection on any trait can be directly predicted as additive genetic covariances (covA) with appropriate components of relative fitness. I used comprehensive life-history, paternity and pedigree data from song sparrows (Melospiza melodia) to estimate covA between a female's liability to produce extra-pair offspring and two specific fitness components: relative annual reproductive success (ARS) and survival to recruitment. All three traits showed non-zero additive genetic variance. Estimates of covA were positive, predicting evolution towards increased EPR, but 95% credible intervals overlapped zero. There was therefore no conclusive prediction of evolutionary change in EPR due to (correlated) selection through female ARS or recruitment. Negative environmental covariance between EPR and ARS would have impeded evolutionary prediction from phenotypic selection differentials. These analyses demonstrate an explicit quantitative genetic approach to predicting evolutionary responses to components of (correlated) selection on EPR that should be unbiased by environmental covariances and unmeasured traits.


INTRODUCTION
Extra-pair reproduction (EPR) by socially monogamous females, and the resulting polyandry, could substantially alter distributions of reproductive success and the evolutionary dynamics of sexually selected traits compared with those arising given strict monogamy [1][2][3]. The evolutionary forces driving polyandry and EPR are therefore of broad interest but remain far from clear [1,2,[4][5][6][7][8]. As for any phenotypic trait, continued evolutionary change (or stasis) in EPR will result from the combined magnitudes of all components of selection coupled with the additive genetic (co)variances underlying EPR and all fitness components and other traits through which selection acts [5,[9][10][11]. Understanding the evolutionary forces driving EPR and polyandry ultimately requires all such selection components and genetic (co)variances to be estimated in populations experiencing natural genetic and environmental variation in reproductive strategy and fitness.
Quantifying the form and magnitude of selection on any trait, and predicting any evolutionary response, is difficult in wild populations where environmental covariation among traits and fitness components may be substantial and all relevant traits, environmental covariates and individuals are unlikely to be measured [11,[12][13][14][15][16]. Some relevant parameters can be estimated for single or multiple traits given sufficient data describing variation in phenotype and fitness components among relatives. These include phenotypic selection differentials and gradients with respect to specific fitness components, heritabilities and additive genetic (co)variances [9,10,17,18]. The evolutionary response to selection over one generation attributable to these fitness components can then be predicted using univariate or multivariate breeder's equations [9,10,17].
The approach of estimating phenotypic selection and genetic (co)variances and then predicting consequent evolution is valuable due to its relative empirical tractability and intuitive interpretation, but has limitations [10 -13,15]. Widely recognized problems are that phenotypic covariances between focal traits and fitness components, and hence estimated selection differentials, can reflect environmental covariances or indirect effects of other correlated traits rather than direct causal effects. Predicted evolutionary responses can consequently be substantially biased [13 -16,19 -23].
These problems can in theory be resolved by measuring all relevant environmental covariates and traits, and estimating multivariate (partial) selection gradients with respect to fitness components of interest [9,11,13,14]. However, all such covariates and traits are unlikely to be measured in practice [9,13,15,16,21]. In particular, correlated selection through individuals that die before trait expression, or on reproductive traits through adults that fail to breed, could substantially bias predictions of evolutionary change (creating so-called 'invisible fractions' [14,24,25]).
One way to eliminate bias in evolutionary predictions due to environmental covariances and unmeasured traits is to predict responses to selection by directly estimating additive genetic rather than phenotypic (co)variation [13][14][15]20,21]. Specifically, the predicted evolutionary response to directional selection over one generation, defined as the change in mean breeding value for focal trait z (DA z ), equals the additive genetic covariance (cov A ) between the trait and relative fitness (w); DA z ¼ cov A (w,z) (the 'secondary theorem of selection' [11,20,26]). Direct estimation of cov A (w,z) can therefore allow unbiased prediction of per-generation evolutionary responses to components of directional selection when environmental variation cannot be controlled or completely measured [14,15,20,21]. Such analyses also allow environmental covariances between traits and fitness components to be explicitly quantified, thereby indicating the degree to which assumptions underlying standard breeder's equation predictions are violated [15,20,21]. An explicit quantitative genetic approach can also allow responses to correlated selection (for example, through individuals that die before trait expression or fail to breed and express reproductive traits) to be predicted because cov A between focal traits and relevant fitness components can be estimated from observed trait values of relatives [14,24].
Despite the potential valuable insights, evolutionary responses to selection have rarely been predicted by estimating cov A between traits and relative fitness components in wild populations, and this approach has not been implemented for female EPR or other forms of polyandry. Even phenotypic selection differentials and gradients on female EPR have not been explicitly estimated with respect to major components of female fitness in the wild. Instead, direct costs and benefits of EPR are typically postulated and estimated in currencies such as decreased paternal care by a female's cuckolded social mate [5] (see §4). No clear consensus regarding such costs and benefits has emerged, not least because observed associations may reflect environmental covariation rather than direct causal effects [27,28]. Explicit estimates of the form, direction and magnitude of components of selection on female EPR in the wild, and consequent evolutionary predictions, therefore remain scant [5,6,8,28].
Estimating additive genetic covariances between female EPR and relative fitness components in wild populations requires appropriate traits to be measured in numerous relatives. I applied 'animal models' to 18 years of comprehensive life-history, paternity and pedigree data from song sparrows (Melospiza melodia) to estimate cov A between female EPR defined as a female's liability (hereafter 'EPR L ') to produce extra-pair offspring (EPO, sired by an extra-pair male) rather than within-pair offspring (WPO, sired by her socially paired male [29]), and two specific fitness components: adult female annual reproductive success (ARS) and survival to recruitment (SR). These covariances, respectively, predict per-generation evolutionary responses to (correlated) selection on female EPR L associated with female reproductive success, including failure, and pre-reproductive mortality. Evolutionary responses attributable to these two components of selection are of specific interest because they are potentially substantial but cannot be adequately predicted through solely phenotypic analyses. I additionally considered the degree to which evolutionary change in EPR L due to selection through ARS predicted by basic application of the univariate breeder's equation would be biased. I thereby demonstrate a quantitative genetic approach to predicting per-generation evolutionary responses to selection on female liability for EPR as encapsulated by additive genetic covariances with two specific fitness components. Key nomenclature is summarized and defined in table 1.

MATERIAL AND METHODS
(a) Study system Song sparrows of both sexes typically breed two to three times per season starting from age 1 year and are primarily socially monogamous with clear social pairings. Females incubate clutches (typically three to four eggs) and both socially paired parents provision hatched offspring. However, song sparrows are genetically polygynandrous, with frequent EPR [30,31].
Mandarte Island, British Columbia, Canada, (approx. six hectares) holds a resident song sparrow population, recently numbering roughly 15-45 breeding pairs, which has been studied intensively since 1975 [32]. Each year, all breeding attempts are closely monitored and all nests located. All offspring surviving to approximately 6 days post-hatch are marked with unique combinations of coloured rings to allow individual identification. The occasional immigrants to Mandarte (1.1 per year on average) are colour-ringed soon after arrival. All social parents of all offspring (those incubating clutches or provisioning chicks) are identified by observation. Immigration is sufficient to maintain neutral allelic diversity [33] and prevent inbreeding from accumulating.
During 1993-2009, 99.4 per cent of ringed offspring and adults were blood sampled and genotyped at 13 polymorphic microsatellite loci to allow assignment of genetic parents [30]. Offspring were sexed using standard molecular methods quantity definition EPR L a female's liability to produce an extra-pair offspring (EPO) rather than a within-pair offspring (WPO) ARS w a female's relative ARS defined as the number of ringed offspring a female produced per year divided by the population mean SR w an individual's relative SR defined as 1 or 0 for individuals that did or did not survive to age 1 year divided by the population mean; SR w was analysed as a binomial trait (SR wB ) and as a Gaussian trait (SR wG ) V A , cov A additive genetic variance and covariance, respectively V Y , cov Y year variance and covariance, respectively V PI , cov PI permanent individual variance and covariance, respectively V R , cov R residual variance and covariance, respectively [34]. Bayesian full probability models that incorporated genetic and spatial information assigned genetic sires to 99.6 per cent of sampled offspring with at least 95 per cent individual-level confidence. Overall, approximately 29 per cent were assigned to males other than a female's socially paired mate and hence were EPO [30] (compared with 24% in a nearby mainland song sparrow population [31]). The probability of excluding a female's social mate as sire averaged 0.9998, and the number of unsampled sires was effectively zero. All genetic mothers matched those assigned by behaviour [30]. Each adult (at least 1 year old) female's ARS was measured as the number of offspring ringed (approx. 6 days post-hatch) per season. Owing to Mandarte's small size and the intensive fieldwork, the probability of resighting an adult during any breeding season is effectively one. Whether or not each ringed offspring survived to recruit at age one (hereafter 'SR') was therefore documented with high confidence [35]. The local recruitment of all ringed offspring and ARS of all adult females in the population were therefore completely described. The high local recruitment rate and scarcity of ringed sparrows on surrounding islands suggest that emigration is relatively rare [32,34]. Since the 'secondary theorem of selection' specifies relative rather than absolute fitness [9,20,26], relative ARS and SR (ARS w and SR w , table 1) were calculated by dividing individual values by the overall across-year population mean with binary SR coded as 1/0. Estimated additive genetic (co)variances were similar when ARS w and SR w were calculated relative to year-specific rather than across-year means.
(b) Quantitative genetic analyses Multivariate 'animal models' were used to estimate additive genetic variances (V A ) and covariances (cov A ) in and among female EPR L , ARS w and SR w . Animal models are mixed models in which pairwise coefficients of kinship (k) estimated from pedigree data define a matrix proportional to the variance-covariance structure of additive genetic random effects [26,36]. Such models estimate genetic parameters for a baseline population that, in practice, comprises individuals in the pedigree that have unknown parents [36] (see the electronic supplementary material).
Focal traits and fitness components may not be normally distributed, and relative fitness follows no obvious statistical distribution [9,12,22]. While deviation from normality may not bias estimates of additive genetic (co)variances, it does impede hypothesis testing [9,12,17,22]. One solution is to consider traits that take discrete values as 'threshold traits' where normally distributed continuous variation in underlying 'liability' is assumed to translate into trait expression at certain threshold values [26,37,38]. Statistically appropriate models can then be fitted and interpreted on underlying liability scales. While considering liabilities has statistical advantages, cov A between a trait and liability for relative fitness has no clear quantitative interpretation in terms of predicted evolutionary change. I therefore fitted multiple models (described below) on appropriate observed and liability scales to provide evolutionarily meaningful estimates of cov A and associated statistical confidence.
Female EPR L , defined as a female's liability to produce EPO rather than WPO, was modelled as a binomial threshold trait with EPO and ARS as numerator and denominator, respectively [29]. Females whose breeding failed (ARS ¼ 0) are consequently uninformative in univariate analyses of EPR L . EPR L , defined as a liability rather than simply the observed ratio of EPO to ARS, is not necessarily correlated with ARS at either phenotypic or genetic levels. ARS w was best approximated by a Gaussian distribution (see §3). SR w was modelled as a binary threshold trait (SR wB ) to estimate cov A with EPR L with interpretable statistical confidence, and as a Gaussian trait (SR wG ) to provide a quantitatively interpretable evolutionary prediction.
Additive genetic variance in traits and fitness components may be small, especially relative to phenotypic variance. This situation is directly informative since cov A (and hence predicted evolutionary responses to selection) must be zero if V A ¼ 0 [16,20]. However, small V A can impede precise estimation of cov A , especially given the limited sample sizes typically available from wild populations [9,15,22]. One approach is to estimate V A in each trait and fitness component of interest and then estimate cov A among traits and components where V A . 0 [15]. In contrast, I fitted pairwise bivariate models between EPR L , ARS w and SR w in order to maximize efficient use of available information to estimate V A in case cov A = 0, and hence maximize statistical power [14].
Models included random effects of breeding year for EPR L and ARS w , or natal year for SR w , and hence estimated among-year variances (V Y , table 1) [29,35]. Models for female EPR L and ARS w included random individual effects to account for repeat observations of females that survived multiple years and estimate 'permanent individual' variance (V PI , which comprises permanent environmental and nonadditive genetic variances [36]). These models also included fixed effects of three age categories defined through preliminary analysis (ages 1, 2 -3 and greater than 3 years). All models included fixed regressions on individual coefficient of inbreeding ( f ) because unmodelled inbreeding depression may inflate estimates of V A , cov A and phenotypic selection [39,40]. The inter-sex genetic correlation for SR is close to one in the study population [35]. Models for SR w therefore used recruitment data from all ringed offspring to maximize statistical power, and included fixed effects of sex because mean SR differs between males and females [34]. Effects were similar when estimated using SR w data for female offspring only. Initial models included random effects of an individual's mother, social father or breeding territory, but estimated variances were close to zero and other parameter estimates were scarcely affected when these terms were removed. Random natal brood effects were not modelled because few same-brood females recruited, and there is little among-brood variance in SR [35].
Model covariance structures were unconstrained, allowing estimation of additive genetic (cov A ) and year (cov Y ) covariances among all traits, and of permanent individual (cov PI ) and residual (cov R ) covariances between female EPR L and ARS w . Since residual variance (V R ) cannot be estimated and V PI is not identifiable for single-observation binary traits (such as SR wB ), V R and cov R for SR wB were fixed to one (by convention) and zero, respectively, and V PI and cov PI were not estimated. Cov R between SR wG and EPR L was also fixed to zero because there is zero phenotypic variation in SR w across females whose EPR was observed.
(c) Analysis implementation Analyses used SR w data for all offspring ringed during 1993-2009 (when paternity was verified), and EPR L and ARS w for all Mandarte-hatched adult females alive during 1993 -2009 (including cases of ARS w ¼ 0). Pedigree data for females and offspring ringed during 1993-2009 were compiled from genetic parentage data taking the most likely sires [29,30,41].
Pedigree data for females and ancestors ringed during 1975-1992 were compiled from observed social parentage and combined with the 1993-2009 genetic data to provide a full pedigree for 1975-2009. Assuming the unobserved EPR rate during 1975-1992 was similar to the approximately 29 per cent observed during 1993-2009 and that all mothers were correctly assigned by social behaviour, approximately 86 per cent of all 1975-1992 pedigree links will be correct. Paternity error for breeders hatched before 1993 will actually be less than this if, as during 1993-2009, EPO were less likely to recruit than WPO [34,35]. The 1975-1992 pedigree data, even though uncorrected for EPR, therefore provide information regarding k between the 1993 breeders and allow the alternative assumption of k ¼ 0 to be relaxed [41]. Since the contribution of a common ancestor to k between two focal individuals decreases by a factor of two per generation of separation, paternity error in the early pedigree introduces little error into estimates of k among recent breeders. Using the full 1975-2009 pedigree data therefore allows full use of all available phenotypic data and provides the most powerful analysis feasible (see [29,41]). Animal model pedigrees were pruned to individuals with phenotypic data and all their known ancestors, and thereby restricted to individuals that are informative for current analyses (see the electronic supplementary material). Individual f-values were calculated from the 1975-2009 pedigree using standard algorithms, and therefore measure the probability of identity by descent relative to pedigree founders [42]. Kinship between immigrants and existing natives, and hence f of offspring of immigrantnative pairings, was defined as zero [42,43]. Phenotypic data from 11 immigrant females were excluded because f is undefined for immigrants (as opposed to their offspring).
Models were fitted using Bayesian methods implemented in MCMCglmm 2.14 in R v. 2.12.2 [44,45] with logit link functions for threshold traits. Fixed effect priors were normally distributed and diffuse with mean zero and large variance (10 8 ). Prior variances were inverse-Wishart distributed with limit variance of one and low degree of belief (0.002). Prior covariances were zero with low degree of belief (0.002). Priors were therefore only weakly informative, and posterior distributions were robust to reasonable prior variation, including covariances of 20.5 to 0.5. Analyses used at least 3 005 000 iterations, burn-in 5000 and thinning interval 3000 to ensure low (less than 0.05) autocorrelation among thinned samples. Posterior means and 95 per cent credible intervals (95% CI) for fixed effects, (co)variances and heritabilities (h 2 ) were estimated across thinned samples. Lower limits on estimates are bounded to zero for variances but not covariances. Posterior distributions for variances were therefore further inspected to assess modes and density close to zero, and hence evaluate the evidence that variances differed from zero. Simulations of null traits with V A ¼ 0 further substantiated these inferences (see the electronic supplementary material). Liability-scale heritabilities of EPR L and SR wB were estimated from binomial models as given logistic variance proportional to p 2 /3 [46]. Heritabilities of ARS w and SR wG were estimated as h 2 To aid interpretation, the per-generation magnitude of evolutionary change in EPR attributable to selection through ARS w or RS w predicted on the liability scale was inverse-logit back-transformed to the observed phenotypic scale.
Full trivariate models for EPR L , ARS w and SR w could not be properly fitted in MCMCglmm because the correct residual covariance structure (with estimated cov R between EPR L and ARS w but zero cov R with SR w ) cannot be specified. A further bivariate model was therefore fitted to estimate cov A between ARS w and SR wB , and hence consider whether estimates of cov A between EPR L and ARS w or SR wB could be biased by unmodelled cov A with SR wB or ARS w , respectively. Trivariate models fitted as rigorously as feasible supported conclusions drawn from bivariate models (see the electronic supplementary material). Data are available at the Dryad Repository: doi:10.5061/dryad.907cv.
(d) Phenotypic selection differential My primary aim was to estimate cov A between EPR L and two specific components of relative fitness (ARS w and SR w ). However, I also compared resulting evolutionary predictions with basic implementation of the univariate breeder's equation, R ¼ h 2 S, where the phenotypic selection differential (S) equals the phenotypic covariance (cov P ) between trait and relative fitness, S ¼ cov P (z,w) [9,15]. Cov P between EPR L and ARS w was estimated as the residual covariance from a bivariate mixed model (assuming binomial and Gaussian error distributions, respectively). This model included random year effects and fixed effects of age and f (as in the animal model used to estimate cov A ), but estimates remained similar when these fixed effects were excluded. Cov P was therefore estimated across female-years where ARS w exceeded zero, and hence from a subset of the data used to estimate cov A (see §3). Cov P between EPR L and SR w is zero because EPR is only observable for females that recruited.

RESULTS
(a) Female extra-pair reproduction and annual reproductive success A total of 224 Mandarte-hatched adult female song sparrows were alive during 1993 -2009, totalling 474 female-years (see the electronic supplementary material). Female ARS varied from zero to 11 ringed offspring (mean 4.7, variance 5.5, median 5, inter-quartile range 3-6, skew 20.1; figure 1). The distribution of ARS w differed from normality (Kolmogorov-Smirnov p ¼ 0.005) but was only slightly skewed. Overall, 211 of the 224 females had at least one offspring in at least one year during 1993 -2009, totalling 446 female-years. Across these cases, the observed proportion of ringed offspring that were EPO ranged from zero to one (mean 0.29, variance 0.10, median 0.20, skew 0.83; figure 1). This proportion did not increase or decrease during 1993 -2009 (correlation coefficient r ¼ 0.001).
(b) (Co)variances in female EPR L and ARS w The bivariate animal model for female EPR L and ARS w estimated non-zero V A and h 2 in both traits with low Genetic covariances and polyandry J. M. Reid 4655 posterior density close to zero, and significant inbreeding depression in ARS w but not EPR L (table 2). Estimated cov A between EPR L and ARS w was positive but small, and the 95% CI overlapped zero (table 2). Estimates of V PI in EPR L and ARS w were small with relatively high posterior density towards zero, and cov PI was small with a 95% CI that overlapped zero (table 2). There was substantial V Y in ARS w but not EPR L , and consequently small cov Y with a 95% CI that overlapped zero (table 2). There was substantial V R in both EPR L and ARS w , and negative cov R with a 95% CI that did not overlap zero (table 2).
Across 446 observations of 211 females with non-zero ARS, the phenotypic covariance (cov P ) between EPR L and ARS w was -0.026 (95% CI -0.117 -0.075). (d) (Co)variances in female EPR L , ARS w and SR w The bivariate animal model for female EPR L and SR wB estimated non-zero V A and h 2 in both traits with low posterior density close to zero, and substantial inbreeding depression in SR wB , but not EPR L (table 2). Estimated cov A between EPR L and SR wB was positive, but the 95% CI was wide and overlapped zero (table 2). There was substantial V Y in SR wB , but small cov Y with a 95% CI that overlapped zero (table 2). Estimated cov A between EPR L and SR wG was also positive but small (table 2). The bivariate animal model for ARS w and SR wB estimated that cov A between these fitness components was small with a narrow 95% CI that overlapped zero (table 2).

DISCUSSION
Despite huge interest in identifying forces driving the evolution of polyandry and female EPR [1,2,4,8], selection differentials or gradients on EPR and corresponding evolutionary responses have not been explicitly estimated or predicted in wild populations [5,6]. Direct costs and benefits of EPR have instead been postulated in terms of decreased paternal care by a female's cuckolded social mate, time, energy or disease costs of mating, or increased fertility or foraging opportunities [1,5,8,47]. However, estimating such effects is difficult because female EPR cannot be readily manipulated in isolation from other behavioural or life-history traits of a female and/or her social mate [47]. Without such experiments, estimated costs and benefits may be biased by environmental covariances between EPR and behaviour, physiology or life history [8,27,28], and may inadequately represent fitness anyway.
Direct selection on polyandry has been experimentally demonstrated in laboratory populations of facultatively polyandrous taxa, particularly where males provide nuptial resources or sperm is limiting, and female fecundity consequently increases with multiple mating and polyandry per se [8] (but see [2]). In contrast, the hypothesis that EPR by socially monogamous females is under positive direct selection through increased female fecundity or reproductive success is often rejected because there is no obvious mechanism rather than through rigorous empirical test. Female EPR is consequently hypothesized to be under positive indirect selection through increased offspring fitness [1,4,28]. However, tests of this hypothesis, which typically quantify whether EPO are phenotypically fitter than their maternal half-sib WPO, have themselves failed to demonstrate any consistent positive effect and may be biased anyway [5,8,34,35,48]. Furthermore, phenotypic selection analyses cannot predict evolutionary responses to correlated selection on female EPR through pre-reproductive mortality or reproductive failure unless all correlated traits are identified and measured [14,24,25]. Purely phenotypic approaches may consequently be inadequate to accurately predict evolutionary responses to (correlated) selection on EPR and polyandry through major fitness components, such as recruitment and reproductive success.
To attempt to circumvent these problems, I estimated additive genetic covariances (cov A ) between female liability to produce EPO (EPR L ) and relative annual reproductive success (ARS w ) and survival to recruitment (SR w ) in song sparrows, thereby parametrizing the 'secondary theorem of selection' with respect to these two specific fitness components. I thereby demonstrate an approach to predicting per-generation evolutionary responses associated with specific components of (correlated) selection on female EPR L that should be unbiased by environmental covariances and unmeasured traits.
(a) Genetic (co)variances There was moderate additive genetic variance (V A ) and heritability (h 2 ) in female EPR L in song sparrows (see also [29]), and small but non-zero V A and h 2 in female ARS w . There was therefore potential for non-zero cov A between EPR L and ARS w , and hence for a positive or negative evolutionary response to selection on EPR L through female ARS. In fact, estimated cov A between EPR L and ARS w was positive, predicting evolution towards increased female EPR L . The posterior mean estimate of approximately 0.09 (table 2) equates to a predicted increase in the proportion of offspring that are EPO of approximately 0.02 per generation on the back-transformed observed scale. However, although the 95% CI was reasonably narrow (20.03 to 0.20), reflecting the completeness of the song sparrow pedigree and the relatively high relatedness among observed females, it overlapped zero. There was therefore no conclusive prediction of an evolutionary increase in female EPR L associated with selection through female ARS. This analysis incorporates correlated selection through reproductive failure (ARS ¼ 0), which cannot contribute to direct selection on EPR L defined as a female's liability to produce EPO rather than WPO or hence be detected by univariate phenotypic selection analysis.
There was also significant V A and h 2 in SR in song sparrows (see also [35]) and hence potential for non-zero cov A between EPR L and SR w . Such covariance, either positive or negative, can be interpreted as a predicted evolutionary response to correlated selection on EPR L through prerecruitment mortality, and hence through individuals that die before EPR L can be expressed [14,24]. Any such selection is by definition indirect, for example reflecting genetic variation underlying physiology or behaviours that affect Genetic covariances and polyandry J. M. Reid 4657 both recruitment and EPR L . Estimated cov A between EPR L and SR w was positive, implying evolution towards increased female EPR L through correlated viability selection. However, the 95% CI for SR wB was wide (20.19 to 0.52) despite the relatively large number of observed individuals and overlapped zero. This reflects reduced power to estimate cov A among traits expressed across rather than within individuals, and intrinsic sampling variance associated with estimating liabilities underlying binary traits such as recruitment. There was therefore no conclusive prediction of evolutionary change in female EPR L owing to correlated selection through pre-reproductive mortality. Estimated cov A between ARS w and SR wB was small with a narrow 95% CI that overlapped zero (20.05 to 0.09). There was therefore no evidence of a substantial genetic trade-off between liability to recruit and female ARS as might be most simply manifested as negative cov A [49] (but see [50]). Estimated cov A values between EPR L and ARS w and SR w estimated from bivariate models are therefore unlikely to be biased by unmodelled cov A with SR w and ARS w , respectively. Indeed, estimates of cov A from trivariate models were qualitatively similar to those from pairwise bivariate models and were therefore robust across the set of three traits considered (see the electronic supplementary material).
Polyandrous female song sparrows were previously shown to produce EPO of lower phenotypic and additive genetic value for fitness components than their WPO [34,35]. There may therefore be weak negative selection against female EPR L owing to low additive genetic value of EPO [35]. In contrast, estimates of cov A between EPR L and ARS w and SR w were positive (although not significantly different from zero). Since combined selection on any trait equals the sum of selection through multiplicative components of relative fitness [10], the response to selection on EPR L associated with SR and subsequent ARS can be predicted as the sum of cov A with SR wG and ARS w , giving roughly 0.18 (table 2). This equates to a predicted evolutionary increase in the proportion of offspring that are EPO of approximately 0.03 per generation on the back-transformed observed scale associated with these two specific fitness components. Combined positive selection on female EPR L associated with ARS and SR may therefore partially balance negative selection through low additive genetic value of EPO and help maintain EPR in this population. However, an appropriate 95% CI for the summed cov A is hard to estimate because SR w is not normally distributed, and there is no clear quantitative evolutionary interpretation of cov A between EPR L and liability for relative recruitment.
(b) Environmental covariances Joint estimation of additive genetic and environmental covariances between phenotypic traits and relative fitness components, as achieved by fitting 'animal models', allows evolutionary responses to selection to be predicted independent of environmental covariances that can confound breeder's equation predictions [15]. Estimating both these distinct covariances provides biological insight and allows key assumptions underlying breeder's equation predictions to be validated [15,20,21].
In song sparrows, the 'permanent individual' covariance (cov PI ) between female EPR L and ARS w , which comprises permanent environmental and non-additive genetic effects that consistently influence individual values of both traits, was positive but small. However, the residual covariance (cov R ), which comprises seasonspecific individual effects on both traits, was significantly negative, and hence opposite in sign to the posterior mean cov A . Negative cov R implies that females with high ARS w due to season-specific individual environmental effects also have low EPR L . This could, for example, reflect correlated environmental effects on female ARS w and male mate guarding (if both vary with resource availability), or negative effects of female EPR on paternal care and hence ARS w , or other similar mechanisms [47].
The estimated phenotypic covariance (cov P ) between female EPR L and ARS w was negative (although not significantly different from zero), reflecting the negative cov R . The evolutionary response to selection on EPR L predicted by simple application of the univariate breeder's equation (R ¼ h 2 S % 0.17 Â 20.03 % 20.005) would consequently be opposite in sign to the posterior mean cov A (contra [21]). The long-standing concern that phenotypic associations between EPR and fitness components may partly reflect environmental covariances is therefore justified [27,28]. However, the 95% CI for cov A between ARS w and EPR L of 20.027 to 0.199 includes the point estimate of R from the univariate breeder's equation. Bias in predicted evolutionary change was consequently relatively small in the current case.

(c) Context and interpretation
Estimating selection on any trait through any fitness component(s), and predicting any consequent evolutionary response, is extremely challenging in wild populations [11,13,21]. Any association-based estimates, whether phenotypic or quantitative genetic, should be interpreted cautiously, and ideally confirmed experimentally [12,16,20,25]. However, for traits such as EPR that cannot be readily experimentally manipulated in the wild, quantitative genetic approaches have advantages over purely phenotypic analyses; namely that evolutionary responses to (correlated) selection can be predicted without bias owing to environmental covariances or unmeasured traits [14,15,20,24]. The current analyses demonstrate this approach with respect to EPR L and two specific fitness components, but some provisos and challenges of application and interpretation remain.
Evolutionary change in female EPR L will ultimately depend on cov A with all components of female and male lifetime fitness from conception [11], including relative adult survival, male ARS (given that cross-sex genetic correlations for fitness components may be less than one [51]) and survival to hatch as well as ARS w and SR w . Indirect genetic effects stemming from male-female interactions, correlated selection through associated mating behaviours (such as extra-pair copulations rather than EPR per se) and non-additive genetic effects could also contribute [6,20,52,53]. Substantial further theory, data and analyses are required to estimate and combine all possible components of selection on EPR L , and achieve the ultimate goal of predicting and understanding the evolution of polyandry. The most tractable and insightful approach to this task is initially to quantify and understand each selection component individually. Accordingly, my current aim was to estimate cov A between EPR L and two specific fitness components (ARS w and SR w ), and hence predict pergeneration evolutionary responses to selection associated with these two specific components. Failure to model further genetically correlated traits and fitness components does not invalidate conclusions regarding the two focal components, but does impede clear distinction between predicted responses to direct and indirect (correlated) selection [20], and means that overall evolutionary responses cannot yet be predicted.
Furthermore, the basic secondary theorem of selection predicts the evolutionary response to directional selection over one generation. In general, 'animal model' analyses that use phenotypic data from multiple (overlapping) generations provide unbiased estimates of key genetic parameters in the baseline population [26,36] (see the electronic supplementary material). Such analyses therefore do not necessarily measure or explain contemporary evolutionary or phenotypic dynamics of EPR L , but rather provide point predictions of evolutionary responses associated with specific components of selection. Moreover, the secondary theorem of selection predicts the evolutionary response to directional selection in terms of a change in mean breeding value for the focal trait (DA z ). This approach can be extended to estimate the response to stabilizing or disruptive selection, defined as a change in additive genetic variance (DV Az ), as cov A between relative fitness and the squared deviation of the trait from its phenotypic mean (DV Az ¼ cov A (w,(z 2 m z ) 2 ), where m z is mean z [20,23]). However, this covariance cannot be readily estimated for liabilities underlying threshold traits (such as EPR L ), meaning that potential responses to nonlinear selection are not easily predictable.
Low power is likely to afflict most field estimates of genetic covariances, particularly those involving fitness components [16,21,22]. The song sparrow dataset provides reasonable power for such analyses despite the relatively small number of observed individuals; the comprehensive pedigree data mean that relatedness between these individuals is estimated relatively accurately and is typically non-zero. Genetic covariances of roughly 0.15-0.25 can consequently be detected [35]. However, even given this dataset, the 95% CIs for the estimated cov A between EPR L and SR wB were wide, indicating low power stemming from high sampling variance associated with estimating liabilities underlying threshold traits. The most robust general conclusions regarding components of selection on EPR, and other traits, may ultimately come through meta-analysis [2,8], if sufficient unbiased estimates of key parameters can be provided. Cov A between EPR L and relative fitness components should therefore be estimated as precisely as feasible in a range of populations. This should become increasingly tractable as multi-generational paternity assignment studies that encompass large sample sizes become available.
I thank everyone who contributed to data collection on Mandarte, Rebecca Sardell for her careful paternity assignment, Peter Arcese and Lukas Keller for their support, Michael Morrissey, Erik Postma and three anonymous reviewers for their thoughtful comments, and the Royal Society for funding.