Maria Moiron
email: mariamoironc@gmail.com

Yimen G Araya-Ajoy

Céline Teplitsky

Sandra Bouwhuis

Anne Charmantier

Understanding the Social Dynamics of Breeding Phenology: Indirect Genetic Effects and Assortative Mating in a Long-Distance Migrant

Keywords: animal model, arrival date, common tern, heritability, laying date, trait-based

published or not. The documents may come

INTRODUCTION

One of the major goals in evolutionary biology is to predict micro-evolutionary responses [START_REF] Falconer | Introduction to Quantitative Genetics[END_REF]. Such responses result from the interplay between natural selection and heritable genetic variation underlying phenotypic trait variation. Variation in the phenotypic expression of a trait is, however, composed not only of genetic, but also of environmental effects.

Since most organisms interact with conspecifics, these environmental effects may often include social effects, e.g. the phenotype of a partner has the potential to influence the phenotype of a focal individual, creating an indirect effect [START_REF] Griffing | Selection in Reference to Biological Groups I. Individual and Group Selection Applied to Populations of Unordered Groups[END_REF][START_REF] Moore | Interacting Phenotypes and the Evolutionary Process: I. Direct and Indirect Genetic Effects of Social Interactions[END_REF]). Interestingly, this indirect effect may itself be partially determined by a genetic as well as an environmental component [START_REF] Griffing | Selection in Reference to Biological Groups I. Individual and Group Selection Applied to Populations of Unordered Groups[END_REF][START_REF] Moore | Interacting Phenotypes and the Evolutionary Process: I. Direct and Indirect Genetic Effects of Social Interactions[END_REF]. Evolutionary theory postulates that the estimation of the rate and direction of any micro-evolutionary change for a given trait will therefore depend on both direct and indirect genetic variation in this trait and on the genetic covariance between the two, as well as on the magnitude and direction of selection.

The timing of reproduction is a trait of special interest in the face of global and climate change, and the potential of reproductive phenology to adapt to rapidly changing environmental conditions has received substantial attention [START_REF] Parmesan | Influences of species, latitudes and methodologies on estimates of phenological response to global warming[END_REF]. In birds, traits such as laying date have been shown to be under strong directional selection (i.e., for earlier laying, Charmantier and Gienapp 2014), and, in some cases, it has been demonstrated that recent climate change is a driver of selection on these traits [START_REF] Marrot | Current spring warming as a driver of selection on reproductive timing in a wild passerine[END_REF]. To predict evolutionary responses to selection on laying date, it may be necessary to quantify both female and male genetic effects, because partner effects may occur when pair members affect each other's pre-breeding environment. For instance, laying date is often considered a female-specific trait (but see [START_REF] Teplitsky | Indirect genetic effects in a sexlimited trait: The case of breeding time in red-billed gulls[END_REF], however males can affect the female pre-laying environment through the timing of mating or through resource provisioning (e.g., via courtship feeding, [START_REF] Nisbet | Courtship-feeding, Egg-size and Breeding Success in Common Terns[END_REF][START_REF] Tasker | A Functional Analysis of Courtship Feeding in the Red-Billed Gull, Larus Novaehollandiae Scopulinus[END_REF]. As such, laying date is likely to be influenced by genetic and environmental effects associated with both the female and the male [START_REF] Brommer | Hidden" reproductive conflict between mates in a wild bird population[END_REF][START_REF] Germain | Direct and indirect genetic and fine-scale location effects on breeding date in song sparrows[END_REF].

