Measures of the Focal Fixed Effects: Inbreeding, Age, and Early Nutrition

We used the pedigree-based inbreeding coefficient F_ped, calculated using the R package pedigree version 1.4 (Coster 2015), as a measure of the degree of inbreeding of an individual (Wright 1922; Knief et al. 2016b); F_ped reflects the proportion of an individual’s genome that is expected to be identical by descent. Hence, F_ped can be used to estimate without bias the slope of the regression of fitness over inbreeding (Howrigan et al. 2011; Knief et al. 2016b). For instance, full-sibling mating produces inbred offspring that are expected to have 25% of the genome identical by descent ($F_{\mathrm{ped}}=0.25$). For practical reasons, all founders were assumed to be unrelated ($F_{\mathrm{ped}}=0$; Forstmeier et al. 2004), even though their true level of identity by descent is likely about 5% (judging from runs of homozygosity; Knief et al. 2015a).

For all birds, we recorded their exact hatch date. Thus, for models of reproductive performance at the level of eggs, clutches, and breeding rounds (as the unit of analysis), we used the exact age (in days) of the female or the male when an egg was laid, a clutch started, or a breeding round started, respectively. At the start of reproduction, individuals were 69–2,909 days old (fig. S1).

On the day of hatching, we individually marked all nestlings on the back using waterproof marker pens (randomly using red, blue, and green and pairwise combinations of these colors if there were more than three nestlings). We checked survival almost daily (daily on weekdays, occasionally during weekends) until offspring became independent (age 35 days). As a measure of early-growth condition, we determined body mass of each nestling to the nearest 0.1 g at 8 days of age (hereafter, condition). Despite the fact that high-quality food was available to all parents ad lib., nestling body mass at this age ranged from about 1.5 to 12.6 g (mean = $7.1\pm 1.7$ SD). For 297 of 6,190 nestlings, body mass was measured on day 6, 7, or 9. For those individuals, we estimated their mass on day 8 as follows. We constructed a linear mixed effects model, with nestling body mass as the dependent variable, actual age of the mass measurement and F_ped as two continuous covariates, and population (1–4; see above) as a fixed factor. We also included the identity of the genetic mother as a random effect. Using the slope of daily mass gain, we estimated mass at day 8 for those 297 individuals by adding or subtracting 0.97 g per day of measuring too early or too late. Because the four populations differ in body mass, we normalized (Z scaled) all measured or estimated values of mass at day 8 within each population before further analysis.

We report effects of inbreeding, age, and early condition always with a negative sign, such that negative values of greater magnitude reflect stronger detrimental effects of being inbred, old, or poorly fed. This allows us to metasummarize the results and to directly compare the strength of the focal fixed effects on reproductive performance.

Measures of Life Span and Reproductive Performance Traits

Table 2 provides an overview of all traits included in this study. To allow direct comparison and easy interpretation of the fixed effects and additive genetic correlations, we scored all traits such that higher positive values reflect better reproductive performance.

Life span was analyzed in the following subset of birds: five generations of birds from the Seewiesen population (referred to as generations P, F₁–F₃, and S₃; $N=\mathrm{1,855}$ individuals) and four generations of birds from the Bielefeld population (F₁–F₄; $N=\mathrm{1,067}$ individuals). Among those birds, we used the four most complete generations, P and F₁–F₃ Seewiesen, for which we recorded the exact life span for all ($N=\mathrm{1,175}$ individuals) as a pool to impute missing life spans. For 219 S₃ Seewiesen birds and for 663 Bielefeld birds, no date of natural death was available (e.g., because individuals were still alive or because their fate was unknown). For these individuals, we used imputed life expectancy in all analyses, defined as the average life span of individuals from the same pool that lived longer than the focal bird when last observed alive.

In aviaries, we identified social pairs by behavior (clumping, allopreening, and visiting a nest together). All parentage assignments were based on conventional microsatellite genotyping using 10–15 microsatellite markers on up to 13 chromosomes (Wang et al. 2017), following Forstmeier et al. (2007a). We assigned every fertilized egg to its genetic mother ($N=\mathrm{11,704}$ eggs). When the egg appeared infertile (no visible embryo; Birkhead et al. 2008), we assigned it to the social female that was attending the clutch ($N=\mathrm{3,630}$ cases). In 36 cases where two females used the same nest to lay eggs, we assigned the unfertilized eggs to the female that laid the most similar eggs (in size and shape), based on eggs that were certainly laid by a given female (e.g., fertilized eggs and eggs in other clutches laid by that female). In cases where birds were not allowed to rear offspring, we quantified female fecundity as the total number of eggs laid by the focal female during the breeding period (see tables S1, S2).

