Inbreeding depression in an outbred stickleback population

Inbreeding depression refers to the reduced fitness of offspring produced by genetically‐related individuals and is expected to be rare in large, outbred populations. When it occurs, marked fitness loss is possible as large populations can carry a substantial load of recessive harmful mutations which are normally sheltered at the heterozygous state. Using experimental cross data and genome‐wide identity‐by‐descent (IBD) relationships from an outbred marine nine‐spined stickleback (Pungitius pungitius) population, we documented a significant decrease in offspring survival probability with increasing parental IBD sharing associated with an average inbreeding load (B) of 10.5. Interestingly, we found that this relationship was also underlined by a positive effect of paternal inbreeding coefficient on offspring survival, suggesting that certain combinations of parental inbreeding and genetic relatedness among mates may promote offspring survival. Our results demonstrate the potential for substantial inbreeding load in an outbred population and emphasize the need to consider fine‐scale genetic relatedness in future studies of inbreeding depression in the wild.


| INTRODUC TI ON
Mating between close relatives often reduces the fitness of offspring through morphological, physiological or reproductive defects, a phenomenon defined as inbreeding depression (Crnokrak & Roff, 1999).
Inbreeding depression may lead in the most extreme case to offspring unviability and eventually to population extinction (Kyriazis et al., 2021). Quantitative and population genomic approaches to study inbreeding depression have become increasingly important as the environmental pressures imposed by global changes are predicted to increase habitat fragmentation, reduce population sizes and consequently increase the probability of inbreeding in the wild (Frankham, 2010;Spooner et al., 2018;van Oosterhout, 2020).
The magnitude of inbreeding depression depends on the interplay between the effective population size (N e , the number of potential breeding individuals and the mean and variance in their reproductive output, Waples, 2022), the rate at which new mutations arise in the population and the efficiency of selection at eliminating harmful mutations (Bataillon & Kirkpatrick, 2000). Populations with large N e are large mutational targets but natural selection is more efficient against deleterious mutations in large than in small populations (Lanfear et al., 2014;Petit & Barbadilla, 2009). As a result, most lethal or nearly-lethal mutations segregating in populations with high N e are recessive, and can subsist in the population at low frequency in the heterozygous state (Charlesworth & Willis, 2009). Conversely in small populations, and particularly those confined to population isolates with no or limited migration, the probability of mating with genetically related individuals is drastically increased. Nonetheless, the cost of inbreeding in such populations should be reduced, as deleterious mutations would be | 3441 FRAIMOUT et al. rapidly lost through genetic drift or purged due to homozygotes' lethality (Bataillon & Kirkpatrick, 2000;Glémin, 2003;Hedrick & Garcia-Dorado, 2016).
Furthermore, the role of inbreeding depression in driving population viability has been primarily investigated in species characterized by long-term small populations sizes (e.g., endangered or domesticated populations; Bozzuto et al., 2019;Brekke et al., 2010;Hasselgren et al., 2021;Leroy, 2014;Norén et al., 2016;Saccheri et al., 1998). In large-sized populations, however, the cost of inbreeding has received less attention (but see Pérez-Pereira et al., 2021), possibly because they are considered to have reduced risk of extinction. Nevertheless, the cost of inbreeding in large populations is expected to be high, as it increases the chance of recessive mutations to be expressed as lethal homozygotes (Willi et al., 2013).
Demographic and genetic parameters specific to marine organisms make them opportune models to test theoretical predictions of inbreeding depression in the wild. Indeed, while several marine species are characterized by large N e , high dispersal and low genetic differentiation among populations (Bailleul et al., 2018;Hoey & Pinsky, 2018;Teixeira et al., 2012;Ward et al., 1994), it has been long recognized that population structure in the sea may be more prevalent than previously thought (Bierne et al., 2016;Cano, Mäkinen, et al., 2008;Cano, Shikano, et al., 2008;Corander et al., 2013;Eldon et al., 2016;Hauser & Carvalho, 2008) and that such structure can have major implications for evolutionary responses to selection Kemppainen et al., 2021).
Moreover, global environmental change and overexploitation of marine fish stocks (Pershing et al., 2015;Worm & Branch, 2012) may lead to rapid declines in population size and further exacerbate the effect of inbreeding depression, even in large-sized populations. In this context, population genomics studies can provide valuable estimates of population-specific inbreeding load and advance our understanding of inbreeding depression in wild outbred populations (Kardos et al., 2016).
The main objective of this study was to investigate whether cryptic relatedness estimated from genome-wide identity-bydescent (IBD) relationships among wild-collected parents influenced offspring survival in a large crossing experiment of the nine-spined stickleback (Pungitius pungitius). Specifically, we tested if relatedness among a set of wild parental fish used to produce F 1 -generation sibs through artificial mating led to inbreeding depression in offspring viability. Our results show that genetic similarity among parents can yield substantial reduction in offspring survival even when the level of relatedness is low.

