Selection of lamb size and early pregnancy in Soay sheep (Ovies aries)

The paradox of stasis – the unexpectedly slow evolution of heritable traits under direct selection – has been widely documented in the last few decades. This paradox is often particularly acute for body size, which is often heritable and where positive associations of size and fitness are frequently identified, but constraints to the evolution of larger body sizes are often not obvious. Here, we identify a trade-off between survival and size-dependent reproduction in Soay sheep (Ovis aries), contributes to selection against large body size. Using recently developed theory on non-linear developmental systems, then decompose total selection of ewe lamb mass along different causal paths to fitness. Larger lambs are more likely to become pregnant, which has a large viability cost. After controlling for this pathway, however, the association between lamb mass and subsequent lifetime fitness is positive. Thus this trade-off does not fully explain stasis of size in tis population, but it does substantially reduce the strength of positive directional selection of size that would otherwise occur. While selection currently favours reduced probability of early pregnancy, largely irrespective of body size, it is likely that the occurrence of early pregnancy could result from adaptation to conditions during a recent period during which population density was much lower.


17
The paradox of stasis (Sogard, 1997) and the pattern of bigger-is-better with respect to selection on 18 body size have received great attention in the last few decades from both macro-and microevolutionary biologists. Positive directional selection on size appears to be widespread (e.g. Kingsolver & 20 Pemberton, 2004a). More than 95% of the individuals living in the study area have been marked with plastic ear tags shortly after birth to enable identification throughout their lifetimes, at which 90 time maternal identity and twin status is ascertained from field observations. Regardless of whether 91 a ewe's lamb is captured (e.g., in cases of neonatal mortality), reproductive status is extensively 92 recorded, allowing assessment of annual and lifetime fitness of ewes. Regular censuses and mortality 93 searches provide death dates for the majority of the sheep. Each year in August a large portion of 94 the study population is captured and phenotyped for multiple traits, including body mass. We focus 95 on phenotypic and life history data on ewe lambs from 1991 to 2015, corresponding to 3457 ewes, as 96 there is no systematic record of pregnancy status from post-mortems of ewes that died during the 97 winter before 1991.

99
The pedigree of this population has been constructed through a combination of observational field 100 data and molecular markers for maternal links, and using molecular markers only for paternal links 101 (Johnston et al., 2013;Bérénos et al., 2014). 315 polymorphic and unlinked SNP markers were used 102 in molecular parentage assignments (for 4371 individuals) with 100% confidence in the R package A linear regression of the form 2) + β β βX ijk X ijk X ijk + a i + u cj + u mk + ijk (1) was adopted to model August body mass of lamb i, born in year j to mother k. α and α t are the 157 intercept and the twinning-specific contrast to the intercept, β d β d β d are the linear (slope) and quadratic 158 terms associated to a second order polynomial function on population density, d, and β β β is a vector 159 with the coefficients for the remaining fixed effects, X X X (see Appendix B). Random intercepts to es-160 timate among-mother (u mk ) and among-cohort (u cj ) variation were also included. Finally, a i is the  Table B1.

164
Early pregnancy 165 The probablibility of early pregnancy was estimated using a binomial regression with a logistic link 166 function, where β m β m β m corresponds to a vector with coeffcients for the linear, quadratic, and cubic terms for lamb pregnancy, p. Although the probability of lamb pregnancy is well described by a binomial distribution, a latent scale variable, p l , can be defined such that p l = ln p 1−p . In that case, both 182 m and p l are assumed to be drawn from a multivariate normal distribution, with a mean vector 183 that includes both mean mass, µ m , and mean pregnancy probability in the logit scale, µ p l , and a 184 variance-covariance matrix P l P l P l , Both µ m and µ p l depend on twin status (singleton or twin), population density, and maternal age 186 at parity, whereas the former also depends on birth and measurement dates. The corresponding 187 parameters and their estimates are listed in Appendix B. Note that the subscript l denotes latent 188 scale, which here applies to early pregnancy only, since August lamb body mass was modelled on 189 its natural scale. The genetic and environmental contributions to P l P l P l were partitioned by including 190 random effects on breeding values (animal model, Henderson 1975). In fact, a G l G l G l matrix and an E l E l E l 191 matrix can be defined such that where subscripts a and e denote the additive genetic and the environmental contribution to the over-193 all (co)variances, respectively. For more information on (co)variance partition, see Falconer (1981).

