Viability selection creates negative heterozygosity–fitness correlations in female Black Grouse Lyrurus tetrix

There is widespread interest in the relationship between individual genetic diversity and fitness–related traits (heterozygosity–fitness correlations; HFCs). Most studies have found weak continuous increases of fitness with increasing heterozygosity, while negative HFCs have rarely been reported. Negative HFCs are expected in cases of outbreeding depression, but outbreeding is rare in natural populations. Negative HFCs may also arise through viability selection acting on low heterozygosity individuals at an early stage producing a skew in the heterozygosity distribution. We tested this idea using survival and clutch parameters (egg mass, egg volume, chick mass, clutch size) in female Black Grouse Lyrurus tetrix and carried out simulations to determine how survival selection may impact the HFCs measured using clutch parameters. We show that survival is positively related to both individual heterozygosity and female body mass. There was a positive effect of body mass on all clutch parameters, but the selective mortality of females with both low heterozygosity and low body mass led to overrepresentation of high heterozygosity-low body mass females and hence a negative relationship between egg volume and heterozygosity. Using simulated data, we show that survival selection acting on both low body mass and low heterozygosity leads to a skew in the quality of breeding females, resulting in negative HFCs with egg volume. Our results indicate that survival selection can strongly influence the strength and direction of any HFCs that occur later in life and that only an integration of all aspects of individual reproductive investment and reproductive success can enable us to fully understand how heterozygosity can shape individual fitness. Selektion auf Überlebensfähigkeit führt bei weiblichen Birkhühnern Lyrurus tetrix zu negativen HFC-Werten Es besteht ein breites Interesse an den Beziehungen zwischen der individuellen genetischen Diversität und fitnessrelevanten Merkmalen (Heterozygosität-Fitness-Korrelationen; heterozygosity-fitness correlations; HFC). Die meisten Untersuchungen ergaben schwache kontinuierliche Fitnesszunahmen bei wachsender Heterozygosität, während von negativen HFC-Werten nur selten berichtet wird. Negative HFC-Werte sind in Fällen von Auszuchtdepression zu erwarten, wobei Auszucht in natürlichen Populationen nur selten vorkommt. Negative HFC-Werte können aber auch dadurch zustande kommen, dass die Selektion der Überlebensfähigkeit bei Individuen mit geringer Heterozygosität bereits in einem frühen Stadium greift und so eine Verschiebung in der Verteilung der Heterozygosität hervorruft, welche dann zu negativen HFC-Werten führt. Wir überprüften diese Theorie anhand von Überlebens- und Gelegeparametern (Eimasse, Eivolumen, Kükenmasse, Gelegegröße) von weiblichen Birkhühnern Lyrurus tetrix und führten Simulationen durch, um zu ermitteln, wie die Überlebensselektion die auf Grundlage der Gelegeparameter gewonnenen HFC-Werte beeinflusst. Wir zeigen, dass das Überleben in positiver Beziehung sowohl zur individuellen Heterozygosität als auch zur Körpermasse der Weibchen steht. Es gibt einen positiven Einfluss der Körpermasse auf alle Gelegeparameter; die selektive Sterblichkeit von Weibchen mit sowohl niedriger Heterozygosität als auch geringer Körpermasse führt zu einer Überrepräsentation von Weibchen mit hoher Heterozygosität und geringer Körpermasse und somit zu einer negativen Beziehung zwischen Eivolumen und Heterozygosität. Anhand von simulierten Daten zeigen wir, dass eine sowohl auf geringe Körpermasse als auch auf niedrige Heterozygosität wirkende Überlebensselektion zu einer Verschiebung in der Qualität der brütenden Weibchen führt, was negative HFCs bezüglich des Eivolumens zur Folge hat. Unsere Ergebnisse deuten darauf hin, dass die Überlebensselektion einen großen Einfluss auf Stärke und Richtung jeglicher zu einem späteren Zeitpunkt im Lebens auftretender HFCs ausüben kann und dass nur die Berücksichtigung sämtlicher Aspekte des reproduktiven Aufwandes und des Bruterfolges es uns ermöglicht zu verstehen, wie Heterozygosität die Fitness eines Individuums prägen kann.