In breeding experiments, we opened all unhatched eggs to check for visible signs of embryo development and classified them as either infertile or embryo mortality. In experiments in which all eggs were incubated artificially for a few days to collect DNA from embryos, we classified eggs as infertile or not but discarded information on embryo viability. Visual inspection of opened eggs has the disadvantage that early embryo mortality may get misclassified as infertility if it occurred before any visible signs of development. Misclassification cannot be avoided entirely, even with more time-consuming examination of eggs, which would be challenging to do for thousands of eggs (Birkhead et al. 2008; Murray et al. 2013). However, genotyping the germinal disk and counting sperm on the perivitelline membranes of 76 freshly laid eggs revealed 22 apparently infertile eggs. Only one of those (5%) had more than 20 sperm on the perivitelline membrane, suggesting early embryo mortality (fig. S2; see also Birkhead and Fletcher 1998). In contrast, among 37 eggs with more than 20 sperm on the perivitelline membrane, 36 (97%) developed diploid tissue. Thus, we expect only a small fraction of misclassification.

In cages, we measured male fertility as a binary trait, that is, whether an egg was fertilized. Because extrapair copulations can be excluded in cages, we only genotyped all surviving offspring with the same set of microsatellites used in aviaries as confirmation (Wang et al. 2017). In 12 cases, one to five eggs (median: one egg) were fertilized by the previous partner of the female, and those were counted as infertile eggs of the focal male. In aviaries, we assessed for each egg whether it was sired by the social male of the female who laid the egg. We refer to this as male within-pair paternity, a trait that reflects a male’s ability to defend his paternity against extrapair males. We also quantified male siring success as the total number of fertilized eggs sired by a focal male. This includes males that remained unpaired (without a social female).

For each fertilized egg that was incubated by the social parents, we recorded whether it hatched (binomial trait for the genetic parents). For each hatched egg that was reared, we recorded whether the nestling survived to independence (day 35; binomial trait for the social parents). We quantified the number of seasonal recruits as the number of genetic offspring that survived to independence within a given breeding season. The number of seasonal recruits was square root transformed to approach normality.

Statistical Models

