IT TAKES TWO: HERITABLE MALE EFFECTS ON REPRODUCTIVE TIMING BUT NOT CLUTCH SIZE IN A WILD BIRD POPULATION

Within-population variation in the traits underpinning reproductive output has long been of central interest to biologists. Since they are strongly linked to lifetime reproductive success, these traits are expected to be subject to strong selection and, if heritable, to evolve. Despite the formation of durable pair bonds in many animal taxa, reproductive traits are often regarded as female-specific, and estimates of quantitative genetic variation seldom consider a potential role for heritable male effects. Yet reliable estimates of such social genetic effects are important since they influence the amount of heritable variation available to selection. Based on a 52-year study of a nestbox-breeding great tit ( Parus major ) population, we apply ‘extended’ bivariate animal models in which the heritable effects of both sexes are modelled to assess the extent to which males contribute to heritable variation in seasonal reproductive timing (egg laying date) and clutch size, while accommodating the covariance between the two traits. Our analyses show that reproductive timing is a jointly expressed trait in this species, with (positively covarying) heritable variation for laydate being expressed in both members of a breeding pair, such that the total heritable variance is 50% larger than estimated by traditional models. This result was robust to explicit consideration of a potential male-biased environmental confound arising through sexually dimorphic dispersal. In contrast to laydate, males’ contribution to heritable variation in clutch size was limited. Our study thus highlights the contrasting extent of social determination for two major components of annual reproductive success, and emphasises the need to consider the social context of what are often considered individual-level traits.