Abstract There is widespread interest in the relationship between individual genetic diversity and fitness-related traits (heterozygosity-fitness correlations; HFCs). Most studies have found weak continuous increases of fitness with increasing heterozygosity, while negative HFCs have rarely been reported. Negative HFCs are expected in cases of outbreeding depression, but outbreeding is rare in natural populations. Negative HFCs may also arise through viability selection acting on low heterozygosity individuals at an early stage producing a skew in the heterozygosity distribution. We tested this idea using survival and clutch parameters (egg mass, egg volume, chick mass, clutch size) in female Black Grouse Lyrurus tetrix and carried out simulations to determine how survival selection may impact the HFCs measured using clutch parameters. We show that survival is positively related to both individual heterozygosity and female body mass. There was a positive effect of body mass on all clutch parameters, but the selective mortality of females with both low heterozygosity and low body mass led to overrepresentation of high heterozygosity-low body mass females and hence a negative relationship between egg volume and heterozygosity. Using simulated data, we show that survival selection acting on both low body mass and low heterozygosity leads to a skew in the quality of breeding females, resulting in negative HFCs with egg volume. Our results indicate that survival selection can strongly influence the strength and direction of any HFCs that occur later in life and that only an integration of all aspects of individual reproductive investment and reproductive success can enable us to fully understand how heterozygosity can shape individual fitness.