| Study species and population
The nine-spined stickleback (Pungitius pungitius) is a small teleost fish with a wide circumpolar distribution and in Northern Europe, the species can be found throughout the Baltic Sea. Recent population genomic work has shed light on the demographic history of Baltic P. pungitius populations and shown that both genetic admixture (Feng et al., 2022) and relatively high levels of isolation by distance (Fang et al., 2020 characterize the distribution of genetic variation in these populations. Aside from these characteristics, Baltic

| Sampling and rearing
The fish used in this study came from a large crossing design used in previous studies of recombination rate variation in P. pungitius (Kivikoski et al., 2021(Kivikoski et al., , 2023. Prior to establishment of artificial crosses, wild individuals were housed in a 1 m 3 plastic aquarium with flow-through freshwater system and fed ad libitum with frozen chironomid larvae twice a day. The crossing scheme consisted of a maternal half-sib/full-sib design: 40 females where each mated to two different males and six females were mated to single males. Individuals used as parents for the artificial crosses were chosen at random (i.e., without distinction of morphological characteristics).
Parental fish we crossed artificially following a standard in vitro approach for sticklebacks (Arnott & Barber, 2000). Eggs were stripped from gravid females by gently squeezing their abdomen and placed in two different petri dishes. Males were over-anaesthetized using tricaine methanesulphonate (MS222) to dissect their testes, which were subsequently minced in the dish containing the eggs from their corresponding females. Fertilization was performed by mixing the eggs and sperm using a sterile plastic pipette and was ensured by looking for signs of post-fertilization membranes under a binocular microscope. In total, 46 females were mated to 86 unique males resulting in a total of 86 half-sib/full-sib families.
Following egg hatching and after larvae started feeding, each half-sib family was thinned to approximately 25 F 1 offspring per family and transferred to two large identical aquaria. Individuals in each family were divided into two groups at random to ensure that ca half of each family was represented in each aquarium. Offspring were mass-reared and their parental identity were subsequently recovered from genotype data as explained in (Kivikoski et al., 2021;and see below). Offspring survival was scored as the number of F 1 fish alive in given family at the end of the experiment. The experiment took place between May 2013 and April 2014, and survival was thus scored after a period of c. 48 weeks.

| Sequencing and genotyping
Parental fish were whole-genome sequenced (WGS) on an Illumina HiSeq platform (BGI, Hong-Kong) at 5-10x sequencing coverage and offspring were sequenced using the diversity arrays technology (DarTseq technology, Pty Ltd). Fastq files were mapped to the contig assembly of the P. pungitius reference genome (Varadharajan et al., 2019) using bwa-mem (Li et al., 2009) and then sorted bams for each file were obtained using SAMtools (Li et al., 2009). Genotype likelihoods were obtained using Lep-MAP3 (Rastas, 2017) which was used also for the linkage mapping and the pedigree construction. Offspring were assigned to parents using the Lep-MAP3's IBD module by calculating the Mendelian genotypic error rate for each combination of parents and offspring and by taking all combinations with error rate <10% (average error rate was 5.5%). Finally, offspring sex was obtained by comparing the sequencing coverage on X and Y specific regions using SAMtools depth. The final genomic data set for both parents and offspring consisted of 49,493 biallelic SNPs.

| Estimation of genomic inbreeding and relatedness
We estimated the genomic inbreeding coefficient (F G ) of all dams and sires from SNP data using the -ibc option implemented in the GCTA software (version 1.93.2, Yang et al., 2011). This metric provides a measure of inbreeding that is unbiased, as homozygote genotypes are scaled by their allele frequencies (Yang et al., 2011). Second, we estimated the pairwise relatedness between each pair of parents by computing the genomic relationship matrix (GRM) in GCTA using the "-make-grm" function. GCTA uses the allelic correlation coefficient (see Yang et al., 2011 p. 76 andGienapp et al., 2017) to construct the pairwise matrix of IBD coefficients among all parental individuals.
To ensure that our results were not influenced by the specific approaches used to estimate F G and IBD in GCTA, we replicated all analyses based on alternative estimates using the software PLINK (version 1.90, Purcell et al., 2007). We used the "-het option" in PLINK to obtain individual inbreeding coefficients and the "-genome" option to estimate the proportion of genome shared IBD between all pairs of parents. Inbreeding coefficients estimated with PLINK correspond to excess homozygosity and deviation from Hardy-Weinberg equilibrium. In PLINK, the GRM is constructed following the approach of VanRaden (2008) and based on the probability of IBD computed from identity-by-state values. We chose these two software as they provide two alternative approaches to the estimation of each parameter of interest and have previously been shown to improve detection of inbreeding depression using genomic data (Yengo et al., 2017) and particularly in populations with large N e (Alemu et al., 2021). For the sake of comparison, we follow the nomenclature of other studies (e.g., Bérénos et al., 2016;Kardos et al., 2016;Nietlisbach et al., 2019;Yengo et al., 2017) and hereafter refer to GCTA estimates of inbreeding coefficient as F UNI and PLINK estimates as F HOM . Similarly, estimates of relatedness using GCTA and PLINK are referred to IBD Yang and IBD VR , respectively.
Finally, to ensure the robustness of our results, all analyses were ran using estimates of F G and IBD from data pruned for minimum allele frequencies (MAF) and linkage disequilibrium (LD) using the R package SNPrelate (Zheng et al., 2012). Specifically, we applied

| Statistical analysis
We ran generalized linear mixed models (GLMMs) to investigate the effect of parental relatedness on offspring survival using the lme4 package (version 1.1-27.1; Bates et al., 2015) in R (version 4.1.1; R Core Team, 2021).
First, offspring survival was modelled using the "glmer" function in lme4 using a logit link function as: where y is the survival rate of a full-sib family coded as a two-column vector indicating the number of live individuals at the start of the experiment, and the number of survivors at the end of the experiment (thus providing the proportion of survivors per family). We included parental relatedness (IBD; the proportion of genome shared IBD between each pair of parents) as a fixed effect so that each full-sib family was assigned the IBD coefficient of its parents (i.e., the relatedness between its mother and father). Finally, we included dam identity as a random effect to account for the maternal half-sib structure of our data. Preliminary investigations indicated that the inclusion of the dam random effect led to a change in the sign of the effect of parental relatedness on offspring survival. Specifically, we found a significant negative relationship between offspring survival and parental relatedness when we did not include dam identity as random effect, but a positive relationship when including dam identity in the model (see Supporting Information Methods). This result suggested that the relationship between parental relatedness and survival at the dam level could not be generalized across the whole data. In other words, within each halfsib families, dams tended to have fitter offspring with the sire most related to them (see Supporting Information Methods). Conversely, between all families, offspring whose parents were on average more related to each other showed lower survival rate (Figure 1, Figure S1). Therefore, to account for these distinct within-and between-family effects, we used a technique of within-subject centring (van de Pol & Wright, 2009) by calculating the between-female (i.e., betweenfamilies) relatedness as the mean relatedness of a given female with each mate (mR) and added it as a fixed effect in our GLMMs. Similarly, we estimated the component of relatedness within half-sib families by normalizing individual values as wR = IBD-mR, for each pair of parents.
This component of within-relatedness can thus be fitted as the deviation in survival rate due to differences between the two sires within a maternal half-sib cross, from the average relatedness of an individual's parents (Fay et al., 2022;van de Pol & Wright, 2009). Finally, we were interested in evaluating the effect of parental inbreeding coefficient (i.e., the degree of inbreeding of each parent) on offspring survival.
To this end, we included the paternal and maternal inbreeding coefficients as fixed effects in our analyses.
Our final model had the following syntax: where F.sire and F.dam correspond to the paternal and maternal inbreeding coefficients, respectively, and y, mR, wR and (1|dam) are as described above.
We ran the model in model 2 for all combinations of SNP data filtering (relaxed, moderate and stringent) and genomic estimator (F UNI /F HOM ) while keeping consistency between software used to estimate F G and IBD (i.e., using F UNI / IBD Yang and F HOM /IBD VR together in the same models). Collinearity among the predictor variables was assessed using the variance inflation factor (VIF) analysis by using the "vif" function in the R package car (version 3.0-11, Fox & Weisberg, 2019). Finally, as the use of a logarithmic link function to GLMMs has been advocated to ensure unbiased estimates of inbreeding load (Nietlisbach et al., 2019), we also ran all models using the "log" link option in lme4 (see Appendix S1).

| Estimation of inbreeding load
Inbreeding depression can be quantified using the inbreeding load B estimated as the negative slope of the regression of offspring inbreeding coefficients on the trait values on the logarithmic scale (Charlesworth & Charlesworth, 1987;Morton et al., 1956;Nietlisbach et al., 2019). When applied to survival or viability data, B represents the predicted mortality in a hypothetical group of homozygote individuals carrying one lethal deleterious allele and is also referred to as the number of lethal equivalents (e.g., Armbruster & Reed, 2005;Nietlisbach et al., 2019;Ralls et al., 1988). For the estimation of B, inbreeding coefficients of offspring are usually obtained from pedigree relationships which allow estimating the expected IBD relationships within individuals (Wang, 2016). Here, as pedigree relationships would not capture any fine-scale variation of realized relatedness among parents (all parents are assumed to be unrelated), we estimated the magnitude of inbreeding load as the slope of the regression of parental relatedness (mR) on offspring survival from model 2 above. Finally, we calculated the mean and 95% confidence interval (CI) of B from 1000 bootstrap iterations for each combination of SNP data and genomic inbreeding estimator.
To ensure consistency of our results across statistical approaches and to make our results comparable to other studies, we also estimated B following the traditional approach using offspring inbreeding coefficients instead of parental relatedness as the main model's fixed effect (see Figure S2).

| Heritability estimates
We estimated the proportion of genetic variation explaining the variance in offspring survival by computing the heritability (h 2 ) of survival rate using an animal model (Lynch & Walsh, 1998)

of the form:
where y is the vector of offspring survival coded as a binary response variable (0 and 1 for dead and surviving individuals, respectively); μ is the population grand mean, Z A is the matrix of relatedness between individuals obtained from pedigree relationships; Z DAM is the random effect of dam identity and ε is the error term.
We further tested if offspring survival variation resulted from inbreeding depression by adding individual F G as a covariate in model 3. Indeed, inbreeding depression is expected to bias the estimation of the genetic variance components from animal models (Wolak & Keller, 2014) and we hypothesized that a downward bias in h 2 resulting from the added fixed effect of F would give further support to the role of inbreeding in survival rate variation in our data. Because genomic data to estimate offspring F G was only available for surviving individuals, we assigned to each offspring the mean F G value calculated from the survivors of each family (see also Supporting Information Methods). Animal models were fitted with the MCMCglmm R package (Hadfield, 2010) using the "threshold" family option appropriate for binomial regressions (de Villemereuil, 2018) and ran for 10,400,000 iterations with a burnin period of 400,000 and sampling every 1000th iteration. We ran standard model diagnostics for MCMCglmm models (de Villemereuil, 2018;Hadfield, 2010) by investigating the trace plots of all variance components and computing their effective sample size using the "effectiveSize" function implemented in the coda R package (Plummer et al., 2006) and checking for autocorrelations between iterations using the "autocorr" function from coda. We computed h 2 both on the latent and observed scales using the "QGparams" function of the QGglmm R package (de Villemereuil et al., 2016).

| Inbreeding coefficients and relatedness among the parents
We found differences in inbreeding coefficients between sexes using both metrics, with males having lower inbreeding coefficients than females (F UNI : t = −5.497, p < .001; F HOM : t = −5.541, p < .001, Figure S3). Relatedness coefficients estimated with PLINK and GCTA were moderately positively correlated (r = .573; p < .001). Overall, relatedness between pairs of mates was low (Table S1), as expected from a random sample of individuals mated at random. The pairwise IBD estimated with PLINK ranged from 0 to 0.237 and from −0.029 to 0.090 with GCTA. Negative values estimated with GCTA indicate that these pairs of individuals were less related than the pairs in the sample on average (Yang et al., 2011).

| Inbreeding depression
There was a negative effect of parental relatedness on offspring survival across all models (Figure 1, Table 1). For all combinations of SNP data and genomic estimators but one (F UNI & IBD Yang , Table 1,   Table S1), we found a significant negative effect of mR on offspring survival. In other words, triads of parents that were on average more genetically related had fewer survivors among their progenies.
Across all models, we found a significant positive effect of paternal inbreeding coefficient on offspring survival (F SIRE , p < .001, Table 1) indicating that more inbred males sired fitter offspring. Maternal inbreeding coefficient tended to be negatively associated with offspring survival, but this result was only significant when using the Overall, different types of data pruning did not change the conclusion of the statistical models but rather, differences appeared to come from the genomic estimator used (Tables 1 and 2, Tables S2   and S3). Estimates of inbreeding depression tended to be higher when based on GCTA estimators compared to PLINK estimators (Tables 1 and 2, Tables S2 and S3). Nevertheless, 10 out of 12 statistical tests performed showed a significant negative effect of parental relatedness on offspring survival (Table 1, Tables S3). Estimates of B varied depending on the genomic estimators rather than on data pruning, with estimates using GCTA showing higher values than PLINK, but also displayed relatively large confidence intervals around the mean ( Table 2, Table S2). The mean B across all models was 10.5. Results based on the clutch-averaged inbreeding coefficients (see Supporting Information Methods) were slightly lower than the above estimates with a mean B = 7.9 across models (Table S2). Finally, results obtained using the "log-link" function gave highly similar results for all GLMMs (Table S3).

| Heritability of offspring survival
We found moderate h 2 in offspring survival with h 2 = 0.17 (0.04;0.32) on the latent scale and h 2 = 0.13 (0.04; 0.25) on the observed scale.
Adding offspring inbreeding coefficients as covariate in the animal model reduced estimates of h 2 with h 2 = 0.05 (1.7e-09; 0.12) and h 2 = 0.06 (2.12e-09; 0.14) on the latent and observed scales, respectively (Figure 1). F I G U R E 1 Inbreeding depression in juvenile survival rate of nine-spined sticklebacks. The survival rate for each family (filled circles) is shown as a function of the relatedness coefficients between pairs of mates estimated with GCTA on the moderately pruned data (see Section 2: Materials and Methods). Blue line represents the fitted regression line from the ggplot "smooth" function and grey shadings are the 95% confidence intervals around the predicted values. Points colour shading corresponds to the inbreeding coefficient of the sire of each family.

| DISCUSS ION
The most salient finding of this study was the fairly severe inbreeding depression observed across all GLMMs, seemingly much more severe than in other taxa for juvenile survival. This result shows that large, outbred populations can harbour substantial mutational load expressed through inbreeding, and more importantly, that inbreeding depression can ensue even when levels of parental relatedness is low. In the following, we discuss the implications of these results for the study of inbreeding depression in the wild.
The magnitude of observed inbreeding depression exceeds that of a wide array of taxa previously reported in the literature based on estimated inbreeding load B (and see Section 2.6: Estimation of inbreeding load). Early estimates of lethal equivalents in several mammals and bird species are on average one order of magnitude lower for survival traits (e.g., Parus major, 0.84; Gazella spekei, 3.75; TA B L E 1 Results from the GLMMs on the effect of parental relatedness on offspring survival. For each model, the slope of regression for each coefficient of the fixed terms is given along with the standard error (SE) and the p-value.  Table 1 in Nietlisbach et al., 2019). Similarly, the estimates obtained in the current study appear to exceed those previously reported for marine fish for the similar traits (Thrower & Hard, 2009;Plough, 2016). Although our estimates present relatively wide confidence intervals and vary among the genomic estimators used, the order of magnitude and range of the inbreeding loads reported here is on par with severe loads previously reported from some natural populations (Brown & Kelly, 2020;Kennedy et al., 2014;Kruuk et al., 2002;Liberg et al., 2005). This has important conservation implications (Willi et al., 2022) as the results suggest that large-sized populations may be at a high risk of fitness loss if rapid decline in population size or habitat fragmentation would increase the probability of inbred matings.
The severity of inbreeding depression is known to depend on the environmental conditions under which it is expressed, with stressful conditions amplifying its severity (Bijlsma & Loeschcke, 2012). Since our study was carried out in a laboratory, we may have underestimated its real magnitude. Indeed, wild P. pungitius encounter parasites (Zander, 2007), In their study, Yuan et al. (2019) propose that such mating strategy (kin avoidance by inbred individuals) would alleviate the fitness cost of inbreeding, particularly in the sex investing more resources in reproduction. In P. pungitius, males provide parental care through nest-building, guarding and fanning of eggs and could thus benefit from this mating strategy. Unfortunately, due to multicollinearity among variables, our current data did not allow us to model the interaction between paternal inbreeding and parental relatedness and test this hypothesis further. However, we found that survival rate was heritable, indicating that genetic differences between families explain some amount of variation in offspring viability. Hence, certain combinations of parental inbreeding coefficients and genetic relatedness yielding higher survival could be favoured by selection and provide basis to the evolution of mate choice in this species (cf. Amos et al., 2001;Roff & Fairbairn, 2015). However, heritability of offspring viability was lowered substantially when accounting for individual inbreeding coefficients, suggesting that inbreeding affected survival to a larger extent than among-family genetic differences.
Experimental studies of mate choice controlling for the effect of parental inbreeding coefficient and parental relatedness should be particularly interesting to address these questions further.
While the results demonstrate that potential for inbreeding depression is present in the study population, the results should not be extrapolated to say that it would be realized in wild P. pungitius populations. This is because all matings in this study were artificial, and therefore eluded all aspects of mate choice, including kin recognition and active inbreeding avoidance. Such behaviours have been studied in the related three-spined stickleback (Gasterosteus aculeatus) in which gravid females prefer to lay eggs in nests built by unrelated males (Frommen & Bakker, 2006). However, such studies of mate choice and inbreeding avoidance often use very different levels of relatedness by contrasting completely unrelated to relatively closely related individuals (e.g., full-sib). In the present study, the range of relatedness among parental fish was narrow, and well below full-sib relationship (the maximum value of IBD VR = 0.237 corresponding to mating with aunt/uncle or among half-sibs). Hence, it is unclear whether kin recognition in sticklebacks would be fine-tuned to the levels of relatedness measured in our study and whether wild females could discern between males with such low levels of relatedness. Furthermore, high levels of local kinship have been reported in other marine fish (Aglieri et al., 2014;Eldon et al., 2016)  This result could be consistent with the possibility that dispersal in nine-spined sticklebacks is male biased as in the case of Baltic Sea three-spined sticklebacks (Cano, Mäkinen, et al., 2008;Cano, Shikano, et al., 2008). Alternatively, recombination rates have been shown to vary between sexes in P. pungitius, with females showing up to three times more autosomal crossovers during meiosis than males (Kivikoski et al., 2023). Whether sex-biased dispersal or sexspecific recombination rates could explain differences in inbreeding coefficient between males and females in our study is, however, unclear.
In conclusion, as predicted by theory, the results show that there is an opportunity for severe inbreeding depression to occur in large outbred marine populations. Future studies of inbreeding depression in the wild should benefit from the inclusion of fine-scale measurement of parental relatedness in models of inbreeding depression.

ACK N O WLE D G E M ENTS
We thank everyone who helped with catching and rearing the fish, and Takahito

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no competing interests.

AUTH O R CO NTR I B UTI O N S
A.F, J.M. L.L. and P.R. conceived the study; P.R. contributed to preprocessing and getting the genotype data; A.F. and L.L. analyzed the data; J.M. and A.F. led the writing of the manuscript; J.M. acquired funding to perform the study. All authors contributed critically to the drafts and gave final approval for publication.

DATA AVA I L A B I L I T Y S TAT E M E N T
Sequence data that support the findings of this study have been de-