with conspecifics. Heterogeneity in the social aspect of the environment likely represents a source of phenotypic variance for many traits. If variation in the influential behavior of associates has a heritable component, then indirect genetic effects (IGEs; also known as associative genetic effects (Muir 2005) or social genetic effects (Ellen et al. 2014)) are present, which means that individuals' phenotypes are influenced by the genes of others. IGEs mediated by conspecifics can greatly reduce the total heritability of a trait if they covary negatively with the direct genetic effects (Dickerson 1947), such that the population's response to selection is slowed or even reversed (Moore et al. 1997;Muir 2005;Bijma and Wade 2008). Conversely, a positive covariance will inflate the effective heritability, and is thus expected to accelerate the response to selection, all else being equal. Detailing how social partners influence each other's phenotype and elucidating the evolutionary impact of the genome being shared among social partners is thus essential if we are to improve our understanding of the evolutionary dynamics operating in nature.
The study of IGEs in wild populations has been dominated by a focus on maternal genetic effects (McAdam et al. 2014). However, social interactions between non-relatives can also give rise to IGEs (Wade et al. 2010;Costa e Silva et al. 2013), although published estimates for populations in natural settings are scant (Hunt et al. 2019). In this context, sexual partnerships have clear potential to facilitate indirect genetic effects, particularly in species where pairs form an association that persists beyond fertilization (Carter et al. 2019;Chakrabarty et al. 2019). In many animal taxa, for example, both sexes exhibit parental care, such that raising offspring is a shared task. Given this interdependence of reproductive success, sexual partnerships are settings in which we might expect to observe IGEs via non-relatives, and relevant data are often collected as part of long-term population studies. Yet for many study populations, observations of marked individuals across multiple breeding events indicate that expression of reproductive behaviors shows consistent variation among females but not among males (van Noordwijk et al. 1980(van Noordwijk et al. , 1981aSheldon et al. 2003;Gienapp et al. 2006;Browne et al. 2007). However, a relative shortage of repeat measures for short-lived species may limit quantification of individual repeatability (Browne et al. 2007), and the smaller effect sizes expected for males would be more strongly affected by statistical power limitations. Consistent with this, studies of longer-lived species report that male heterogeneity contributes to observed variation in breeding behavior, as for the sparrowhawk (Accipiter nisus: Newton and Marquiss 1984), tawny owl (Strix aluco: Brommer et al. 2015), and mute swan (Cygnus olor: Charmantier et al. 2006).
A quantitative genetic approach-and, in particular, the animal model (Henderson 1953(Henderson , 1975-is able to infer the heritable component of individual heterogeneity by leveraging the (partial) genotypic replication inherent in sampling known relatives. Indeed, in another passerine bird, the song sparrow (Melospiza melodia), a heritable contribution of male heterogeneity to variation in reproductive timing has been described, and a positive covariance of the male and female genetic effects means the estimated total heritable variance in laydate is doubled once the contribution via males is explicitly considered (Germain et al. 2016). However, while the pair-level organization of reproductive behaviors makes them well-suited for quantification of IGEs (Bijma 2010) and their direct relevance to individual fitness means that their evolutionary dynamics are of particular interest to biologists, the findings of the small number of published studies relying on animal modeling are inconsistent (Brommer and Rattiste 2008;Caro et al. 2009;Teplitsky et al. 2010;Liedvogel et al. 2012;Germain et al. 2016). For example, in contrast to Germain et al.'s (2016) study of song sparrows, Brommer and Rattiste (2008) report a negative genetic covariance between female and male effects on laydate in common gulls (Larus canus), yielding a lower total heritability estimate than when additive genetic variance via males is ignored. Furthermore, whether male genetic effects extend to reproductive quantum (i.e., clutch or litter size), the other principal axis of within-population variation in reproductive output, has been more neglected still (though see: van Noordwijk et al. 1980;Gibbs 1988;van der Jeugd and McCleery 2002;van Noordwijk et al. 1981a), despite being linked to both reproductive timing and lifetime reproductive success in many avian study populations (Klomp 1970).
The animal model, which is a linear mixed-effects model in which a population pedigree is linked to an individual identity matrix, is well suited to studying the quantitative genetics of wild populations (Kruuk 2004). Besides inferring the contribution of additive genetic effects to phenotypic (co)variance, the animal model can accommodate explicit consideration of particular environmental contributions that might otherwise bias additive genetic statistics. For example, the shared natal environment of siblings contributes to the phenotypic resemblance of relatives and can upwardly bias additive genetic (co)variance estimates (Kruuk and Hadfield 2007;Hadfield et al. 2013). Yet relatives can continue to share environments beyond independence through a combination of natal philopatry (Regan et al. 2017) and local-scale environmental variation (Wilkin et al. 2007a). If dispersal distances differ between the sexes, as is typical of the mammalian and avian populations used in quantitative genetic research of wild populations (Greenwood 1980), then joint consideration of male and female contributions to heritable variation may be subject to sex-dependent environmental bias. In the case of birds, females generally disperse further from their natal origin (Greenwood 1980), such that environmentally driven spatial heterogeneity in the expression of the focal trait is expected to exert a greater bias on the estimation of male genetic effects. However, the phenotypic variance attributable to such positive spatial autocorrelation can itself be modeled (Costa e Silva et al. 2001;Dutkowski et al. 2002), as has been applied in quantitative genetic studies of free-living animal populations (Stopher et al. 2012;Germain et al. 2016;Regan et al. 2017).
We use data from a study of nestbox-breeding great tits (Parus major) that spans more than half a century to quantify the heritable contributions of females and males to phenotypic variation in the two major components of the reproductive behavioral syndrome: laydate and clutch size. While egg laying is performed by the female, male great tits play an important role in the reproductive process by maintaining a breeding territory that excludes conspecifics, the quality of which is correlated with laydate and clutch size (Wilkin et al. 2007a). Heterogeneity among males is known to influence the quality of territory they hold (Krebs 1971), so it seems reasonable to suggest that male heterogeneity may influence reproductive behavioral outcomes. If part of this male heterogeneity has a heritable component then the potential for male genetic effects will be realized. Fatherson regressions have been used to estimate male heritability of both laydate and clutch size in great tits but statistical uncertainty was large (van Noordwijk et al. 1980(van Noordwijk et al. , 1981avan der Jeugd and McCleery 2002). Animal modeling offers advantages over parent-offspring regression in this respect, since it leverages the known relatedness structure of the population in full. Additionally, sources of bias arising through the common environments of relatives (van der Jeugd and McCleery 2002) can be explicitly accommodated. Here, we constructed bivariate animal models of laydate and clutch size to accommodate their covariation, comparing a 'classic' animal model structure that ignores male contributions to expression of these traits with an 'extended' animal model in which the effects of heritable male heterogeneity are explicitly quantified. We incorporated quantification of spatial autocorrelation into a third model to assess the extent to which those heritable male effect estimates are prone to inflation as a result of the more limited natal dispersal of male great tits (Greenwood et al. 1979;van der Jeugd and McCleery 2002).

Methods
We used long-term records of great tits (Parus major) breeding in Wytham Woods, near Oxford (UK), in which 1020 nestboxes are in place across the 385 hectares of mixed deciduous woodland. The number and position of nestboxes have been fixed since 1965 so we used breeding records from then until 2016, giving us a dataset spanning 52 years. Previous analyses indicate that a large majority of the breeding population of great tits adopts nestboxes as nesting sites (Greenwood et al. 1979). As described elsewhere (Perrins 1965;McCleery and Perrins 1989), nestboxes were visited regularly from early April onward to check for signs of nest-building and to count eggs. Laydate is defined as the day on which the first egg was laid, assuming that females lay one egg per day without interruption, with each egg laid early in the morning (i.e., prior to any human observation on that day). Thus, if eggs are observed in a nestbox prior to completion of the clutch, laydate is inferred by back-counting. Laydate is recorded on a calendar-based scale in which the 1st April is defined as day 1. Clutch size represents the total number of eggs observed in the nest. To overcome convergence difficulties, we excluded data with values for either trait that were beyond the 0.5th and 99.5th percentiles (representing 1.5% of data), which reduced the range of laydates from 0-76 to 5-55 and the observed range of clutch sizes from 1-21 to 4-13 ( Fig. 1). Nests subjected to experimental manipulations were also excluded, as were known repeat and second clutches, nests for which the identity of either provisioning adults was unknown, and nests for which laydate or clutch size was unknown. Provisioning adults are presumed to be the biological parents of the chicks in the nest, and are caught at their nestbox during chick-feeding and identified using uniquely numbered, metal leg-rings. Extra-pair paternity will cause errors in the population pedigree but the proportion of chicks sired by extra-pair males in this population is of such a magnitude (13%: Patrick et al. 2012) that its impact on estimates of additive genetic (co)variance is expected to be minor (Firth et al. 2015a;Charmantier and Réale 2005;Perrier et al. 2018). In relying on a social (i.e., non-genetic) pedigree, we assume that extra-pair sires exert a negligible influence on the laydate and clutch size of the cuckolding female: the strong breeding territoriality of this species (Krebs 1971) likely means that an extra-pair sire's interaction with the resident pair is very limited. In the pedigree, we assigned a unique 'dummy identity' to unidentified parents (i.e., those that were not caught at the nest) to ensure sibships among their offspring were recognised in the pedigree.
As for many avian populations (Klomp 1970), laydate and clutch size covary in this population . We, therefore, adopted a multivariate animal modeling approach to quantify the contributions of additive genetic and environmental effects to phenotypic (co)variation in these two principal axes of reproductive behavioral variation, which involved three successively more complex bivariate animal models. The first of these followed a "classic" animal model structure, in which only contributions to phenotypic (co)variance via females were explicitly considered, thus ignoring the potential for heritable or environmental effects mediated by males. Alongside the additive genetic effect, "year" and "plot" (the study site is spatially divided into nine nestbox plots, ranging in size from 19 to 121 hectares, which originally served to define areas of responsibility for individual fieldworkers but which also approximate coarse-scale environmental heterogeneity within the site) were also fitted as random effects to model these sources of temporal and spatial variance. Our initial models included maternal identities of both females and males as random effects, to quantify environmentally derived phenotypic resemblance of full siblings and maternal half-siblings ("brood identity" contained insufficient within-subject repeats to be informative). However, these did not converge and univariate models revealed natal environmental effects to be minor (Tables S1 and S2), so maternal effects were not fitted in the presented models. The first model decomposes the (conditional) phenotypic variance-covariance matrix of the two traits, P, into environmental matrices (year, Y; plot, Plot; female permanent environmental, ♀PE; and residual, R, effects) and an additive genetic matrix, G ♀ , that describes the contribution of heritable variation among females: Selection on these traits is typically estimated on a yearspecific basis (van Noordwijk et al. 1995;Charmantier et al. 2008;Visser et al. 2015), so we used the estimate of within-year phenotypic variance (the sum of all variance effects (conditioned upon the fixed effect(s) of individual age: Wilson 2008) except year) in estimating the corresponding narrow-sense heritability (i.e., h 2 = V ♀A VP (within−year) ). Additionally, the long-term shift in reproductive behavior in our study population (McCleery and Perrins 1998) would likely elevate our estimate of annual variance beyond what individuals would experience (the oldest individuals in our dataset were last observed breeding at 9 years of age and the annual survival rate of female breeders is 49%: Bouwhuis et al. 2012).
Our second model followed an "extended" animal model structure, in which the identity of the male was also included as a random effect linked to the pedigree, such that both female and male additive genetic effects (and all six possible genetic covariances) were estimated for each trait, as were both male and female permanent environmental effects: for which G social can be considered as a block matrix composed of the male and female additive genetic matrices, i.e., and B (and its transposition B T ) is composed of the intersexual genetic covariances, i.e., The total additive genetic variance for each trait, V A (total ) , was estimated as where V ♀A and V ♂A are the female and male additive genetic variances, respectively, and COV ♀A,♂A is their covariance (Bouwman et al. 2010;Bijma 2014). This assumes that paired individuals are unrelated (Szulkin et al. 2007). Note that the intersexual additive genetic covariance does not represent within-pair covariance of the breeding value of a female with that of her mate but rather is the covariance between the heritable effects on the phenotypic trait of male and female relatives (weighted by their relatedness by descent: Wolak and Reid 2016).
In systems where limited dispersal creates spatial aggregation of relatives, environmental heterogeneity can potentially bias additive genetic (co)variance estimates by contributing to phenotypic resemblance of relatives, as demonstrated in this population with respect to laydate heritability estimates based on mother-daughter regression (van der Jeugd and McCleery 2002). In great tits, males settle closer to their nest-of-origin (Dhondt 1979;Greenwood et al. 1979) so bias arising through environmental covariance among relatives may be particularly important to consider when quantifying heritable variation mediated via males, and may not be adequately described by a coarse-grained approach based on a limited partitioning of the study site. To better assess the impact of spatial heterogeneity within the study area, we constructed a third bivariate animal model that modeled spatial autocorrelation. The spatial position of all nestboxes within Wytham Woods was recorded previously using a differential GPS system (Wilkin et al. 2007a), and using these coordinates we allocated nestboxes to rows and columns of width 50 m (aligned with the Ordnance Survey National Grid), yielding 71 column identities and 60 row identities (Fig. 2). To model two-dimensional spatial autocorrelation, we applied a first-order separable autoregression function that modeled distinct row and column autocorrelations, as well as fitting a variance term (Gilmour et al. 1997) that is the direct product of a variance matrix and a correlation matrix, i.e., SA = V col col ⊗ row , where V col is the column variance matrix, and col and row are the correlation matrices for the column and row models, respectively (Butler et al. 2017). This approach deals well with the high spatial heterogeneity that may be present within a forest system (as compared to other settings, such as field trials: Dutkowski et al. 2002). Following Stopher et al. (2012), we did not include spatial . Univariate models showed that corresponding spatial correlations of the two traits were of the same sign (Tables  S3 and S4), such that we felt using averaged spatial autocorrelations in our bivariate model was a reasonable solution to the difficulty of specifying autocorrelations for each trait (Gilmour 2010). Convergence issues forced us to drop the plot effect from this model, such that this model decomposes the phenotypic matrix, P, as follows: and represents variation attributable to spatial autocorrelation. All analyses were conducted using ASReml-R version 4.1.0.106 (VSN International, Hemel Hempstead) within the R framework (version 3.6.1) (R Core Team 2019) and were based on restricted maximum likelihood estimation. While clutch size is discontinuous, we do not consider it to be determined by a, for example, Poisson process, and it is thus modeled here assuming a Gaussian process (de Villemereuil 2018), as is laydate. Total narrow-sense heritability was estimated as the ratio of the total additive genetic variance relative to the within-year phenotypic variance. Fixed effects were tested by performing approximate F tests using conditional Wald statistics. For random effects, we specified a correlation structure using the corgh() function and tested correlation and variance statistics using log-likelihood ratio tests, in which the full model is compared to one in which the focal effect(s) is excluded, with the test statistic defined as twice the difference in log-likelihood between the two models and assumed to follow a chi-squared distribution. Degrees of freedom were calculated following Self and Liang (1987). Models specified using the us() function were initially used and give quantitatively similar results but restricted models in which year variances were fixed to zero (for testing) did not converge. The code for all analyses is archived in the Dryad Digital Repository, along with the primary data (Evans et al. 2020).

Results
Our dataset is based on observations from 8189 nests (with a mean of 158 per year [range: 52-306]) and represents 5638 individual females with a mean of 1.45 (range: 1-8) observations per female and 5649 individual males, likewise with a mean of 1.45 (range: 1-6) observations per male. Our "classic" bivariate animal model of laydate and clutch size, which considered only female contributions to phenotypic (co)variation, returned heritability estimates (±SE) of 19 ± 3% and 24 ± 3% for laydate and clutch size, respectively, with an inter-trait genetic correlation of -0.22 ± 0.09 (Table 1). Permanent environmental effects also contribute to individual variation in laydate (12 ± 3%) and particularly clutch size (22 ± 3%), such that individual-level heterogeneity among females is estimated to account for 31 ± 2% and 46 ± 2% of the (conditional, and within-year) phenotypic variance in laydate and clutch size, respectively. Annual variation was considerable for both traits but for laydate it was greater in magnitude than all other variance components combined (Table 1). Plot-level spatial variation received strong statistical support for both traits, indicating that this segregation of the study site captures spatial heterogeneity in reproductive behaviors, although the proportion of phenotypic variance accounted for was modest.
The extended animal model suggested that individual heterogeneity among males contributes to phenotypic variance in laydate, with a male heritability estimate of 4.6 ± 2.0% and a (statistically non-significant) male permanent environmental effect accounting for 3.2 ± 2.5% of within-year phenotypic variance. The additive genetic correlation between the female and male effects was moderate and positive, but with large standard error (r A = 0.33 ± 0.22). If we nonetheless use this best estimate of the additive genetic covariance to calculate the total additive genetic variance in laydate, we have an estimated total heritability of 28 ± 4%, which is 1.5 times greater than the estimate from the unilateral model (Fig. 3).
With respect to clutch size, the "classic" animal model considering only female contributions to phenotypic variance returned a heritability estimate of 24 ± 3%, with temporal and permanent environmental sources of variance also of considerable  Inter-trait covariances were fitted for all environmental random effects with the exception of the spatial autocorrelation effect (third model), for which terms (column and row correlations, and variance; marked with an asterisk) are averaged across the two traits and tested as a collective unit (i.e., comparing the full model to a model in which the variance and the two correlation statistics were not estimated). Variance-standardized effect sizes are relative to the within-year phenotypic variance ( V P (w−y) ), defined as the sum of all variance components except annual variance (see Methods). Trait-specific female age (first-year or older) was fitted as a fixed effect in all models; the second and third models also included male age as trait-specific fixed effects.    magnitude (Table 1). Extending the model to explicitly consider sources of variation mediated by males found little evidence for either permanent environmental (1.9 ± 2.2%) or additive genetic effects (1.3 ± 1.7%) of males on clutch size. The intersexual genetic correlation for clutch size was also moderate and positive, but estimated poorly ( r A = 0.30 ± 0.40). In contrast to laydate, there is little support for male effects-be they of heritable or environmental origin-on clutch size in this population. Our third bivariate model of laydate and clutch size replaced the plot effect with explicit quantification of the spatial autocorrelation in reproductive behavior across the study site. While the presence of spatial autocorrelation was strongly supported, this model indicated that spatial autocorrelation does not account for the heritable male effects on reproductive timing we report, since the additive genetic (co)variance estimates for laydate from this third model closely matched those of the preceding model (Table 1).

Discussion
Previous studies concluded that the contribution of individual heterogeneity among males to reproductive behavioural variation in great tits is sufficiently small that it can be ignored (van Noordwijk et al. 1980(van Noordwijk et al. , 1981aBrowne et al. 2007). However, the lack of repeat-sampling in populations of shortlived animals renders estimation of individual-level variation challenging (Browne et al. 2007), and estimation of additive genetic (co)variances can benefit from the additional replication provided by sampling of known relatives. We used a 52-year study of great tit breeding behavior to assess whether males contribute to heritable variation in laydate and clutch size, finding support for heritable male effects on laydate. That the additive genetic (co)variance estimates were robust to the inclusion of the spatial autocorrelation term suggests that they are not strongly biased by fine-scale environmental effects shared by relatives. Our analyses thus suggest that seasonal reproductive timing in great tits can be considered a joint trait (Queller 2014), subject to the influence of both members of a breeding pair. A moderate, positive correlation between the male and female genetic effects, although non-significant, suggests that heritable genetic variation in laydate may be considerably higher than indicated by 'classic' animal models that overlook the contributions of males, such that the response to selection on this trait will be more rapid than otherwise predicted, since influential female and male behaviors will co-evolve in mutually reinforcing directions. Although beyond the scope of this study, these findings indicate that future quantitative genetic studies of the contemporary evolution of reproductive timing should consider the contribution of male genetic effects if they are to fully quantify the phenotypic change attributable to heritable factors (Teplitsky et al. 2010); correspondingly, estimates of selection should incorporate any social aspects. However, an empirical comparison of the observed rate of evolution to rates predicted using female-only versus dyadic inheritance models is likely to be insightful only when estimates of both observed and predicted rates are sufficiently precise, and when the intersexual genetic correlation is so strong that predictions diverge enough to be statistically discriminable.
Clutch size represents the second major axis of variation in breeding behavior. We find that males make minor contributions to phenotypic variation in clutch size, as reported previously for this population (Browne et al. 2007) and other bird species (Sheldon et al. 2003;Charmantier et al. 2006;Auld et al. 2013;Brommer et al. 2015). Thus, the two principal axes of reproductive behavioral variation differ with respect to the extent of social determination of phenotypic expression. Nonetheless, the moderate intrasexual genetic correlations between laydate and clutch size (non-significant for males) demonstrate the integrated nature of the genetic architecture of reproductive traits and suggest that bivariate, two-sex models are necessary for accurately predicting their evolutionary dynamics (Jensen et al. 2008): any evolutionary change in one trait is expected to be accompanied by a change in the other. We note that our "classic" animal model returned a similar estimate for the inter-trait correlation in female genetic effects, such that the female additive genetic (co)variance statistics appear to have been little biased by the omission of male effects. The inter-trait genetic correlation for male effects was estimated poorly but a negative relationship, if real, would be interesting, given that any male influence is free of the physiological constraints we expect for females. However, a very large sampling effort would be required to permit definitive conclusions about the genetic correlation of male genetic effects on laydate and clutch size in our population.
The animal model readily facilitates incorporation of terms describing spatial heterogeneity. There are a number of ways this can be approached, each of which will be best suited to particular study systems, data structures, and research questions. In our case, we progressed from a typical approach to quantifying spatial heterogeneity (partitioning the study site into discrete areas that demarcate relatively homogeneous areas) to one that explicitly estimates variance attributable to spatial autocorrelation and has recently been applied in quantitative genetic studies of freeliving animal populations (Stopher et al. 2012;Regan et al. 2017). This latter approach relies on aggregating spatial locations of nestboxes within a grid, mimicking the 'block' design typically applied to agricultural and forest genetics trials (Dutkowski et al. 2002). A previous quantitative genetics study of this population, relying on mother-daughter regression, reported a marked bias in the heritability of laydate due to spatial autocorrelation (van der Jeugd and McCleery 2002). While we did not follow the same methodology in exploring this process, our results are consistent with the suggestion that the animal model is less sensitive to environmental bias than parent-offspring regression (Postma and Charmantier 2007). Nonetheless, the introduction of a term explicitly quantifying spatial heterogeneity does not guarantee that quantitative genetic parameters are estimated without bias: habitat choice may be heritable, such that part of the heritable variation in reproductive behaviors might be erroneously attributed to shared environmental effects (Postma 2006;Stopher et al. 2012).
We benefit from the accumulation of more recent data compared to earlier studies that quantified individual-level male effects on laydate and clutch size (Browne et al. 2007) and heritable male effects on laydate (van der Jeugd and McCleery 2002;Liedvogel et al. 2012) in this population. Indeed, we relied on data collected over more than half a century. Despite the scale of our sampling, intersexual genetic correlations were estimated with only moderate precision, limiting our ability to reach firm conclusions regarding the total heritable genetic variance of each trait. This is a recognized challenge for quantitative genetic studies of wild populations ) but the fact that it applies to a study relying on such an extensive dataset serves as a cautionary note for assessment of the impact of IGEs on the evolutionary dynamics of wild populations. Study populations with far more limited immigration rates (e.g., the Mandarte Island song sparrow system: Smith et al. 2006) may be preferable in this respect: if we adjust our definition of phenotypic variance to include annual variance then the standard errors of our female and male heritability estimates for laydate are very close to those of Germain et al. 2016; Table 2), despite our much larger sample size (8189 laydate observations in our study versus 1040 laydate observations for Germain et al. 2016) study). This may reflect the lower immigration rate of the insular song sparrow population, which means a greater proportion of the elements of the relatedness matrix will be non-zero. Nonetheless, the two studies report concordant results: our estimates of male and female heritabilities in laydate correspond closely with those reported by Germain et al. 2016) (allowing for their inclusion of annual variance in the phenotypic variance estimate). As also reported by Germain et al. 2016), we found in preliminary models that the natal environmental effects via females and males were marginal despite our broadly inclusive definition (relying on maternal identities to define factor levels; Tables S1 and S2).
Our analyses are consistent with a heritable influence of males on the seasonal timing of egg laying and suggest that this effect may covary positively with the female genetic effect, such that heritability estimates that do not consider male sources of heritable variation risk substantially underestimating the level of standing genetic variation for a trait that is subject to strong selection (McCleery et al. 2004). While studies published to date suggest that male heritabilities of reproductive phenology are generally small in magnitude (Table 2), the impact of male genetic variance on the total amount of heritable variation available to selection is dependent on the intersexual genetic correlation: if strongly negative then the total genetic variance will be lower than the female genetic variance (Bijma 2014). However, we found only five attempts to estimate the intersexual genetic correlation for reproductive phenology in wild populations (Table 2), and it was inestimable in one of these owing to negligible male genetic variance (Caro et al. 2009). Of the other four, Brommer and Rattiste's (2008) study of laying date in common gulls (Larus canus) is exceptional in reporting a negative intersexual genetic correlation. In a congener species, the intersexual genetic correlation was similar in magnitude but positive (Teplitsky et al. 2010), highlighting the difficulty of reaching general conclusions at present. This uncertainty occurs despite these studies focussing on a form of IGE that manifests within breeding pairs, which maximizes the number of "groups" available for analysis and is thus optimal for estimation of social genetic statistics (Bijma 2010). However, phenotypic variance is inconsistently quantified  (Brommer and Rattiste 2008;Teplitsky et al. 2010) but this is not considered in any of the studies on passerine species. The impact of model structure on additive genetic (co)variance estimation is illustrated by Liedvogel et al.'s (2012) study of laydate among great tits in Wytham Woods (including data from early years excluded in our sample), which reported female and male additive genetic variances of 2.62 ± 0.67 and 0.75 ± 0.41, respectively, based on an animal model that included breeding altitude and local breeding density as fixed effects, but of 16.41 ± 1.06 and 15.15 ± 1.16 (V ♀A and V ♂A , respectively) when fixed effects and other random effects were removed. Heritable variation is spatially non-random in our population (Garant et al. 2005) and, with respect to reproductive timing, heritable variation could covary with the fine-scale environmental variation described by altitude and breeding density (Wilkin et al. 2007b;Wilkin et al. 2006). Indeed, Liedvogel et al. 2012) report estimates from their main model (i.e., including fixed effects) for both female and male additive genetic variance that are considerably lower than those we report here, consistent with such an explanation for the discrepancies between analyses with broadly overlapping sample populations. Clearly, incorporating the impact of the social environment into empirical description of evolutionary dynamics under natural settings is a considerable challenge, and it is perhaps unsurprising that published estimates of direct-indirect genetic correlations for free-living populations are exceedingly few (Fisher and Pruitt 2019;Hunt et al. 2019), despite the increasing application of quantitative genetic analyses to long-term study populations more generally (Postma 2014). Nonetheless, while the extent to which indirect genetic determination applies to behavioral traits remains unclear, recent results from laboratory studies suggest it is potentially widespread (Carter et al. 2019;Dewan et al. 2019). Of course, expression of a host of behaviors is sensitive to the phenotype of unrelated social partners, as testified by studies of mate choice (Bateson 1983), resource competition (Hardy and Briffa 2013), and social structure (Krause et al. 2015). Indeed, social partners have been shown to impact on an individual's own phenotype in our study population (Firth et al. 2015b), suggesting a potential for IGEs to also operate across larger and more dynamic social networks (Firth et al. 2017), and demonstrate the potential for heritable effects of residential neighbors in a population of North American red squirrels (Tamiasciurus hudsonicus). Agonistic encounters, which can often be reduced to dyadic interactions, are another social context in which we should expect IGEs to be prominent (Santostefano et al. 2017), since contest outcomes are typically determined by the ability of protagonists relative to their opponents (Wilson et al. 2011;Sartori and Mantovani 2013). As already discussed, quantification of IGEs is data demanding (Ellen et al. 2014), such that the availability of suitable datasets may be the limiting factor in expanding our empirical knowledge, particularly as the members of large groups are unlikely to interact uniformly with all groupmates. That said, the adoption of automated tracking technologies to study social interactions within animal groups is increasingly generating data of the type that would be required for assessing social genetic parameters within a single, large group (Hamede et al. 2009;Rutz et al. 2012;Godfrey et al. 2014;König et al. 2015), which represents a promising development for efforts to improve our understanding of indirect genetic effects in wild populations. Certainly, more work is needed if we are to deliver an informed perspective on the extent to which the evolutionary dynamics of wild populations are socially regulated.