Introduction
There is widespread empirical evidence that individual multilocus heterozygosity can be correlated to fitness traits (David 1998;Hansson and Westerberg 2002;Coltman and Slate 2003;Chapman et al. 2009;Szulkin et al. 2010). The impact of heterozygosity on fitness (heterozygosity-fitness correlations; HFCs) is commonly explained by inbreeding effects across the whole genome (referred to as the general effect hypothesis) or by localized effects at single loci (referred to as the local/direct effect hypotheses) (David 1998;Hansson and Westerberg 2002). HFCs are typically positive linear relationships but they are generally weak and explain a low amount of the variance in fitness (Coltman and Slate 2003;Chapman et al. 2009). Negative or non-linear HFCs reflecting quadratic relationships between heterozygosity and fitness (Marshall and Spalton 2000;Neff 2004) are considerably less common in the literature (Chapman et al. 2009 but see e.g. Olano-Marin et al. 2011a). As inbreeding depression is generally expected to explain positive HFCs, outbreeding depression is generally expected to explain negative HFCs (LeBas 2002;Chapman et al. 2009;Jourdon-Pineau et al. 2012). However, outbreeding depression is rare in wild populations (Pusey and Wolf 1996) and positive and negative HFCs may co-occur through differing local or general effects at functional or neutral loci (Olano-Marin et al. 2011a, b).
Individual fitness measures overall genetic contribution to the next generation(s) and results from an individual's ability to extract and allocate limited resources to competing functions (growth, survival, reproduction). Fitness components such as age-specific survival probabilities and reproductive success are not independent (Williams 1966), meaning that individuals' performance at each age/stage may have consequences on the subsequent ones directly (e.g. early death) and indirectly (e.g. negative correlations of reproductive success in consecutive years). Therefore, viability selection, whereby an individual's survival probability depends on a given trait's expression, may influence all measured evolutionary processes occurring in later ages/stages (Hadfield 2008;Nakagawa and Freckleton 2010; see e.g. Mojica and Kelly 2010). In studies investigating HFCs, fitness components are often considered separately with, for example, heterozygosity and shortterm survival (Bean et al. 2004;Canal et al. 2014) or heterozygosity and clutch size (Ortego et al. 2007;Wetzel et al. 2012). When more complex approaches are undertaken, the traits considered (heterozygosity and morphological traits such as body mass) typically act simultaneously on a single component of fitness such as survival (Richardson et al. 2004). But the degree to which viability selection directly or indirectly associated with individual fitness influences the strength and direction of downstream HFCs is not yet well known (Fig. 1).
We explored this idea in female Black Grouse Lyrurus tetrix, a short-lived lekking species with reproductive skew to older dominant males (Alatalo et al. 1991;Lebigre et al. 2007;Kervinen et al. 2016). In this species, males are philopatric and female-biased dispersal is effective at reducing inbreeding (Lebigre et al. 2008(Lebigre et al. , 2010. However, inbreeding does occur and females have an increasing likelihood of mating with kin as they age (Soulsbury et al. 2011(Soulsbury et al. , 2012. Inbreeding and heterozygosity are known to be important in Black Grouse; offspring from inbred matings show reduced hatching mass (Soulsbury et al. 2011) and male reproductive success is positively related to individual heterozygosity (Höglund et al. 2002). We initially set out to explore the relationship of heterozygosity with two components of female fitness: survival, and clutch parameters (egg mass, egg volume, chick mass, clutch size). As there is little overlap between these two datasets, we generated a dummy dataset of individual heterozygosities to which we applied viability selection. We then used the resulting relationship to determine the degree to which viability selection influences the HFC measured subsequently using female clutch parameters.

Study sites and capture data
We used clutch data from Black Grouse captured at from one to seven winter feeding sites in Central Finland during 2001-2006. Adults were captured using walk-in traps, ringed, and aged as yearling (1 year old) or older (C2 years old) based on colouration of the outermost primary wing feathers (Helminen 1963). Females were then weighed and several morphometric measurements taken (including wing length, tarsus length). All females were fitted with metal rings and combinations of coloured plastic rings to aid long-distance identification. A subset of females was fitted with 16-to 20-g necklace radio transmitters (approximately 1.5-2.7% of body mass) with a mortality sensor, and lifetime of 104 weeks (Holohil RI-2D); recaptured females were fitted with new transmitters.
After the mating season (end of April-early May), the nest sites of radio-tagged females were localized and the hatching date of the brood was estimated by floating the eggs in warm water (Lebigre et al. 2007). Eggs were weighed to the nearest 0.1 g and measured (length and width) to the nearest millimetre. Egg volume was calculated using Hoyt's formula (Hoyt 1979 where k is a function of egg shape [in Black Grouse, k = 0.51 (Lindén 1983)]. Egg mass is more variable than egg volume because of the effect of embryo development on the measured egg mass. Eggs that have been incubated for a longer time are lighter. Females were followed for multiple years if they survived; however, female lifespan is short (average lifespan = 1.5 years; maximum age = 6) and often we did not recapture females in another year, thus lacked further body mass measurements for them. We only used clutch data where we had a body mass measurement for that year. Females may lay a second clutch when their first nesting attempt fails early in the season (typically because of predation). These replacement clutches are fundamentally different from first clutches in terms of size/mass/sex ratio (Soulsbury et al., unpublished data) and were excluded from the datasets.

Survival analysis and behavioural data
Transmitters did not influence female survival as recapture rates of ringed females with and without transmitters were similar (Pekkola et al. 2014). Monthly female survival was monitored by searching the study areas within an approximate 5-km radius by car and scanning the signals at elevated locations. During the incubation, a few females in more remote summer ranges or hidden by natural obstacles were located from an airplane within an approximate 10-km radius. We assessed total lifespan using females of exact known age (i.e. females caught as yearlings). As yearling females may continue to disperse up until their first breeding season (Marjakangas and Kiviniemi 2005), which may lead to biases in our measures of female survival, we assessed total lifespan (weeks) of yearling females located from calendar week 24 (n = 75 females 5391 bird-weeks monitored) onwards from their first breeding season.

Multilocus heterozygosity
Genomic DNA was extracted from the red blood cells using reagents from the BioSprint 15 DNA Blood Kit (reference 940017; Qiagen) and a Kingfisher magnetic particle processor. All individuals were genotyped at 11 autosomal microsatellite loci (BG6, BG15, BG16, BG18, BG19, TTD2, TTD3, TTT1, TUD6, TUT3, TUT4) and one-sex linked locus (BG10) (see Lebigre et al. 2007 for details). The software package CERVUS (Kalinowski et al. 2007) provided summary statistics for the data set including heterozygosity, number of alleles per locus and v 2 goodness-of-fit tests for deviation from Hardy-Weinberg equilibrium (HWE). Inbreeding occurs in this population [*13% of all clutches were produced by close relatives (Lebigre et al. 2010)]. We used standardised multilocus heterozygosity (sMLH) as a measure of female MLH; sMLH was calculated from all adult genotypes using the Rhh package (Alho et al. 2010) run in R version 3.2.1 (R Core Team 2014). We tested whether our markers captured genome-wide heterozygosity by calculating heterozygosity-heterozygosity correlations and tested for identity disequilibrium with g 2 using the inbreedR package (Stoffel et al. 2016). However, our statistical power to capture genome-wide heterozygosity is likely to be biased by our low numbers of loci (e.g. Chapman and Sheldon 2011).

Data analyses
We first tested the relationship between female total lifespan (weeks) in relation to sMLH and female mass at first capture using linear models. This meant we had estimates for both sMLH and body mass that we could subsequently use in our dummy dataset. Using data for these females, we also correlated female mass and sMLH.
Using a separate dataset of females (n = 130-135 females/142-147 clutches depending on clutch parameter), we then tested the relationship between sMLH and female mass on five clutch parameters (mean egg mass, mean egg volume, mean chick mass, clutch size, clutch sex ratio). Clutch parameters are often collinear, so we analysed them separately as they may show different responses to female reproductive investment. We predicted that egg size (egg mass, egg volume) and clutch size would correlate with sMLH because previous work showed chick mass was related to pairwise relatedness of parents (Soulsbury et al. 2011). In all models, sMLH and (log-transformed) female body mass were set as fixed factors and female ring number was included as a random factor. We chose not to use female age (i.e. yearling or C2 years old) as a cofactor because (1) female age and body mass are interdependent (yearling females are 6-7% lighter than older females; mean individual increase in mass ± SD = 57.0 ± 1.3 g) and because body mass is more likely to directly influence female clutch parameters than age per se. Egg mass, egg volume and chick mass were analysed using the lmer function from the lme4 package (Bates et al. 2011) and lmerTest package (Kuznetsova et al. 2013). To account for the underdispersion of clutch size data (Devenish-Nelson et al. 2013), we used a generalized linear mixed-effects model using the glmmPQL function with a quasi-Poisson family from the mass package (Venables and Ripley 2002). These models used clutch size as the dependent variable, sMLH and female body mass as fixed factors, and female ring number as a random variable. Lastly, we tested whether there could be a trade-off between egg mass, volume, chick mass and clutch size using Pearson's correlations. All analyses were run in R version 3.2.1 (R Core Team 2014).

Dummy dataset generation and analysis
Our survival dataset and our clutch datasets are essentially separate, with very little overlap in data. To examine whether the combined effect of sMLH and mass were impacting survival, we created a dummy dataset. We used the standardised heterozygosity values of 479 females. Each individual sMLH was assigned a body mass of 740-1110 g at 10-g increments which represented the lightest and heaviest female in our dataset. In total, there were 17,723 females in our dummy dataset.
Predicted lifespan was calculated using coefficients derived from the linear model of female lifespan and sMLH and body mass (see Results section). For each individual we then calculated predicted egg mass, egg volume, chick mass and clutch size based on linear mixedeffects models of clutch parameters and log female body mass only. To avoid large sample sizes impacting the results, we random resampled 1000 individuals from this initial dummy dataset. We then tested the relationship between clutch parameters, and sMLH on (1) all data, and (2) individuals predicted to survive C24 weeks. We generated 1000 randomly sampled datasets, for which we extracted the coefficients from the model. We then compared estimates with and without viability selection using two-sample t-tests and tested differences in frequency distributions using a two-sample Kolmogorov-Smirnov test.

Genetic profile of population
A total of 1657 full-grown Black Grouse were genotyped leading to a total number of alleles for each locus which ranged from five to 23 (mean = 12 alleles; Table S1). Deviations from the Hardy-Weinberg expectations were found at one locus (BG15 ; Table S1), and may be partly explained by the weak population subdivision observed in this study system [the Wahlund effect (Lebigre et al. 2008)] and the pooling of samples collected in different years within each lek (the presence of transgenerational close relatives might strengthen deviations from the HWE). There was a non-significant positive heterozygosityheterozygosity correlation (r = 0.01; 95% CI = -0.03/ 0.05). Identity disequilibrium estimated from these loci did not differ significantly from zero in the complete dataset (g 2 = 0.006; bootstrap confidence interval = -0.01-0.02, P = 0.24).
Moreover, female body mass was positively related to all clutch parameters (Table 1; Fig. 3a-c) while sMLH was only negatively related to egg volume (Table 1; Fig. 4).
There was no relationship between clutch size and mean egg mass (r p = 0.10, P = 0.205), mean egg volume (r p = -0.02, P = 0.836), and mean chick mass (r p = -0.02, P = 0.759) suggesting that the females were not trading-off investment in eggs for a greater number of eggs. Similar to the lifespan data, there was a non-significant negative relationship between female mass and sMLH (linear mixed model, estimate ± SE = -0.01 ± 0.03, t = -0.46, P = 0.649).

Modelling linked HFCs
Without applying the effect of survival selection, coefficients of the relationship between sMLH and clutch parameters were close to 0 but consistently positive (Table S2; Fig. 5). After applying survival selection, there was a significant reduction in all measured coefficients (all models, t 1998 = 14.41, P \ 0.001) which became consistently negative (Fig. 5). There was significant difference in the frequency distributions (all models, D = 0.294, P \ 0.001).

Heterozygosity and fitness
Neither identity disequilibrium (g 2 ) nor heterozygosityheterozygosity correlations were significantly different from zero. This is not surprising, as most reported g 2 are non-significant [mean g 2 = 0.007 with only 26 of 129 estimates of g 2 significantly different from zero (Miller and Coltman 2014)]. Without variance in inbreeding (i.e. when g 2 = 0) HFCs cannot arise (Szulkin et al. 2010). However, inbreeding (Lebigre et al. 2010) and inbreeding depression occur in our study population [e.g. chick mass (Soulsbury et al. 2011)]. Non-significant values of g 2 do not contradict the observation of HFC, and low levels of inbreeding are typically easier to detect through its effect on the phenotype or fitness, rather than its effect on heterozygosity. At the same time, our genotyping data are from adult birds in the population. It is likely that early viability selection (e.g. at the chick stage or prior to sampling) lowers g 2 and weakens heterozygosity-heterozygosity correlations.
Viability selection results from the differences in individual phenotypic and genetic quality leading to the selective mortality of lower quality individuals. The nonrandomness of viability selection leads to a skew in the distribution of traits of interest and may strongly influence their observed relationship with fitness (Hadfield 2008;Nakagawa and Freckleton 2010). Among phenotypic traits, body size and mass are often linked to juvenile and adult survival as, across multiple taxa, heavier and larger individuals have higher survival rates [e.g. mammals (Milner et al. 1999); birds (Haramis et al. 1986;Monaghan and Metcalfe 1986)]. This effect is, however, often stronger in juvenile individuals than adults (Festa-Bianchet et al. 1997), though this probably depends on the life history of the species involved. Consistent with these previous findings, female Black Grouse with longer lifespans were heavier. In this species, females have relatively short lifespans, and the main source of mortality is Goshawk Accipiter gentilis predation (Angelstam 1984;Pekkola et al. 2014). Most predation occurs during the breeding season (Pekkola et al. 2014), and lighter females may have an increased need to feed during the active periods of predators leading to increased predation risk (e.g. Selås 2003;Wirsing et al. 2002). Mortality of lighter females may also occur after the breeding season if the energetic investment of rearing offspring reduces their ability to escape predation. In addition to body mass, heterozygosity also positively influenced the survival of female Black Grouse. This is similar to the findings of many other studies where genome-wide heterozygosity and heterozygosity at specific loci have been associated with increased survival in some vertebrate species [mammals (Banks et al. 2010;Forcada and Hoffman 2014); fish (Evans et al. 2010); birds (Worley et al. 2010;Cézilly et al. 2016)]. Similar to body size, the relationship between survival and heterozygosity has been mainly found in juveniles (e.g. Bean et al. 2004;Acevedo-Whitehouse et al. 2006;Cohas et al. 2009), but there is growing evidence that these effects continue throughout adulthood (Velando et al. 2015;Cézilly et al. 2016). We do not have data directly linking individual chick mortality and heterozygosity, but it is likely to be significant because of the strong relationship between inbreeding depression and chick mass (Soulsbury et al. 2011), which in itself is a major determinant of chick survival (Ludwig et al. 2010). Our data show that the effects of individual heterozygosity persist throughout adulthood and that, combined, both heterozygosity and body size determine survival in female Black Grouse.

Heterozygosity, body size and reproduction
Consistent with other studies, we found that clutch parameters were positively associated with body mass (Christians 2002). We also found a negative relationship between heterozygosity and egg volume, with weak negative trends for egg mass. While several studies have reported a positive relationship between clutch size and heterozygosity (Foerster et al. 2003;Ortego et al. 2007;Wetzel et al. 2012), few studies have examined the relationship between heterozygosity and other clutch parameters. Indeed, the only other clutch parameter examined was egg size, and some studies found a positive effect of heterozygosity (Forstmeier et al. 2012;Wetzel et al. 2012) or no effect (Ortego et al. 2007) probably because individuals trade-off egg size with clutch size (Smith et al. 1989;Nager et al. 2000). Even so, a negative relationship is unexpected. The negative HFC was weaker for egg mass, but this is not surprising since it is impacted by the timing of measurement and seems to be less responsive to female mass (see Christians 2002). The lack of any relationship between female heterozygosity and chick mass probably reflects the overriding effect of inbreeding depression acting on this trait (Soulsbury et al. 2011).

Viability selection impacts HFCs
Our empirical data show that viability selection is strongest on females with both low body mass and low sMLH. We would expect this to create a skew in the quality of breeding female in which low body mass-high sMLH females would be overrepresented leading to a negative relationship between female mass and sMLH, which could lead to negative HFC.
Our simulated data showed that applying viability selection lowered the estimated relationship between sMLH and clutch parameters, in all cases turning them Fig. 3 Relationship between log female body mass and a egg mass (g), b egg volume (cm 3 ) and c chick mass (g) Fig. 4 Relationship between egg volume (cm 3 ) and standardised multilocus heterozygosity (sMLH) Fig. 5 Mean ± SD regression coefficients for the relationship between clutch parameters and sMLH. Coefficients before (open square) and after (filled squares) survival selection are show. SD values for clutch size are too small to be seen on the graph negative. Therefore, viability selection acting at any age class is likely to impact subsequent HFCs leading to negative HFCs, as in this study, but it can also equally lead to a weakening of the already weakly positive HFCs [effect sizes = 0.05-0.10 (Chapman et al. 2009)]. A key future direction for study would be to combine multiple agents of selection and assess how they impact HFC direction and strength.

Conclusion
Our results indicate that positive and negative HFC can cooccur for genome-wide heterozygosity measures. This is supported by a strong relationship between heterozygosity, body mass and survival. This creates variance in female quality, meaning that breeding females with low heterozygosity tend to have greater body mass. Hence, these individuals have more resources available for reproduction, which creates a negative HFC. Fitness is affected by multiple factors, and our results indicate that the effect of heterozygosity on different fitness components cannot be considered solely in isolation. Viability selection is likely to have important downstream effects on reproduction and the subsequent relationships with HFC eventually leading to negative HFCs.