Several studies of male effects on female breeding phenology have simultaneously estimated female and male effects on laying date at the genetic level. Non-zero male additive genetic variances were found in two seabirds with courtship feeding behavior (h 2 =0.048 in common gulls, Larus canus, Brommer and Rattiste 2008; h 2 =0.134 in red-billed gulls, Larus novaehollandiae scopulinus, [START_REF] Teplitsky | Indirect genetic effects in a sexlimited trait: The case of breeding time in red-billed gulls[END_REF]) and in a songbird (h 2 =0.02, song sparrow, Melospiza melodia, [START_REF] Germain | Direct and indirect genetic and fine-scale location effects on breeding date in song sparrows[END_REF]. Two other studies on songbirds, however, estimated male additive genetic variances that were not different from zero (blue tits, Cyanistes caeruleus, [START_REF] Caro | Local adaptation of timing of reproduction: Females are in the driver's seat[END_REF]; great tits, Parus major, [START_REF] Liedvogel | Integrating candidate gene and quantitative genetic approaches to understand variation in timing of breeding in wild tit populations[END_REF]. As such, results to date are mixed and quantifying male additive genetic variance in more systems is needed to evaluate the conditions under which such effects are important. Additionally, the mechanisms and causal processes underlying the observed male effects have remained unexplored. No study has investigated which male trait might underpin the reported male genetic effects and we advocate moving from purely descriptive work to a more integrative framework to not only describe quantitative genetic patterns, but also unveil the potential pathways driving them in nature.

Migratory species, which are specifically vulnerable to the global and climate changes mentioned above, face particular challenges since breeding can only start once the two partners have arrived (and recovered) from their migration. Hence, arrival date strongly determines when pairs can form and initiate reproduction. Additionally, arrival date of the male can also influence laying date in the female because those individuals that arrive early from the wintering grounds are birds of high quality (i.e. that live longer, Zhang et al. 2015a; and provision at a higher rate, [START_REF] Limmer | The relative role of age and experience in determining variation in body mass during the early breeding career of the Common Tern (Sterna hirundo)[END_REF], that may be better able to secure good foraging territories. The foraging ability of the male could then affect when its partner attains the necessary breeding condition (via courtship feeding behavior), and, ultimately, affect the onset of reproduction for the pair. This tight link between timing of migration and timing of breeding creates special circumstances for indirect genetic effects (IGEs) in migratory species so that laying date might not only be determined by the female and male phenotype and genotype with respect to reproductive phenology (direct and indirect effects), but also by their migratory phenotypes and genotypes, as well as the covariance between them. Importantly, the covariance between the phenotype of interacting individuals can arise because of different mechanisms, complicating the estimation of IGEs. This is because such a genetic covariance could arise due a plastic response of the trait of the female to the genotype of the male (IGEs), but also due to a non-random assortment of genotypes (assortative mating).

Assortative mating for arrival date indicates that females that arrive early from the wintering grounds breed with males that also arrive early [START_REF] Ludwig | Supply and demand: causes and consequences of assortative mating in common terns Sterna hirundo[END_REF]. Assortative mating might confound the estimation of male effects on female breeding timing because it generates a covariance between the breeding values of the pair members. This is problematic because a variance-partitioning model fitted to test for IGEs is only explicitly modelling the covariance between direct and indirect effects; however, in the presence of assortative mating, some of the covariance caused by assortative mating will spill over to the component estimated as an indirect effect, ultimately inflating the estimate of IGEs. We made use of data simulations to investigate this added layer of complexity with the aim of exploring the magnitude of the potential bias that assortative mating might introduce to estimates of IGEs.

In this study, we used 26 years of data from a pedigreed population of a long distance migratory seabird, the common tern (Sterna hirundo) to investigate the extent to which male environmental and genetic effects influenced female laying date. We expected such partner effects, because (i) common terns are income breeders, such that the energy females require to produce eggs must be acquired from current resources, (ii) males provide their partners with courtship food during the pre-laying stage [START_REF] Nisbet | Courtship-feeding, Egg-size and Breeding Success in Common Terns[END_REF]) and (iii) females in good condition are known to breed earlier and more successfully [START_REF] Wendeln | Body mass change in breeding Common Terns Sterna hirundo[END_REF]. As such, the first objective of the study was to explore the evolutionary genetics of female laying date and the associated partner effects.

To do so, we quantified the existence of male (environmental and genetic) effects on female laying date using a variance-partitioning approach. These analyses were then complemented with a phenotypic trait-based approach investigating whether male arrival date could (partially) explain the observed partner effects on female laying date. Applying this phenotypic trait-based approach also enabled us to quantify sources of (co)variation in the arrival dates of the breeding pair as well as the direct and indirect genetic effects on female laying date that are independent of arrival date.

As a second objective, we investigated whether assortative mating could influence the observed partner effects. We first quantified patterns of assortative mating in our population of common terns, and then investigated its influence on male variance components using data simulations.

Using these complementary IGE approaches plus data simulations, we aim to effectively dissect the evolutionary causation of a complex social trait in a long distance migratory species.

METHODS

Study System

The data came from a long-term study population of common terns. The colony site is located at the Banter See on the German North Sea coast (53°36´N, 08°06´E) and consists of six concrete islands, each of which is surrounded by a 60 cm wall. The individual-based study was initiated in 1992, when 101 adult birds were caught and marked with individually numbered subcutaneously injected transponders. Since 1992, all locally hatched birds have similarly been marked with a transponder shortly prior to fledging and the presence and reproductive performance of marked individuals has been monitored following a standard protocol [START_REF] Becker | A New Application for in Population Ecology of the Common Tern[END_REF]. Since 1994, the arrival of marked terns has been monitored with antennae on resting platforms affixed to the walls of the islands.

Common terns arrive from their wintering grounds (in subequatorial Africa; [START_REF] Becker | Common Terns on the East Atlantic Flyway: temporal-spatial distribution during the non-breeding period[END_REF] from April, lay their first eggs from early May, and return to their wintering grounds from the end of August [START_REF] Becker | Common Terns on the East Atlantic Flyway: temporal-spatial distribution during the non-breeding period[END_REF]. Thanks to the antenna system, it is possible to record the exact arrival date of every marked bird in the population. Arrival date is thus defined as the day of first return to the breeding grounds (January 1 st =1) (Zhang et al. 2015b). As part of the standard protocol, the colony is checked for new clutches every 2-3 days throughout the breeding season.

These checks allow identification of clutch initiation, and laying date is defined as the date of appearance of the first egg in each female's first clutch of the season (January 1 st =1) (Zhang et al. 2015b). Parents are identified using portable antennae placed around each nest for 1-2 days during incubation, which is shared by both partners. In common terns, second clutches are extremely rare [START_REF] Moore | The Production of Second Clutches in the Common Tern: Proximate Effects of Timing and Food Supply[END_REF], even though replacement clutches are regularly produced if the first clutch or brood is lost [START_REF] Wendeln | Parental care of replacement clutches in common terns (Sterna hirundo)[END_REF][START_REF] Becker | Renesting of Common Terns Sterna hirundo in the life history perspective[END_REF]. Data are deposited in the Dryad Digital Repository: https://doi.org/10.5061/dryad.5mkkwh73c [START_REF] Moiron | Data from: Understanding the social dynamics of breeding phenology: indirect genetic effects and assortative mating in a long distance migrant[END_REF]).

The social pedigree was constructed from the observations of parents and their fledged offspring. Due to the low levels of extra-pair paternity in our population, (97.1% of true siblings among 22 broods, [START_REF] González-Solís | Courtship feedings, copulations and paternity in common terns, Sterna hirundo[END_REF], the social pedigree is a good approximation of the genetic pedigree. For the period 1992-2017, the pruned pedigree comprised 1553 individuals. The maximum depth was 5 generations, the number of paternities and maternities 539 and 566, respectively.

Quantitative Genetics

We used a variance-partitioning approach to estimate female and male environmental and genetic effects on laying date. To do so, we fitted two linear mixed-effects models to our data (Table 1).

The first model estimated female and male effects at the individual level ("phenotypic model", model 1), while the second model decomposed these effects into genetic and environmental components ("genetic model", model 3).

For the "phenotypic model", we built a univariate model that simultaneously estimated individual variance components attributed to female and male individual effects on female laying date (VI-female and VI-male, respectively), reflecting permanent environmental and/or genetic effects.

To build the phenotypic model, we fitted female laying date (mean centered and variance standardized) as a response variable, and random intercepts for year (to control for temporal variation), female identity (to control for repeated measures of females across years) and male (social partner) identity. We did not include pair identity as a random effect because of the overlapping variation explained by female and male identities with pair identity (common terns have, on average, less than two partners in their lifetime, see below). As fixed effects, we fitted linear and quadratic female age effects (both mean centered and variance standardized; [START_REF] Ezard | Correlations between age, phenotype, and individual contribution to population growth in common terns[END_REF]Zhang et al. 2015b). We removed those individuals from which the age was not known (n = 190 observations on 27 females first caught as adults). Results from fixed effects estimates are reported in the Supplementary Material (Table 2 andTables S1-S2). Thus, for the "phenotypic model" we assumed that the laying date y of female i paired with male j in year k was described by the model

𝑦 𝑖𝑗𝑘 = 𝛽 0 + 𝛽 1 𝑎𝑔𝑒 𝑖 + 𝛽 2 𝑎𝑔𝑒 𝑖 2 + 𝐼 𝑖 + 𝐼 𝑗 + 𝑌 𝐾 + 𝑒 𝑖𝑗𝑘 (1a)
where 𝛽0 was the overall population mean, 𝛽 1 the linear effect of the age (𝑎𝑔𝑒 𝑖 ) of female i and,

𝛽 2 the quadratic effect of female age 𝑎𝑔𝑒 𝑖 2 . We also modeled individual deviations from the population mean 𝐼 for female i and the individual deviations of the effects of her male social partner 𝐼 𝑗 . 𝑌 𝐾 reflected yearly variation in the mean phenotype of the population and 𝑒 𝑖𝑗𝑘 , the residual (unexplained) variation. These parameters were fitted as random effects and variances were estimated assuming a normal distribution with a mean of zero.

For the "genetic model", we split the individual effects of males and females into additive genetic and permanent individual effects. To do so, we built a univariate animal model that decomposed the total phenotypic variance of laying date into female and male additive genetic variances (VA-female, VA-male), female and male permanent individual variances (VPI-female and VPImale), and overall year (VYR) and residual (VR) variances. The simultaneous estimation of additive genetic variances in female and male laying date also allowed us to quantify the additive genetic covariance and correlation for the female and male effects (CovA-female,male and rA-female,male, respectively). We did not fit the covariance between male and female permanent individual effects or residuals because there can be no cross-sex covariance for individual identity or residual components in sex-limited traits. As before, we fitted female laying date (mean centered and variance standardized) as a response variable and random intercepts for year, female and male identity, while adding an additional female and male identity linked to the pedigree. Linking the identities of both female and male to the pedigree was possible by adding a second "animal" term relating the identity of the social mate to variation in the phenotype of the focal individual. Thus, this second "animal" term estimates how much of the variation in the focal phenotype is explained by additive genetic effects of the social partner. As fixed effects, we fitted linear and quadratic female age effects (both mean centered and variance standardized; [START_REF] Ezard | Correlations between age, phenotype, and individual contribution to population growth in common terns[END_REF]Zhang et al. 2015b). Thus, for the "genetic model" we assumed that the laying date y of female i paired with male j in year k was described by the model

𝑦 𝑖𝑗𝑘 = 𝛽 0 + 𝛽 1 𝑎𝑔𝑒 𝑖 + 𝛽 2 𝑎𝑔𝑒 𝑖 2 + 𝑎 𝑖 + 𝑝𝑒 𝑖 + 𝑎 𝑗 + 𝑝𝑒 𝑗 + 𝑌 𝐾 + 𝑒 𝑖𝑗𝑘 (1b)
where 𝛽0 was the overall population mean, 𝛽 1 𝑎𝑔𝑒 𝑖 and 𝛽 2 𝑎𝑔𝑒 𝑖 2 the linear and quadratic effects of the age of female i, respectively. We also modeled female and male additive genetic components (ai and aj, respectively), and female and partner permanent individual components (pei and pej, respectively). 𝑌 𝐾 reflected yearly variation in the mean phenotype of the population and 𝑒 𝑖𝑗𝑘 , the residual error. These parameters were fitted as random effects and variances were estimated assuming a normal distribution with a mean of zero. Unstandardized variance components for the "genetic model" are reported in Table S6.

In our population, common terns often maintain their partner across years, although mismatched arrival, divorce and mortality sometimes result in the formation of new pairs [START_REF] Rebke | Better the devil you know: Common terns stay with a previous partner although pair bond duration does not affect breeding output[END_REF]. In our data, individual common terns had slightly less than two partners (average ± SD: females = 1.87 ± 1.21 [range: 1-8], males = 1.74 ± 1.00 [range: 1-8]), while 381 (52%)

females and 435 (53%) males had only a single partner. We omitted those pairs that had only one partner from the data to more reliably separate male and female effects on laying date. Including all available data (i.e., individuals with one partner or more) in the analyses, however, yielded qualitatively similar estimates of VA-male (although estimates of VI-male and VPI-male were lower, Table S3).

Once we quantified environmental and genetic effects using a variance-partitioning approach, we proceeded to test the hypothesis that male arrival date underlies (part of) the observed partner effects on female laying date. To do so, we applied a complementary approach: the trait-based approach at the phenotypic level [START_REF] Bijma | The quantitative genetics of indirect genetic effects: a selective review of modelling issues[END_REF]. This approach is used to investigate which traits in a social partner are causing indirect effects on the phenotype of the focal individual, but can only be applied when the individual-level characteristics hypothesized to cause indirect effects are measured. We fitted male arrival date as a standardized linear fixed effect to the two previous univariate (phenotypic and genetic) models (models 2 and 4) and assessed whether this reduced the male individual, environmental or genetic effect, as expected if male arrival date is indeed an important driver of variation in female laying date. Because of the strong phenotypic correlation between male and female arrival dates (see below) and the knock-on effects of female arrival date on female laying date, we also fitted female arrival date as a standardized linear fixed effect.

Additionally, we modelled female and male age (linear and quadratic effects, mean centered and variance standardized) as fixed effects. This way, we could effectively test for any independent effect of male arrival date on female laying date (models 2 and 4) and quantify the female additive genetic variance in laying date after correcting for female arrival date, potentially detecting other biological processes, such as genetic variation in the speed with which females can lay their first egg after arriving from the wintering grounds. In principle, it would have been desirable to model a trait-based indirect genetic effects approach (McGlothlin and Brodie 2009). However, we did not apply this approach because it is conceptually difficult to treat male laying date as a phenotype, when ultimately laying date is only expressed by the female. Additionally, fitting such a model would result in a redundant estimation of some of the indirect effect components (i.e., direct variance of male laying date would equal to indirect variance of male in female laying date). Thus, modelling a trait-based indirect genetic effects approach would be insightful in a general sense, however, in our specific case, its implementation would not have changed the interpretation of our results.

All models were fitted using a Bayesian framework implemented in the statistical environment R (v. 3.6.1, Team R Core 2019) and the R-package MCMCglmm [START_REF] Hadfield | MCMC Methods for Multi-response Generalized Linear Mixed Models : The MCMCglmm R Package[END_REF]. For all models, the minimum effective sample size per model was 2000. We used parameter-expanded priors with an inverse Gamma prior distribution (chosen to be non-informative). We modelled female laying date with a Gaussian error structure. Posterior point estimates and 95% Credible Intervals (95% CI) were estimated for (co)variances and correlations. We considered a random effect to be supported if the lower limit of the posterior distribution was not zero or leaning towards zero. For fixed effects, estimates for which the 95% CI did not include zero were considered to support the existence of an effect.

Assortative mating analysis and data simulations

To explore the degree to which assortative mating for arrival date might influence the estimation of male effects on female laying date, we first quantified VI-male in laying date with a univariate linear mixed-model considering laying date as a male trait, to then compare the estimates of variance components with the ones obtained from the model where both VI-female and VI-male were estimated simultaneously (model 1). We did not investigate the potential bias in VA-male due to assortative mating because our estimate of VA-male was close to zero (see below). We hypothesized that if the VI-male estimate from the 'separate model' was larger than that estimated in the 'simultaneous model', our estimate of VI-male could be biased due to unmodeled processes such as assortative mating [START_REF] Reid | Female and male genetic effects on offspring paternity: Additive genetic (co)variances in female extra-pair reproduction and male paternity success in song sparrows (melospiza melodia)[END_REF].

Second, we investigated whether our estimate of VI-male in laying date could still be affected by assortative mating even after controlling for both effects in a single model (model 1). To do so, we quantified patterns of assortative mating for arrival date. We built a bivariate model in which female and male arrival date (both mean centered and variance standardized) were fitted as response variables. We included random effects of year and pair identity [START_REF] Class | A statistical methodology for estimating assortative mating for phenotypic traits that are labile or measured with error[END_REF], and sex-specific linear and quadratic components of age (both mean centered and variance standardized) as fixed effects.

In addition, we simulated datasets where partner effects were always set to zero but with varying levels of assortative mating, and ran the same statistical models to quantify the existence of partner effects. Assortative mating was simulated as the (positive) correlation between individual-specific values across mated pairs and ranged from random (r = 0) to perfectly assortative (r = 1) at an interval of 0.2 (including the parameter estimate extracted from our model, Table 3). To model a scenario that mimicked our own empirical question, we used the empirical maximum annual number of breeding pairs (729), number of years (26), and overall survival probabilities (0.75) (Zhang et al. 2015a).

As a first step, we simulated individual-specific data (we implemented our data simulation at the individual-level because we found male additive genetic effects in laying date close to zero, Table 1, model 3). We started by generating the individuals forming pairs. We simulated data so that the number of males and females was constant over the years: from one year to the other, 75% of males and females survived and males and females that did not survive were replaced by the same number of new males and females. Thus, given that both members were alive, pairs re-mated with a 50% probability. For the remaining males and females, half died and half re-paired with new immigrants. In our simulations, two divorced individuals could not re-mate with each other, i.e., divorced individuals could only mate with new individuals. As a second step, we created an individual's phenotype by adding observation-specific environmental effects. These environmental effects were uncorrelated between pair members. It is important to note that individuals mated based on their individual-specific values, resulting in a covariance between individuals-specific values for the mated pair (instead of a covariance between the observed phenotypic value of the female and the individual specific value of the male). All simulated effects followed a Gaussian distribution. As a third step, we estimated partner effects by implementing the same statistical model (model 1) used in our main analysis to the simulated data, and compared

the effects obtained to investigate potential biases in our partner effect estimates. We replicated the data simulation process 100 times to obtain 100 independent datasets and, therefore, estimates of male effects per level of assortative mating. We then quantitatively assessed the bias caused by assortative mating given that partner effects were set to zero. The R-code for data simulations is available from GitHub: https://github.com/MariaMoiron/Moiron_etal_IGES_common_terns.

RESULTS

The data comprised of 1106 breeding events by 544 breeding pairs (omitting those pairs that only had one partner in their lifetime) consisting of 209 females paired with 406 males of known identity. The overall mean (±SD) Julian arrival date across the 26 years (1992-2017) included in our dataset was 120 ± 14.04 (April 30), and laying date was 142 ±14.14 (May 22).

Male effects on female laying date and underlying mechanisms

Common terns showed repeatable differences in laying date: female identity accounted for 23.9%

of the total phenotypic variance, while their partners' identity accounted for 16.9% (Table 1). The proportion of the total phenotypic variance in laying date explained by additive genetic effects was 12.3% for females and 3.1% for males (Table 1). The covariance and correlation between female and male additive genetic components were close zero, with wide 95% CI (Table 1).

The phenotypic trait-based analysis demonstrated that male arrival date strongly influenced female laying date (Tables 1): a late arrival of the male was linked to a late laying date of the female (Table 2). Furthermore, adding male arrival date as a fixed effect led to a four-fold decrease in male individual variance, i.e. from 0.12 to 0.03 (VI-male in models 1 versus 2, Table1). While the decrease was substantial, the 95% CI of both estimates slightly overlapped, calling for caution when interpreting its significance of the reduction in the variance explained. At the genetic level, adding male arrival date did not affect the male genetic component (because it was already close to zero), but it did reduce the male permanent individual component from 0.09 to 0.02 (the 95% CI of these estimates were slightly overlapping, VPI-male in models 3 versus 4, Table 1).

Furthermore, adding female arrival date as a fixed effect allowed the estimation of female genetic effects that were independent from female and male arrival date, reducing the magnitude of the variance in female laying date(VI-female in models 1 versus 2; VA-female in models 3 versus 4, Table 1).

Assortative mating for arrival date

The bivariate analysis of arrival date of breeding pairs showed a positive covariance between female and male arrival dates (Table 3). Scaling the covariance by VI-female and VI-male, we found the among-pair correlation to be r = 0.77 (95% CI = 0.71-0.81), demonstrating the existence of strong assortative mating. We also observed a strong, positive correlation at the within-pair level (within-pairs r = 0.28, 95% CI = 0.25 -0.33). This within-pair (residual) correlation is not indicative of assortative mating and likely arises due to shared environmental effects (e.g., when the male and female react similarly to varying annual conditions, such as the temperature of the sea surface in the wintering grounds).

The model that simultaneously estimated VI-female and VI-male in laying date returned a more than two-fold lower estimate of male effects than the model that modelled them separately (VI-male estimated in a separate model = 0.13, Table S4; VI-male estimated from the full dataset model = 0.06, Table S3). This suggests that our estimates of male effects could be biased by assortative mating. However, simultaneously estimating male and female effects does not ensure that we fully control for the effect of assortative mating. Indeed, when we simulated data in which male effects on laying date were set to zero, but assortative mating occurred at the empirically observed level, we found a male effect that was within the range of magnitude as the one we estimated from the empirical data (empirical VI-male = 0.12, data simulation VI-male = 0.07) (Table 1, Fig. 1).

DISCUSSION

We demonstrated that female common terns show repeatable differences in their laying date, and that these differences are underpinned by both permanent individual and genetic components (Table 1). Our heritability estimate of female laying date (h 2 = 0.12 ± 0.08) was lower than that previously reported for the same common tern population based on a different subset of data (0.27 ± 0.09; [START_REF] Dobson | Plasticity results in delayed breeding in a long-distant migrant seabird[END_REF]), but in line with that reported for other bird species (e.g., h 2 = 0.19 in collared flycatchers, [START_REF] Sheldon | Natural selection and inheritance of breeding time and clutch size in the collared flycatcher[END_REF]0.16 in great tits, McCleery et al. 2004;0.15 in common gulls, Brommer and Rattiste 2008).

We expected males to affect the female pre-laying environment, mainly due to (repeatable and/or heritable) variation in arrival date. In line with this prediction, we showed that male effects explained variation in female laying date, and that these effects were underpinned by a permanent individual, and, to a lower extent, genetic, component. As we argue below, caution is, however, required when interpreting these findings, since we cannot exclude the possibility that the observed effects are mostly caused, or at least partially inflated, by assortative mating. Our estimate of the male genetic effect on female laying date (h 2 =0.03 ± 0.04) is somewhat lower than that of two other studies in seabirds (h 2 =0.05 in common gulls, Brommer and Rattiste 2008; h 2 =0.13 in redbilled gulls, [START_REF] Teplitsky | Indirect genetic effects in a sexlimited trait: The case of breeding time in red-billed gulls[END_REF]), but in line with a study in passerines (h 2 =0.02 in song sparrows, [START_REF] Germain | Direct and indirect genetic and fine-scale location effects on breeding date in song sparrows[END_REF], overall supporting the notion that traits are often heritable only in the sex in which they are expressed (e.g. [START_REF] Kruuk | Antler size in red deer: heritability and selection but no evolution[END_REF][START_REF] Sheldon | Natural selection and inheritance of breeding time and clutch size in the collared flycatcher[END_REF][START_REF] Charmantier | Evolutionary Response to Selection on Clutch Size in a Long-Term Study of the Mute Swan[END_REF].

Given that laying date in common terns is mostly a sex-limited trait and the male-female genetic covariance was close to zero, we can conclude that its evolutionary trajectory should almost exclusively be determined by the genetics of the female and selection acting on females. Hence, while accounting for partner effects may be important to quantify the evolutionary potential in some species, in our case accounting for partner effects should not affect the estimates of the evolutionary potential of laying date.

The importance of indirect genetic effects (IGE) has long been acknowledged in breeding programs. However, empirical IGE studies of diverse phenotypes in wild animal systems have only emerged in recent years. Among the studied traits, the timing of reproduction has attracted large attention, mostly in birds and mammals. Overall, these studies of IGEs in breeding timing have demonstrated that there is varying potential for male effects (including a case where laying date was not heritable in females but significantly so in males, [START_REF] Teplitsky | Indirect genetic effects in a sexlimited trait: The case of breeding time in red-billed gulls[END_REF]). However, most previous studies focused on species that are resident, or at least, not long distance migrants.

Migratory species pose a special challenge, as timing of breeding is tightly linked to timing of migration. Hence, understanding the quantitative genetics of breeding timing in migratory species requires quantifying the female and male phenotype and genotype for laying date, their migratory phenotypes and genotypes, as well as the covariance between them. Although it may be challenging, adding this layer of complexity will help to further our understanding of the whole complexity of partner effects, particularly in the case of migratory species. More generally, the mechanisms and causal processes underlying the observed male effects have remained unexplored.

As such, the time is ripe to start investigating the complex causal pathways underlying the existence of IGEs. This task is not trivial, since it requires not only collecting phenotypic and genetic data on the focal individual, but also on their partners.

In this study, we considered male arrival date to be a good candidate to underlie partner effects on female laying date, because birds that arrive earlier from the wintering grounds are known to be birds of high 'quality'. They, for example, live longer (Zhang et al. 2015a) and provision their offspring at a higher rate [START_REF] Limmer | The relative role of age and experience in determining variation in body mass during the early breeding career of the Common Tern (Sterna hirundo)[END_REF], which suggests they may be better able to locate and explore good foraging areas, which might also be reflected in better courtship feeding that gets their partners into breeding condition more efficiently. As expected, we found that female common terns differed in their laying date in response to the phenotypic value of their males' arrival date (Table 2): an earlier arrival of the males was linked to an earlier laying date of their females. This effect of male arrival date was corroborated by the decrease in variance explained by male identity when arrival date was fitted as a covariate (even though 95% CI of the effect sizes are slightly overlapping, Table 1). As we did not find that laying date was substantially affected by the genetic component of the partner's identity, it indicates that female common terns are mainly influenced by environmental components of male arrival date. Indeed, the permanent individual component of male arrival date explained 24.3% of the total variance in the trait, while the additive genetic component explained only 4.3% (Table S5).

In our phenotypic trait-based analysis, we also modelled the effect of female arrival date on female laying date. By doing so, we simultaneously confirmed the expected knock-on effect (an earlier arrival of the female was linked to an earlier laying date, Table 2) and estimated the female genetic effect on laying date that is independent from female and male arrival date. After correcting for differences in arrival date, females still differed in when they start breeding. This means that other mechanisms, independent of migratory phenology, are underlying variation in laying date among females. The potential behavioral or physiological mechanisms underpinning these differences remain unknown.

We then asked whether unmodeled environmental effects could underlie or inflate estimates of male effects. In case of assortative mating, for example, we presumed estimates of partner effects on female traits to be biased because of the variance-partitioning approach being unable to disentangle the true partner component and the female component. Therefore, we investigated patterns of assortative mating in our population, with the aim of quantifying the degree to which non-random mating may influence the observed partner effects. We confirmed that individuals were strongly assortatively mated with respect to their arrival date (Table 3; also see [START_REF] Ludwig | Supply and demand: causes and consequences of assortative mating in common terns Sterna hirundo[END_REF]. Additionally, we observed that the model that simultaneously estimated female and male individual effects on laying date returned a lower estimate of male effects than did separate models (considering laying date as a male trait, Table S3 v. S4, [START_REF] Bijma | Multilevel selection 2: Estimating the genetic parameters determining inheritance and response to selection[END_REF][START_REF] Bouwman | Maternal and social genetic effects on average daily gain of piglets from birth until weaning[END_REF][START_REF] Reid | Female and male genetic effects on offspring paternity: Additive genetic (co)variances in female extra-pair reproduction and male paternity success in song sparrows (melospiza melodia)[END_REF][START_REF] Wolak | Is Pairing with a Relative Heritable? Estimating Female and Male Genetic Contributions to the Degree of Biparental Inbreeding in Song Sparrows ( Melospiza melodia )[END_REF]. This suggests that unmodeled environmental effects (possibly assortative mating) were indeed influencing our estimate of partner effects, rendering it necessary to simultaneously model female and male effects. While the results of our model comparison can be taken as an indicator of assortative mating inflating our estimates of partner effects, the evidence is not conclusive, because the existence of other unmodeled environmental factors cannot be excluded. We therefore compared our empirical results with those obtained from analyses of simulated data and showed that, given our data structure and the observed among-pair correlation for arrival date, assortative mating could largely explain our estimate of male individual variance on female laying date (Fig. 1).

Our simulation approach aimed to illustrate that quantifying male effects on female phenology becomes complicated when assortative mating occurs. The current lack of empirical work on the link between IGEs and assortative mating, however, hampers the possibility to draw general conclusions on how widespread this bias might be. It also makes it difficult to speculate how our findings might retrospectively change conclusions from previous studies on IGEs (i.e., whether other studies might also suffer from an overestimation of IGEs due to the existence of assortative mating or other shared environmental effects in the population). Importantly, assortative mating might not be the only mechanism that could inflate the estimation of IGEs. Any common environmental effects that is not explicitly accounted for in the statistical model might be captured by another term in the model, as, for instance, the IGE component (e.g., [START_REF] Germain | Direct and indirect genetic and fine-scale location effects on breeding date in song sparrows[END_REF]. Such scenarios will need to be explored if we are to understand the relevance of indirect (genetic) effects in nature.

Altogether, our study provides timely evidence for female and male genetic and environmental variation in an important phenological trait, laying date, in a long distance migrant.

Importantly, we also investigate the mechanism underlying these observed partner effects, identifying male arrival date as an important trait, while exploring the evolutionary consequences of a covariance between the migratory phenotype/genotype of pair members. Data simulations showed that assortative mating can indeed generate partner effects of similar magnitude as our empirical estimate of male individual variance in female lay date. As such, our findings highlight that, although it may be challenging, to further our understanding of the nature of partner effects, there is a need to evaluate the potential bias that can be caused by assortative mating, or any other common environmental effect that increases the phenotypic resemblance among pair members.

FIGURE

Figure 1. Male individual effects on female laying date as a function of the strength of assortative mating. Models were fitted to simulated datasets that parameterized our data structure and where VI-male was set to zero.

ACKNOWLEDGEMENTS:

We are indebted to Peter H. Becker for setting up the long-term common tern population study and thank the numerous field workers that have contributed to compiling the data. We thank Barbara Class, Jon Brommer, Niels Dingemanse, Jane Reid, and Jarrod Hadfield for discussions.

We also two anonymous reviewers, Tim Connallon and Daniel Bolnick for comments that greatly improved the manuscript. M.M. was funded by a Marie Curie Individual Fellowship (PLASTIC TERN, Grant Agreement Number 793550), A.C. by the European Research Council (Starting grant ERC-2013-StG-337365-SHE).

AUTHOR CONTRIBUTIONS:

The study was conceived by M.M. and Y.A-J. and developed with input from all authors. S.B. managed the data collection and database. M.M. and Y.A-J. designed the analysis protocol. M.M. performed the statistical analyses, and M.M. and Y.A-J, the data simulations. M.M. drafted the manuscript with input from all co-authors.

TABLES

Table 1. Variance component estimates (posterior point estimates with associated 95% Credible Intervals) from "phenotypic" and "genetic" models. Variance in female and male individual effects are indicated by VI-female and VI-male; variance in female and male additive genetic effects by VA-female and VA-male; covariance and correlation between female and male additive genetic effects by COVAfemale,male and rA-female,male; and variance in female and male permanent individual effects by VPI-female and VPI-male, respectively. VYR is the among-year variance, and VR, the residual variance. Sample sizes (i.e. the number of observations) for models 1 and 3 were 1106, and 1035 for models 2 and 4. Terms not fitted are marked with "-".