194
To simplify the notation, we do not use any extra subscript to denote phenotypic (co)variances (see 195 matrix P l P l P l in Equation 4). Besides residual variances and covariances, the matrix E l E l E l also includes 196 variances and covariances associated with cohort and maternal random effects. Since the overdis-197 persion variance of a generalised linear mixed model (GLMM) is unobservable for binary data, the 198 residual variance for p l was set to one. to expected and data scales, the difference between them being the random noise that distinguishes 215 a model from the data itself. Here, we adopt the expected scale. Following their work, the expected 216 values for body mass and the probability of early pregnancy is given by where g g g −1 corresponds to the inverse link functions used to model body mass and the probability 218 of early pregnancy. As the adopted link functions for these variables were the identity and the logit 219 functions, respectively, their inverses are the identity and the logistic functions, the latter correspond-220 ing to exp(p l ) exp(p l )+1 . l l l are the latent values for lamb body mass and early pregnancy. f MVN (l l l, µ µ µ, P l P l P l ) is 221 the probability density of a multivariate normal distribution with vector mean µ µ µ and variance P l P l P l .

222
Likewise, the phenotypic variance-covariance matrix on the expected scale is given by 223 P P P = g g g −1 (l) (l) (l) −z z z 2 f MVN (l l l, µ µ µ, P l P l P l ) dl.
The G G G matrix on the expected scale, can also be derived (Villemereuil et al., 2016), where Φ Φ Φ is the average derivative of the expected values 225 with respect to the latent values, Details on how derivatives and integrals were solved are provided in Appendix C. The heritability 227 on the expected scale is obtained directly through its definition, as the proportion of the additive 228 genetic variance relative to the phenotypic variance, using the derived values. The results obtained 229 using the above expressions can also be found in Table B3. The heritability of early pregnancy on 230 the expected scale (probability) is 0.40 (95% CrI 0.11; 0.63), higher than the approximation shown 231 above (but note that this approximation is derived for the data scale), and corroborating the genetic 232 basis of this trait. Using the derived values on the expected scale, we also obtained the conditional 233 genetic variance of lamb body mass. This measure allows us to evaluate the proportion of the additive 234 genetic variance in lamb body mass that is independent of early pregnancy. Following Hansen et al.

235
(2003), this quantity, σ 2 m|pa , is given by which in this particular case equals 0.40 (95% CrI 0.06; 0.85). The analogous metric for early preg-237 nancy corresponds to 0.03 (95% CrI 0.00; 0.05). Since σ 2 ma and σ 2 pa were estimated to be, respectively, 238 0.78 (95% CrI 0.20; 1.27) and 0.04 (95% CrI 0.00; 0.08), the additive genetic variance in lamb body 239 mass and early pregnancy that are independent from one another corresponds to 59% (95% CrI 29%; 240 100%). Substantial proportions of the genetic variances in lamb body mass and early pregnancy are 241 not independent, and therefore selection for each one of these traits necessarily results in selection of 242 the other.

244
First-year survival 245 We investigated whether early pregnancy has fitness costs by evaluating first-year survival in ewe 246 lambs as a function of pregnancy and body mass. First-year survival was defined using the first day 247 of May as the cut-off date to make sure that dying during labour was considered in the first annual 248 cycle of the ewes (89% of births are before 1 May). Ewes with known first-year survival status and 249 that survived until the rut were included in the analysis (n 2 = 947). A binomial regression with logit 250 link function of the form was adopted to model the probability of lamb i, born in year j to mother k, survives its first 252 year of life. α and α p are the intercept and the pregnancy-specific contrast to the intercept, β m 253 is the slope for lamb body mass, and β β β is a vector with the coefficients for the remaining fixed ef-254 fects. Although not pictured in Equation (10), a pregnancy-specific contrast to the slope for mass 255 is also estimated (Appendix B). As first-year survival, and survival in general, is highly dependent 256 on cohort effects due to variability in winter conditions and food availability (Clutton-Brock et al., 2004), we included cohorts as random effects (u cj ), as well as maternal environmental effects (u mk ).

258
Variance in residuals ( ijk ) was set to one, as overdispersion is unobservable in binomial mixed models.