All mixed effects models were run in R version 4.0.0 (R Core Team 2020), using the R package lme4 version 1.1-23 (Bates et al. 2015). All animal models were run using VCE6 (Neumaier and Groeneveld 1998) because (a) it allows running a 12-trait multivariate animal model that consists of 2,346 individuals with at least one trait value per individual and (b) it has a reasonable running time. To check the consistency of model outputs, we repeated all animal models in the R packages pedigreemm version 0.3-3 (Vazquez et al. 2010; univariate animal models only) and MCMCglmm (Hadfield 2010; univariate and bivariate animal models). All model details, with the supporting data and R scripts, have been deposited in the Open Science Framework (https://doi.org/10.17605/OSF.IO/TGSZ8; Pei 2020). Model outputs of all methods are given in the supplemental PDF. The heritability and additive genetic correlation estimates were highly correlated between methods ($r>0.65$, $P<.002$). We report the VCE6 estimates, unless otherwise stated. Figure S3 shows the exact range of each focal fixed effect and each performance trait value. Here, we Z transformed all covariates and response variables across populations to allow direct comparison of the effect sizes for inbreeding, age, and condition across all models. The 95% confidence intervals of fixed effects from mixed effects models were calculated using the function glht from the R package multcomp version 1.4-13 while controlling for multiple testing (Hothorn et al. 2008). Data analysis involved four consecutive steps (fig. 1).

Step 1: estimation of fixed effects and variance decomposition. The goal of step 1 was to estimate (a) all fixed effects on reproductive performance and (b) individual repeatability of performance traits (fig. 1). All fixed and random effects of models used in step 1 are listed in tables S3 and S4. In brief, we first fitted all models with a Gaussian error distribution to compare and metasummarize the estimated effect sizes of the fixed effects and to estimate the variance components for the random effects. We used all observations with information on the three fixed effects (age, F_ped, and early condition of the male, female, and the individual egg if applicable) and included population (fixed effect) and female, male, and pair identity (random effects). We analyzed traits that were measured at either egg, clutch, or season level. As applicable, we fitted as fixed effects the laying sequence of eggs within a clutch, the order of hatching of offspring within a brood, the order of the clutches that were laid by a female over the course of a season, the sex ratio in the aviary, and the duration of the season (table S1). For models of embryo survival, we also controlled for whether the eggs were incubated in a nest that still contained offspring from a previous brood (7% of embryos). For models of nestling survival, we added as fixed effect pair type (pair formed through mate choice or through force pairing; Ihle et al. 2015). For models of egg-based fertility, within-pair paternity, and embryo and nestling survival, we also tested the effect of egg volume on egg fate (we calculated volume as V = (1/6) × π × width² × length, where egg length and width had been measured to the nearest 0.1 mm). For this analysis, we fitted the mean egg volume of each female and the centered egg volumes (centered within individual females) to distinguish between the effects of between- and within-female variation in egg size (van de Pol and Wright 2009). We estimated the variance components for male, female, and pair identity and further controlled for clutch identity and identity of the setup (see tables S1, S2), as applicable, by adding them as random effects. Life span had no repeated measurement; therefore, we included only individual identity as a dummy random effect for practical reasons when running the model and extracting estimates in R. For this lm model, the correlation between the residuals and the dummy random effect equals 1, and the fixed effect estimates were unaffected by the dummy variable. Table 2 shows for which group of individuals, that is, female, male, or the offspring itself, we tested which focal fixed and random effects.

To allow direct comparison of the magnitude of fixed effects at the same level of measurement, we also aggregated data within clutches (e.g., proportion of infertile eggs within a clutch) and within individuals over the course of a season. Models on aggregated data were weighted by the number of eggs within a clutch or by the number of eggs or clutches for an individual within a season (fig. 1). As expected, the proportion of variance explained by male, female, and pair identity increased from the egg level to the season level (see “Results”). However, the relative proportions explained by female, male, and pair identity did not change notably. Therefore, we focus on the analyses of fixed effect estimates at the breeding season level.

To compare the overall effect sizes between the focal fixed effects, we metasummarized the estimated effect sizes for inbreeding, age, and condition using the weighted lmer function from the R package lme4 (fig. 1, step 1, metasummarization of estimated effect sizes). The uncertainty of each estimate was accounted for by using the multiplicative inverse of the standard error (1/SE) of the response variable as weight. In this metamodel, we used effect size estimates from models that had been aggregated at the season level as the dependent variable. Note that effects of inbreeding of the egg on fertility in cage breeding and nestling survival were taken from egg-based models because they cannot be aggregated by clutch or season. Additionally, we tested whether effect sizes differed among males, females, and offspring (fixed effect with three levels) or among traits (random effect with 11 levels; as listed in Table 2).

Additionally, we tested for early-starting aging effects by selecting reproductive performance data for males and females that were <2 years old when reproducing. We then metasummarized the mean age effect estimates using the R function lm, weighted by the multiplicative inverse of the standard error.

We calculated the amount of variance explained by each fixed effect (Nakagawa and Schielzeth 2010) as the sum of squares of the fixed effect divided by the number of observations ($N-1$; Henderson 1953). In weighted models, we divided the variance components of the fixed effects and the residual by the mean weight value (Bates et al. 2015).

Step 2: estimation of heritability of fitness-related traits. The goal of step 2 was to estimate the heritability of reproductive performance traits using univariate Gaussian animal models (fig. 1). Because quantitative genetic models require large amounts of data, we restrict our analyses to the populations Seewiesen and Bielefeld. Note that the pedigrees of our four captive populations are not connected, so it was not useful to analyze them jointly.

We kept the general model structure from step 1 but excluded the fixed effects of egg volume on male fertility, embryo, and offspring survival (to avoid removing biological variation that is potentially heritable and hence of interest; note that the effect sizes of egg volume are small; see “Results”). For the embryo survival model, we excluded the nonsignificant fixed effects of male age, inbreeding, and condition. For the model on male fertility from cage breeding, we excluded the nonsignificant effect of the level of inbreeding of the egg itself. To most effectively use the available information on reproductive performance, we included individuals with missing values for condition ($N=231$ founder individuals and $N=23$ individuals of the F₂ generation; i.e., 7% of Seewiesen birds). These missing values were replaced by the population mean. Individual identity was fitted twice, once linked to the individual correlation matrix (pedigree) to estimate the amount of variance from additive genetic effects (V_A) and once to estimate the remaining amount of variance from permanent environmental effects (V_PE; Kruuk and Hadfield 2007). Animal models on nestling mortality were run twice, once for the mother and once for the father. We calculated heritability based on the total phenotypic variance, V_Ph, as $h^2=\left(V_{\mathrm{A}}/V_{\mathrm{Ph}}\right)$, and we also quantified V_A relative to individual repeatability as ($V_{\mathrm{A}}/\left(V_{\mathrm{A}}+V_{\mathrm{PE}}\right)$).

We compared the estimates of heritability (and V_A relative to the individual repeatability) between the domesticated population Seewiesen and the recently wild-derived population Bielefeld using the R function lmer. We used the multiplicative inverse of the standard error as weight to control for variation in uncertainty of each estimate. We used the estimates of heritability as the response variable and fitted population as a fixed effect (two levels) and trait as a random effect (nine levels, including only traits that were measured in both populations).

Step 3: calculation of mean individual fitness-related trait values using best linear unbiased predictions (BLUPs). The only goal of step 3 was to extract individual estimates of reproductive performance needed for step 4. We kept the model structure from step 1, except that we used a binomial error structure for binary traits, that is, male fertility in cages and aviaries and embryo and nestling survival. Missing values for condition (mostly founders of each population; 6% of all birds of the four populations) were replaced with population means as in step 2. For the embryo survival model, we again excluded the nonsignificant effects of male inbreeding, age, and condition. We also excluded (a) effects of egg volume from all egg-based models and (b) the effect of the level of inbreeding of the egg itself from the model of male fertility measured in cages (see step 2). We extracted the BLUPs for female or male identity (as applicable) as the estimated life-history trait value of that individual (Table 2) for step 4.

Step 4: estimation of additive genetic correlations. The goal of step 4 was to estimate additive genetic correlations between different performance traits using multivariate animal models. Before fitting a 12-trait animal model that estimates for each matrix (genetic and residual) all 12 variances and 66 covariances simultaneously, we aggregated the raw data to one phenotypic value per individual for each trait (fig. 1, step 3). This was necessary because we are not aware of software that can handle the full complexity of the underlying raw data (involving more than 26 different fixed effects). Because simple averages of multiple measures can result in outliers when sample size is small, we used the phenotypic BLUPs described above. BLUPs do not produce outliers and account for all considered fixed and random effects (Robinson 1991; Houslay and Wilson 2017). Breeding values (genetic BLUPs) suffer from nonindependence because the phenotype of one individual influences the breeding values of all its relatives (Hadfield et al. 2010). Note that this is not the case for the phenotypic BLUPs we use here. However, the uncertainty that is inherent to each BLUP is not taken into account, which may lead to underestimation of standard errors (Houslay and Wilson 2017). To check the robustness of our results, we compared our estimates with those obtained (a) using a smaller data set from another population (Bielefeld) with the same method and (b) using bivariate animal models in MCMCglmm version 2.29 (Hadfield 2010; population Seewiesen). The latter approach is presumably less powerful than a full 12-trait animal model.

For each of the 12 traits, we fitted an intercept and the pedigree as the only random effect to separate additive genetic from residual variance. We ran these models for the largest and most comprehensive data set (population Seewiesen; $N=\mathrm{2,346}$ individuals with at least one trait value, BLUPs for 12 traits, and 66 covariances) and for the more limited data set (population Bielefeld; $N=\mathrm{1,134}$ individuals, BLUPs for 9 traits, and 36 covariances; see “Results”; fig. 1, step 4, estimate additive genetic correlations).

We used the weighted lm function in the R package stats to summarize the estimated additive genetic correlations within and between the major categories of traits, that is, female, male, offspring traits, and life span, for each population separately (Table 2; fig. 1, step 4, metasummarization of estimated additive genetic correlations). We fitted the estimates of additive genetic correlations (for each pair of traits, weighted by the multiplicative inverse of the standard error of each estimate) as the dependent variable, with trait class combination as a predictor with seven levels. We removed the intercept to estimate the mean additive genetic correlation for each pairwise combination of classes. We then computed the eigenvectors of the additive genetic variance-covariance matrix of traits, using the R function eigen, and visualized the orientation of the traits in the additive genetic variation space defined by the principle components PC1 and PC2 (Berner 2012). The proportion of variance explained by the first two principle components was calculated using the functions summary and prcomp in the R packages base and stats, respectively.

Effects of Laying and Hatching Order, Clutch Order, and Egg Volume on Egg and Embryo Fate

The fate of an egg and its embryo depended on the order of laying within a clutch, the order of hatching within a brood, and the order of consecutive clutches within a breeding season (fig. S4; table S3, models at the egg level; see also fig. 1, step 1). First-laid eggs in a clutch were significantly more likely to be infertile or to contain a dead embryo. Fertility and embryo viability were the highest for the third egg (fig. S4). Male fertility significantly increased over the first three clutches and stayed high afterward. In contrast, clutch order did not affect the probability of embryo and nestling survival.

The average effect of egg volume on measures of egg fate was small (mean: $r=0.040\pm 0.016$ SE; fig. S5). Effects of egg volume were largest for nestling survival after hatching and smallest for embryo survival (table S3; fig. S5). Despite large sample size ($N=\mathrm{9,785}$ eggs), embryo survival was not significantly influenced by egg volume (between-female variation: $r=0.015\pm 0.017$ SE, $P=.37$; within-female variation: $r=0.018\pm 0.010$ SE, $P=.08$; table S3). Additionally, embryos in clutches that were incubated in the presence of nestlings from previous breeding attempts were more likely to die before hatching ($b=0.192\pm 0.048$ SE, $P<.0001$; table S3). Overall, the total amount of variance explained by laying and hatching order, clutch order, and egg volume on egg fate was less than 5% (table S4).

Effects of Inbreeding, Age, and Early Condition

Individuals performed worse in virtually all studied reproductive traits when they were more inbred, as they became older, and when they weighed less at 8 days of age (figs. 2, S3, S6; table S3; see also fig. 1, step 1). Interestingly, reproductive performance did not show an initial increase at a young age (metasummarized effect size of age among birds younger than 2 years: $r=\mathrm{-0.013}\pm 0.011$ SE; figs. 2C, 2F, 3, A3). Inbred eggs were equally as likely to be infertile as outbred eggs, while inbred embryos and offspring were more likely to die (fig. 3C). Together, this suggests that most infertile eggs were not cases of undetected early embryo mortality. Individuals lived shorter lives when they were inbred and when they had low weight at day 8 (fig. 3; table S3). However, the fixed effects of inbreeding, age, and condition together explained, on average, only 2% of the variance across all traits (fig. 4; table S5).

Metasummarized effect sizes of inbreeding ($r=\mathrm{-0.117}\pm 0.024$ SE) and age ($r=\mathrm{-0.132}\pm 0.032$ SE) were similar in magnitude and were about twice as large as the remarkably small effect of early condition ($r=\mathrm{-0.058}\pm 0.029$ SE; fig. 3; table S4; see also fig. 1, step 1, metasummarization of estimated fixed effects). There was no significant difference among males, females, and offspring in how strongly they were affected by these three factors ($b\le 0.012\pm 0.028$ SE, $P=.63$; table S4). Fitting trait (fitness component, 11 levels) as a random effect explained 1.5% of the variance in effect sizes ($P=.02$; table S4), suggesting that some components might be less sensitive than others (fig. 3; table S3). Female traits significantly predicted offspring survival and male fertility (independent of whether they were measured in a cage or in an aviary), whereas male traits showed no effect on offspring survival (fig. 3).

Variance Components and Heritability

Variance components for all reproductive performance traits are shown in figure 4 (see also table S4; fig. 1, step 1, estimate repeatability). Overall, individual reproductive performance traits were significantly repeatable (median $R=0.28$, range: 0.15–0.50). Female reproductive performance traits (clutch size, fecundity, and female seasonal recruits) showed reasonably high repeatability for individual females ($R\sim 0.26–0.40$). Likewise, male fertility, male siring success, and male seasonal recruits were highly repeatable for individual males ($R\sim 0.24–0.50$). Female reproductive traits from aviary breeding were analyzed independently of whether the focal female had a partner (Table 2), but female clutch size measured in a cage showed no contribution from the male partner or from pair identity. In contrast, male fertility depended on all three random effects and was repeatable for males ($R>0.23$, $P<.0001$) but less so for females ($R<0.18$, $P<.1$) or for the particular pair combinations ($R<0.23$, $P<.05$). The model on embryo survival showed significant female and pair identity (genetic parents) effects that were similar in size (both $R=0.20$, $P<.0002$), while genetic male identity explained no variance (fig. 4). In contrast, social female ($R=0.17$, $P=.017$) and social male ($R=0.15$, $P=.039$) identity explained significant amounts of the variance in nestling survival, while the effect of pair identity (parents that raised the brood) was less clear ($R=0.14$, $P=.11$).

Reproductive performance traits and life span in general had low narrow-sense heritability ($V_{\mathrm{A}}/V_{\mathrm{Ph}}$; Seewiesen: median $h^2=0.07$; Bielefeld: median $h^2=0.11$) and explained only a limited amount of the individual repeatability ($V_{\mathrm{A}}/(V_{\mathrm{A}}+V_{\mathrm{PE}})$; Seewiesen: median = 0.29; Bielefeld: median = 0.32; see all heritability estimates in tables S6 and S7; fig. 1, step 2). Heritability estimates from the recently wild-derived population Bielefeld were similar to those from the domesticated Seewiesen population (for nine traits measured in both populations; mean difference in $h^2=0.02$, range: −0.10 to 0.13, metasummarized difference after controlling for the uncertainty of each estimate: $\mathrm{\Delta }b<0.0001$; mean difference in $V_{\mathrm{A}}/\left(V_{\mathrm{A}}+V_{\mathrm{PE}}\right)=0.20$, range: −0.13 to 0.68, metasummarized difference: $\mathrm{\Delta }b=0.0002$; table S8).

Additive Genetic Correlations

Reproductive performance traits were grouped into three classes: (1) aspects of male reproductive performance, (2) aspects of female reproductive performance, and (3) aspects of offspring survival (Table 2). Traits within each of these classes were, on average, positively correlated with each other at the additive genetic level (for the Seewiesen population, female traits: mean $r_{\mathrm{A}}=0.66$, $P<.0001$; male traits: mean $r_{\mathrm{A}}=0.67$, $P<.0001$; offspring survival traits: mean $r_{\mathrm{A}}=0.36$, $P=.09$; fig. 5A; see also fig. 1, steps 3 and 4). Results for the Bielefeld population are shown in figure S7. The metasummarized results are given in table S9, and all additive genetic correlation estimates are listed in tables S10 and S11 (fig. 1, step 4). Estimates of the additive genetic correlations from bivariate animal models using MCMCglmm are shown in figures S8 (Seewiesen) and S9 (Bielefeld).

Male and female reproductive performance traits were weakly negatively correlated at the additive genetic level (mean $r_{\mathrm{A}}=\mathrm{-0.14}$, $P=.04$; see MF in figs. 5A, S8A). Accordingly, the eigenvectors for male and female fitness traits were pointing in different directions (figs. 5B, S8B). This pattern was somewhat consistent between the Seewiesen and Bielefeld populations (see figs. S7 and S9 for the Bielefeld population). However, the negative correlation between male and female fitness traits was no longer significant when estimated by the bivariate animal models in MCMCglmm and disappeared in the Bielefeld data set (table S9). The orientation of offspring survival traits relative to male and female fitness traits was less consistent. In the Seewiesen population, female fitness traits were positively correlated with offspring survival traits at the additive genetic level (mean $r_{\mathrm{A}}=0.61$, $P<.0001$), while male fitness traits were not aligned with offspring survival traits (mean $r_{\mathrm{A}}=\mathrm{-0.11}$, $P=.24$; fig. 5). In contrast, in the Bielefeld population, both female and male fitness traits were positively correlated with offspring survival traits (fig. S7). Life span tended to be positively correlated with all reproductive performance traits (Seewiesen: mean $r_{\mathrm{A}}=0.19$, $P=.02$; Bielefeld: mean $r_{\mathrm{A}}=0.60$, $P=.0006$; figs. 5, S7; table S9).

Effects of Inbreeding, Age, and Early Condition

Many studies have shown that inbreeding depression significantly influences morphological, behavioral, and fitness-related traits in zebra finches (Bolund et al. 2010a; Forstmeier et al. 2012; Hemmings et al. 2012; Opatová et al. 2016) and in other species (Amos et al. 2001; Reed and Frankham 2003; Williams et al. 2003; Michaelides et al. 2016). This study confirms that inbreeding negatively influenced all phases of offspring survival, reproductive performance, and life span. We found that the level of inbreeding of both genetic parents negatively influenced egg fertility, suggesting that this is a matter of not only sperm functionality (Opatová et al. 2016) but also female reproductive performance (e.g., egg quality or copulation behavior). Male and female fitness estimates (seasonal recruits) were most strongly affected by inbreeding (fig. 3), presumably because the successful rearing of offspring to independence requires proper functionality at every step of reproduction.

Age effects on reproductive performance typically show an initial increase in performance in both short- and long-lived species (e.g., in great tits Parus major [Bouwhuis et al. 2009], wandering albatrosses Diomedea exulans [Lecomte et al. 2010], Houbara bustards Chlamydotis undulata [Preston et al. 2011], Langur monkeys Presbytk entellus [Harely 1990], and red deer Cervus elaphus [Pemberton et al. 2009]). Interestingly, in our captive zebra finches we found that reproductive performance (especially male fertility, female clutch size, fecundity, and female effects on embryo survival) did not show an initial increase after birds reached sexual maturity at about 100 days of age (figs. 2C, 2F, A3). This could be because zebra finches are short-lived opportunistic breeders that reach sexual maturity earlier compared to most other birds (Zann 1996). Thus, zebra finches might have been selected to perform best early on. Alternatively, this effect may not be present in the wild, where experience might play a more important role in determining reproductive success.

Over the past decades, numerous studies focused on how early developmental conditions affect behavior, life history, and reproductive performance later in life (Tschirren et al. 2009; Rickard et al. 2010; Boersma et al. 2014). Here we show that even dramatic differences in early growth conditions of surviving offspring (see range of X-axis in fig. 2B) have remarkably small (though statistically significant) effects on adult reproductive performance.

Overall, the proportion of variance explained by inbreeding, age, and early condition (characteristics of conditions) was less than 3% (fig. 4; table S4). This indicates that individuals’ robustness against poor conditions appears more noteworthy than their sensitivity. As will be discussed below, the majority of the individual repeatability in reproductive performance cannot be explained by such individual characteristics.

Repeatability and Heritability of Reproductive Performance

Individual zebra finches were remarkably repeatable in their reproductive performance. Our variance-partitioning analysis showed that infertility is largely a male-specific trait, whereas embryo and offspring survival are mostly related to female identity (fig. 4; table S4). In contrast, in polyandrous crickets, egg hatching (primarily a matter of embryo survival) was mostly influenced by male identity (García-González and Simmons 2005; Ivy 2007). The effects of pair identity on infertility and offspring mortality in zebra finches may reflect behavioral incompatibility, while the pair effect on embryo mortality more likely reflects genetic incompatibility (Ihle et al. 2015).

Although male and female zebra finches are highly repeatable in their reproductive performance, the heritability of fitness traits was low. Heritability estimates were similar between the recently wild-derived Bielefeld population and the domesticated Seewiesen population. This contradicts the idea that ongoing adaptation to captivity would result in a higher heritability of fitness traits. Overall, our findings indicate that there are some additive genetic components underlying zebra finch reproductive performance.

Evidence for Sexually Antagonistic Pleiotropy and Other Potential Causes of Reproductive Failure

Some of the standing additive genetic variance in reproductive performance could be maintained by intralocus sexual antagonism between male fitness traits and female (and offspring) fitness traits (Cox and Calsbeek 2009). This has, for example, been suggested in quantitative genetic studies on Drosophila (Innocenti and Morrow 2010), red deer (Foerster et al. 2007), and the bank vole Myodes glareolus (Mills et al. 2012). We found that male fertility, siring success, and seasonal recruitment were overall negatively correlated with female fitness and offspring survival traits, suggesting that alleles that increase male fitness tend to reduce female and offspring fitness (fig. 5). In contrast, life span and reproductive performance tended to be positively correlated at the additive genetic level, which is suggestive of some overall good gene variation in our population (fig. 5). Some words of caution should be added to these observations. VCE6 (figs. 5, S7) yielded higher absolute values of estimates than those calculated with the R functions PedigreeMM (heritability estimates only) and MCMCglmm (see figs. S8, S9; also see tables S6, S7, S10, S11). Nevertheless, the additive genetic correlation estimates are highly correlated between the two methods ($r>0.7$, $P<.0001$; see tables S10, S11). Estimating genetic correlations between traits with low heritability requires large data sets, especially on additive genetic correlations of between-sex reproductive performance where the traits of male fertility and female clutch size in cages are missing (N performance traits: Seewiesen = 12, Bielefeld = 9; N birds have at least one entry of reproductive performance data: Seewiesen = 2,346, Bielefeld = 1,134; hence, these results are presented in fig. S7). Despite this lack of power in our second-largest data set of population Bielefeld, its overall orientation of traits in the additive genetic variation space of the principle components PC1 and PC2 is very similar to population Seewiesen (note that life span is in the center of all fitness traits and that aspects of female fitness do not align with male fitness in figs. 5B, S7B, S8B, and S9B).

Individual repeatability of fitness-related traits could arise from permanent environmental effects (e.g., early developmental conditions and long-lasting diseases) or from genetic effects. However, although food shortage experienced during early development (reflected in body mass at 8 days old) strongly predicted nestling mortality (Pei et al. 2020), it explained only <1% of variation in reproductive performance later in life (mean $r=\mathrm{-0.058}$; figs. 2B, 2E, 3, 4). Additionally, our captive zebra finches were raised and kept in a controlled environment with no obvious diseases detected. Additive genetic effects explained only about 30% of the large remaining unexplained individual repeatability in fitness-related traits, suggesting that reproductive performance might be (predominantly) dependent on genetic effects of local over- or underdominance and epistasis, that is, incompatibility between loci. For instance, high levels of reproductive failure could be maintained when alleles show nonadditive effects, with selection favoring the heterozygous genotype (see, e.g., Sims et al. 1984; Grossen et al. 2012). In zebra finches, males that are heterozygous for the inversion on the Z chromosome produced fast-swimming sperm and sired more offspring (Kim et al. 2017; Knief et al. 2017), while heterozygous males for inversions on chromosomes Z and 13 produced slightly more dead embryos, likely caused by unbalanced crossover during spermatogenesis (Knief et al. 2016a). However, these phenomena explain only a small fraction of infertility and embryo mortality. Overall, there is little over- or underdominance for fitness related to the major inversion polymorphisms that segregate in wild and captive zebra finch populations (Knief et al. 2016a).

Epistatic effects that involve several genes (e.g., incompatibility between nuclear loci or between mitochondrial and nuclear genomes; Zeh and Zeh 2005) could be evolutionarily stable when certain combinations of genotypes perform better than others, especially when combined with overdominance. Examples of incompatibilities are mostly known from hybrid systems (Arntzen et al. 2009; Hermansen et al. 2014; Eroukhmanoff et al. 2016), but they could also be segregating within a species after the mixing of two lineages that have evolved weak incompatibilities. Some studies on inbred lines in invertebrates found evidence of mitonuclear incompatibilities. For example, in the spider mite Tetranychus evansi, the eggs of F₁ hybrid females of two genetic lineages showed higher hatching failure compared to the pure parental lines (Knegt et al. 2017), and in Drosophila melanogaster, the interaction of mitonuclear background explained a small but significant amount of variation in female fitness (Dowling et al. 2007).

Infertility, as one of the main and puzzling sources of reproductive failure, behaved as a male-specific trait that may also depend in part on behavioral compatibility between pair members (reflected in copulation behavior) and in part on the male’s genotype at sexually antagonistic loci. The intrinsic male fertility, measured in a cage, that is, in the absence of sperm competition, correlated negatively with all female and offspring survival traits at the additive genetic level (sexual antagonism; median $r_{\mathrm{A}}=\mathrm{-0.30}$, range: −0.45 to −0.01; fig. 5; table S10). In contrast, in the presence of sperm competition (aviary breeding), high male within-pair paternity, siring success, and seasonal recruitment should also be influenced by the competitive ability of the individual, and this could explain why these traits correlated positively with life span and trade off less with female traits and offspring rearing ability at the additive genetic level (figs. 4, S7; tables S10, S11).

Embryo mortality, another main source of reproductive wastage, mostly depended on the identity of the genetic mother and the identity of the genetic pair members. A previous study using cross fostering of freshly laid eggs also showed that embryo mortality is a matter of the genetic parents rather than the foster environment (Ihle et al. 2015). The female component suggested an overall female genetic quality effect yet with limited heritability (pointing toward dominance variance or epistasis). The effect of the combination of parents on embryo mortality might reflect an effect of the genotype of the embryo itself, possibly involving multilocus incompatibilities (Corbett-Detig et al. 2013).