260
We first establish that first-year survival is size-dependent by adopting a similar model to the one 261 in Equation (10) where α is the model intercept, β m is a slope for lamb body mass, and β β β is a vector with the coeffi-288 cients associated to the covariates in X X X, which includes capture and birth dates, as well a linear and 289 a quadratic terms for population density. The variance in the model residuals, i , was set to one.

290
Due to sample size limitations, we did not consider any random effects.

291
Only 2.2% (95% CrI 0.57%; 6.17%) of the offspring born to ewes getting pregnant as average-sized 293 lambs in average-density years survive until the winter (App. B, Tab. B5). It is interesting to note 294 that the slope for body mass is positive and significantly different from zero, showing that larger ewe 295 lambs are more likely to have surviving offspring, and suggesting that, as for first-year survival, if a 296 ewe lamb is to get pregnant, than the larger the better. Overall, AReS not only increases with ewe 297 body mass, but also significantly decreases as population size increases (App. A, Fig. A1c).

299
Subsequent lifetime rearing success (LReS) 300 We also investigated the effect of lamb body mass and early pregnancy on LReS. As the LReS of ewes where µ is the model intercept and µ t and µ p are the twin status and early pregnancy specific con-305 trasts to the intercept, respectively. β m is the slope for lamb August mass, β d β d β d the linear (slope) and 306 quadratic terms for density, and β β β a vector with the slopes of the remaining fixed effects, X X X. Although 307 not shown, a pregnancy-specific contrast to the slope for mass is also estimated (Appendix B). u c u c u c 308 and u m u m u m are random effects, assumed to be drawn from independent normal distributions, for cohort 309 and maternal identity. Finally, ijk is the residual of individual i. To study LReS, the data were 310 further constrained to only include ewes from cohorts that were completely phenotyped for this trait, including very few LReS records associated with ewes that were still alive by the cut-off year (2015).

313
While these surviving ewes' lifetime fitness is underestimated, they have reared most the offspring 314 they will rear in their lives. Excluding these records would result in stronger bias because ewes that 315 lived longer had potentially higher rearing success.

317
We applied the model in Equation (12)  lambs, respectively (Fig. 2b). As no strong statistical support is found for these differences, nor is 332 any biological justification known for such an effect of early pregnancy, LReS will be assumed to be ing not only that LReS is highest in ewes born at low population densities, but also that the effect 336 of lamb body mass on LReS is stronger in ewes that were also born in years of low population density. being a twin is required so that the quantities such as mean phenotype or selection gradients can be 343 average over this probability. For this purpose, a binomial mixed model was fitted according to where E[t ijk ] are individual probabilities of twinning, and other terms are defined analogously to 345 other models. Detailed information on the estimated quantities is provided in Appendix B, Table B7.

347
The analyses shown so far establish that (1) early pregnancy is size-dependent, (2) there is a fairly Equation umbers given in parentheses reference the statistical models used to parameterise each 362 component. As before, although corresponding to different quantities, we use the same Greek let-363 ters for equivalent coefficients in all models: α represents the intercepts, α t and α p the twinning and 364 pregnancy-specific contrasts to the intercept, respectively, β m and β m β m β m the slope and the linear (slope), 365 quadratic and cubic terms for mass, β d β d β d the linear and quadratic terms for population density and, 366 finally, β β β is a vector containing the effects associated with the remaining covariates. β β β also includes 367 interaction terms for body mass and population density, interaction terms between early pregnancy 368 status and population density, and between these two and body mass, in all models where the respec-369 tive predictors are included. u cj and u mk are random effects associated with cohort j and mother k, Mean phenotype,z z z, and the phenotypic variance-covariance matrix, P P P, on the expected scale were 382 obtained using Equations (5) and (2). We estimatez z z as , and P P P as shown in This selection gradient thus represents the ultimate effect of small changes in phenotype (possibly, 394 latent values for phenotypes that are modelled probabilistically, as is pregnancy) on fitness, scaled to 395 relative fitness, averaged over the range of conditions experienced by individuals in the population.

396
The algorithm for implementing the calculation described by equation 15 is detailed in appendix C.

400
Importantly, the conflict between viability and fecundity is density-dependent, and is almost nonexis-401 tent when population density is very low. In such circumstances, virtually all ewe lambs survive, and 402 therefore early pregnancy does not have a significant cost. As a consequence, η p is higher (i.e., less 403 strongly favouring low values) when population density is at its lowest and decreases with population 404 density (Fig. 4, Tab. 3). On the other hand, at high densities, when the probability of first-year 405 survival is lowest, selection favouring larger body size is stronger.

407
To investigate the effect of pregnancy on the total effect of mass on fitness, we calculated a direct 408 selection gradient, β m , as opposed to the extended directional selection gradient, η m . The direct 409 selection gradient describes that part of the effect of a focal trait (mass) on fitness that occurs inde-410 pendently of other focal traits (in this case, pregnancy). We calculated the direct selection gradient 411 using the twin model approach described in , which amounts to re-calculating the selection gradient 412 for mass, but using a modified path diagram (as Figure ?? very high (Fig. 4), fecundity governs selection on lamb body mass, as mothers will successfully rear 438 several offspring. As population size increases and first-year survival decreases (Fig. 4), viability se-439 lection becomes proportionally more important as it occurs before any opportunity for reproduction.

440
Under exceptional high population density, the probability of survival is very low. As most ewes 441 die, fecundity selection becomes slightly more important. That is not observed when the correlation 442 between lamb body mass and early pregnancy is dismissed, because the pregnancy-related viability 443 costs are not considered.
a version of Equation (2) where variance in breeding values is estimated, associated with both the 454 intercept (a int ) and the slope for body mass (a sl ). As for assessing selection, we (1)  lambs would be very high, regardless of pregnancy status ( Figure A1). Additionally, the survival of 541 offspring of ewes that became pregnant as lambs is extremely strongly density dependent (Table B5)      (c) size-dependent probability of early pregnancy. The line in teal corresponds to the prediction of a binomial regression with logit link function (see Appendix B for full specification and parameter estimates), whereas the grey bars correspond to 95% confidence intervals on the binomial probability obtained from the raw data; (d) lamb body mass (upper panel) and probability of early pregnancy (lower panel) across time. The black lines correspond to predictions of a linear and a binomial model (logit link function), respectively, correcting for within-year birth and measurement dates, maternal age and population density (linear and quadratic terms for the latter two covariates). Random effects for cohort and maternal identity were also considered. The grey areas illustrate the variation in lamb body mass and probability of early pregnancy attributed to varying population density in 1 standard deviation from the mean; (e) non-parametric survival curves of Soay sheep according to mother's age at conception. The p value corresponds to a Peto & Peto test for differences among Kaplan-Meier curves.  (a, b), and both lamb body mass and early pregnancy status (c, d). Dots correspond to observed data and curves to model predictions. Binomial regressions for first-year survival with logit link function were fit to the data. LReS was estimated in ewes that survived their first year of life and Poisson regressions with exponential link function were fit to the data. In all models, population density and maternal age were included as covariates (linear and quadratic terms were estimated for both), as well as birth and measurement dates (see Tab. B4 and B6 in App. B). Variance among cohorts and maternal identities were also estimated. included in all regressions, as well as the following linear interaction terms: between lamb body mass 724 and population density in models including lamb body mass as predictor, between early pregnancy status and population density in models including early pregnancy as a predictor, and between these 726 three traits in models including both lamb body mass and early pregnancy status as predictors.

727
Likewise, random effects to estimate among-mother and among-cohort variation were included in all 728 regressions. This fixed and random structure was applied to all statistical models except the one 729 modelling AReS (due to a lack of statistical power).         (14). In such circumstances, solving integrals using Monte Carlo simulations is particularly useful 735 as the error is independent from the integrand dimensionality, but mainly because other numerical 736 methods are likely to fail. Monte Carlo integration relies on the fact that the integral L g(l)f (l)dl (or 737 L g −1 (l)f (l)dl if applied to an inverse link function) corresponds to the expectation of g, E f [g(L)].

738
As a consequence, a random sample of latent values, L, generated using its density f can be used to 739 approximate the empirical average of g, where n is the number of draws from L.ḡ is also an approximation to L g(l)f (l)dl. This approx- where h is some small value. So, the extended selection gradient for mass would be obtained by 754 evaluating W for a certain value of latent mass l and at for l + h. This procedure was applied to