Genomic and fitness consequences of inbreeding in an endangered carnivore

Reduced fitness through genetic drift and inbreeding is a major threat to small and isolated populations. Although previous studies have generally used genetically verified pedigrees to document effects of inbreeding and gene flow, these often fail to capture the whole inbreeding history of the species. By assembling a draft arctic fox (Vulpes lagopus) genome and resequencing complete genomes of 23 additional foxes born before and after a well‐documented immigration event in Scandinavia, we here look into the genomic consequences of inbreeding and genetic rescue. We found a difference in genome‐wide diversity, with 18% higher heterozygosity and 81% lower FROH in immigrant F1 compared to native individuals. However, more distant descendants of immigrants (F2, F3) did not show the same pattern. We also found that foxes with lower inbreeding had higher probability to survive their first year of life. Our results demonstrate the important link between genetic variation and fitness as well as the transient nature of genetic rescue. Moreover, our results have implications in conservation biology as they demonstrate that inbreeding depression can effectively be detected in the wild by a genomic approach.


| INTRODUC TI ON
Anthropogenic activities are causing many wild populations to rapidly decline in numbers (Dirzo et al., 2014), resulting in depleted genetic diversity and exposing them to the negative effects of genetic drift and inbreeding (Hedrick & Kalinowsky, 2000;Keller & Waller, 2002). Inbreeding can lower individual fitness and hamper population growth (i.e., inbreeding depression) which is generally attributed to increased expression of recessive deleterious alleles (Charlesworth & Willis, 2009). However, inbreeding depression can be circumvented through gene flow from unrelated individuals by masking partially recessive deleterious alleles, resulting in increased genetic variation, individual fitness and population growth (i.e., genetic rescue; Ingvarsson, 2001;Tallmon et al., 2004).
Although inbreeding depression and genetic rescue have been of scientific interest for centuries (Hasselgren & Norén, 2019), studying these phenomena in the wild is challenging, as obtaining genetically verified pedigrees linked to life history traits requires detailed individual monitoring and DNA-collection over long time series (Pemberton, 2004(Pemberton, , 2008. During the last two decades, inbreeding depression has been studied through pedigree analyses in a range of wild animal taxa, although studies on mammals and | 2791 HASSELGREN Et AL. birds predominate (e.g., Åkesson et al., 2016;Hogg et al., 2006;Jamieson et al., 2006;Kruuk et al., 2002;Liberg et al., 2005;Marr et al., 2006;Nielsen et al., 2012). Even though genetically verified pedigrees are the foundation for studying inbreeding depression, they have several shortcomings. Pedigrees generally fail to capture historical inbreeding and relatedness among founders or immigrants (Kardos et al., 2016), resulting in downwardly biased inbreeding coefficients (Hedrick & Garcia Dorado, 2016;Kardos et al., 2018;Robinson et al., 2019;Wang, 2016). Moreover, individuals with the same level of relatedness obtain identical pedigree based inbreeding coefficients (F PED ) but might vary in reality due to Mendelian segregation, linkage and limited recombination (Taylor et al., 2015).
Recent advances in genomics enable accurate quantification of inbreeding without requirement of detailed pedigree information (Kardos et al., 2016). By sequencing complete or near complete genomes with high enough coverage to call heterozygous genotypes, inbreeding can be measured continuously across the genome by identifying regions in runs of homozygosity (ROH, i.e., long chromosomal stretches that are identical by descent [IBD] and a signature of inbreeding; Allendorf, 2017). The distribution of ROH length is informative of the time to coalescence of IBD haplotypes. Very long ROH are found when parents share a recent common ancestor. Short ROH indicate coalescence in more distant ancestors, and mutations in these regions have been subject to selection over longer time frames (Browning & Browning, 2015;Li et al., 2006). For instance, extremely long ROH stretching over entire chromosomes were found in severely inbred populations of grey wolves (Canis lupus) in Iberia, Scandinavia and North America (Gómez-Sánchez et al., 2018;Kardos et al., 2018;Robinson et al., 2019). Similar patterns were also found in Island foxes (Urocyon littoralis; Robinson et al., 2016) and North American pumas (Puma concolor; Saremi et al., 2019). Empirical support for a link between genomic inbreeding and fitness was for a long time rare for wild populations. However, developments within the field have recently enabled such studies (Bérénos et al., 2016;Chen et al., 2016;Harrisson et al., 2019;Hoffman et al., 2014;Huisman et al., 2016;Niskanen et al., 2020).
In an isolated subpopulation of the endangered Scandinavian arctic fox (Vulpes lagopus), both inbreeding depression and genetic rescue have recently been documented through pedigree analyses linked to life history traits (Hasselgren et al., 2018;Norén et al., 2016). Isolation at an extremely small population size led to a rapid increase in inbreeding, with lower survival and reproduction in inbred individuals (Norén et al., 2016). Recently, an event of natural immigration by three outbred males, of which two were brothers, resulted in genetic rescue through elevated fitness in F1 offspring between immigrants and native individuals. However, no increase in fitness was detected in F2 or F3 offspring with immigrant ancestry (Hasselgren et al., 2018;Lotsander et al., 2021). Moreover, no significant difference in microsatellite heterozygosity before versus after immigration was found (Hasselgren et al., 2018).
Here, we investigate the genomic consequences of inbreeding and outbreeding over multiple generations by generating a draft arctic fox genome and resequencing high-coverage genomes of 23 individuals collected before and after the immigration event.
Specifically, we quantify realized inbreeding as genome-wide heterozygosity and genomic inbreeding coefficients (F ROH ) and compare these measures between native foxes, immigrant F1 (i.e., crosses between immigrants and native individuals) and immigrant F2 + F3 (i.e., immigrant F1 or F2 that have backcrossed with native individuals). Furthermore, we investigate the relationship between inbreeding coefficients (F ROH ), genome-wide heterozygosity and juvenile survival.

| Study population
During the end of the 19 th century, the Scandinavian arctic fox population experienced a major demographic and genetic bottleneck due to intensive hunting for fur (Lönnberg, 1927). In the 1920's, the population was protected by law, but has nevertheless recovered slowly, which can primarily be explained by increased competition with red fox (Angerbjörn et al., 2013;Elmhagen et al., 2017) and reduction in rodent abundance (Angerbjörn et al., 2013) in the mountain tundra. Loss of genetic variation due to inbreeding and genetic drift is also probably an important factor hampering population recovery (Hasselgren et al., 2018;Norén et al., 2016).
The southernmost population in Sweden, residing in Helags mountains (3,400 km 2 ; Figure 1) was functionally extinct from the 1980s, but was naturally refounded by four individuals of unknown origin in the year 2001 (Norén et al., 2016). Three additional individuals of unknown origin established during the following years and the population was thus founded by seven individuals (four males and three females; Norén et al., 2016). Until 2010, the population was isolated at extremely small population size (N = 0-26 reproducing individuals) and inbreeding depression was documented through lower juvenile survival and reproduction in inbred individuals (Norén et al., 2016). However, in 2010 and 2011, three male arctic foxes (two brothers and one unrelated individual), released from the Norwegian captive breeding programme (Landa et al., 2017) immigrated into the population (Figure 1). The parents of the unrelated immigrant originated from Reisa Nord and Saltfjellet in northern Norway. The parents of the immigrant brothers were born within the captive breeding programme but had ancestry from five different subpopulations through grandparents and great grandparents (Figure 1). The immigrants produced at least 63 offspring in total (Lotsander et al., 2021), and 5 years after immigration, genetic rescue was evident through 1.9 times higher survival and 1.3 times higher breeding success in immigrant F1 offspring compared to native ones (Hasselgren et al., 2018).
During the same time period, the population increased from 26 to 58 breeding individuals. In 2015, 89% of the litters produced were descendants of immigrants (Hasselgren et al., 2018).

| Juvenile survival
Known arctic fox den sites were visited every summer to document survival and reproduction. This is an effective method to detect survivors as arctic foxes utilize the same dens every year rather than creating new ones (Dalerum et al., 2002). Juveniles were ear-tagged with unique colour combinations and first year survival was based on visual observations and genetic assignment (Hasselgren et al., 2018;Lotsander et al., 2021;Norén et al., 2016). Individuals that were visually or genetically verified at least 1 year after they were born were classified as first year survivors. Those that were not observed in the population again and had not appeared in the subpopulations (which are also closely monitored) were assumed to be nonsurvivors.

| De novo assembly
We generated a de novo reference genome for arctic fox (Vlagopus_ NRM_v1.fa; Genbank accession number: JAFLEK000000000) by sequencing a combination of Hi-C, mate-pair and short insert libraries from high molecular weight DNA extracted from one male and one female from the study population, born in 2015 and 2013, respectively. The female was a native individual with F PED = 0.10 and the male was an immigrant F2 with F PED = 0.05. These individuals were only used for the de novo assembly and were not included in the study specimens. DNA was extracted using the Thermo Scientific KingFisher Cell & Tissue DNA kit. A short-read assembly was generated from the female material and Hi-C scaffolding was accomplished using the male material. The following short-read libraries were prepared: one Illumina TruSeq PCR-free sheared to 180 bp insert size and three Nextera mate-pair libraries with insert sizes of 3 Kb, 5-8 Kb and 20 Kb. The short-read data was assembled using Allpaths-LG (Gnerre et al., 2011) to a draft genome assembly with an N50 scaffold value of 22 Mb. Frozen muscle tissue from the male individual was prepared using the Dovetail Hi-C kit (Dovetail F I G U R E 1 Arctic fox distribution in Scandinavia with study population, captive breeding station, release site and origin of three immigrants marked out. Two of the immigrants were brothers and their origin is represented by blue stars. White stars represent origin of the unrelated immigrant. The parents of the unrelated immigrant were wild caught from two different subpopulations. The parents of the immigrant brothers where born within the captive breeding programme and originate from five subpopulations through grandparents and great grandparents [Colour figure can be viewed at wileyonlinelibrary.com] Genomics) to generate a Hi-C proximity ligation library. The final genome assembly was done using the HiRise pipeline (Putnam et al., 2016), using the short-read assembly as input. To identify sexlinked chromosomes, the genome was blasted against the red fox X-chromosome (Kukekova et al., 2018) using BLAST+2.5.0 (Camacho et al., 2009). Repeats and transposable elements were de novo predicted from the genome assembly using RepeatModeler v.1.0.8 (http://www.repea tmask er.org/Repea tMode ler.html). The assembly was masked in RepeatMasker v.4.0.7 using the repeat library from RepeatModeler as input (http://www.repea tmask er.org). Finally, we identified CpG sites (all sites where a C nucleotide is followed by a G nucleotide in the reference genome) using a custom script (https://github.com/tvdva lk/find_CpG). The final size of the de novo assembly was 2.36 Gb and comprised 3,731 scaffolds with an N50 of 132 Mb.
We called variants using bcftools mpileup v1.6 (Li, 2011;Li et al., 2009) and bcftools call v1.6 using a minimum depth of coverage of 8× (i.e., one third of the average coverage) and maximum of 50× (i.e., twice of the average coverage). We then filtered SNPs by base quality QV ≥ 30 and those within 5 bp of indels. For all downstream analyses, we then excluded the X chromosome, mitogenome and hard masked repeat regions using BEDtools v2.27.1 (Quinlan & Hall, 2010). After merging all individual vcf files we used PLINK v1.9 (Purcell et al., 2007) to filter out SNPs out of Hardy-Weinberg equilibrium (HWE) and obtained 6,792,779 high quality SNP calls. For all downstream analyses, we only kept variants called in all individuals making for a total of 4,579,222 high quality SNPs.

| Heterozygosity and genomic inbreeding coefficients
We first estimated the number of heterozygous sites per 1 kb with the population mutation rate (θ), which approximates the expected heterozygosity per site under the infinite sites model and uses bam files as input using mlRho (Haubold et al., 2010). We filtered out bases with quality below 30, reads with mapping quality below 30 and positions with root-mean-squared mapping quality below 30.
Because high or low coverage in some regions, resulting from structural variation, can create erroneous mapping to the reference genome and false heterozygous sites, we also filtered out sites with depth lower than 8× (i.e., one third of the average coverage) and higher than 50× (i.e., twice of the average coverage) across all our specimens.
Secondly, to estimate levels of inbreeding based on runs of homozygosity (F ROH ), we used PLINK v1.9 (Purcell et al., 2007). The high coverage of genomes (19×-40× with 95% of genomes covered at a depth larger than 10×) in combination with high contiguity of the reference genome (N50 = 132 Mb) makes our data well suited to obtain realistic ROH estimates. We used a sliding window size of 200 SNPs (0.2 kb; homozyg-window-snp 200). A window was then defined as homozygous if there were not more than three heterozygous sites per window (homozyg-window-het 3). If at least 5% of all windows that included a given SNP were defined as homozygous, the SNP was defined as being in a homozygous segment of a chromosome (homozyg-window-threshold 0.05). This threshold was chosen to ensure that the edges of a ROH are properly delimited. A homozygous segment was defined as a ROH if all of the following conditions were met: the segment included ≥100 SNPs (homozygsnp 100) and covered ≥100 kb (homozyg-kb 100). Furthermore, the minimum SNP density was one SNP per 50 kb (homozyg-density 50) and the maximum distance between two neighbouring SNPs was ≤1,000 kb (homozyg-gap 1,000). In order to assess the sensitivity of our estimates to the parameters chosen, we also identified ROH using alternative parameters, varying the window size from 200 to 100 SNPs and maximum number of SNPs within a window from 3 to 1 ( Figure S1). We found that length of ROH did not change much when varying the window size, but decreasing the maximum number of heterozygous SNPs within a window to 1 tended to break up ROH.
Genomic inbreeding coefficients (F ROH ) were calculated as the fraction of the total genome consisting of runs of homozygosity. We explored number of generations back to the common ancestor for each ROH by using the following equation: g = 100/(2rL), where g is the number of generations, L represents length of ROH in Mb and r is the recombination rate. We used the recombination rate for silver fox (Vulpes vulpes; 0.6 cM per Mb; Kukekova et al., 2007). F ROH was calculated both for all ROH larger than 100 kb (expected to arise from ancestors approximately 850 generations back in time), for ROH larger than 2 Mb (reflecting coalescence ~45 generations back in time, corresponding to the time period the arctic fox became endangered in Scandinavia) and for ROH larger than 8 Mb (reflecting coalescence ~10 generations back in time, corresponding to the time period when the study population was re-established). Worth noting is that the inbreeding history retrieved from the equation should not be considered an exact measure as the recombination rate for arctic fox is not known. The variance is also expected to be large since we have measured ROH in Mb rather than genetic map length (cM).
Additionally, we used an average generation time of 2 years (Dalén et al., 2005), but in reality the generation time of Scandinavian arctic foxes varies with the phase of the rodent cycle. Foxes born during the increase phase often reproduce as 1-year-olds but those born during the peak phase often wait until they are 3 years old (Hasselgren et al., 2018;Norén et al., 2016). To rule out that depth of coverage influenced estimates of F ROH and heterozygosity, we plotted these parameters against each other. The R 2 -values varied between 0.01 and 0.02 ( Figure S2). Finally, we calculated inbreeding coefficients from a pedigree (F PED ) based on 11 autosomal microsatellite loci that were retrieved from previously published data (Hasselgren et al., 2018;Norén et al., 2016).

| Statistical analyses
We fitted linear regressions to investigate whether heterozygosity (measured as number of heterozygous SNPs per kb), F ROH and length of ROH differed between the three groups: native individuals, immigrant F1 offspring and immigrant F2 + F3 offspring. We fitted generalized linear models with binomial distribution and logit links to test if heterozygosity and F ROH influenced the probability of juvenile survival. We performed the same test to investigate whether pedigree inbreeding (F PED ) gave a similar estimate of inbreeding depression as F ROH . To test how well F PED correlated with genome-wide heterozygosity and F ROH we performed correlation tests and to find out whether F PED underestimated inbreeding we performed paired t tests. Tests on F ROH were accomplished first considering ROH larger than 100 kb (arising from ancestors <850 generations back) and thereafter considering only ROH larger than 2 Mb (<45 generations back) and 8 Mb (<10 generations back), respectively. All tests were performed in R version 3.6.1.

| De novo assembly, depth of coverage and runs of homozygosity
We generated a de novo assembly sequencing a combination of Hi-C, mate-pair and short insert libraries from two individuals from the study population. The final size of the de novo assembly was of 2.36 Gb and comprised 3,731 scaffolds with an N50 of 132 Mb. 99% of the genome was comprised within 24 scaffolds >30 Mb which probably correspond to chromosomal size levels. We identified 12 sex-linked scaffolds when we blasted the genome against the Xchromosome of red fox (Vulpes vulpes).
We performed whole-genome resequencing of 23 arctic foxes and obtained a depth of coverage of 19×-40× with an average of 25× (Table S1). We found that many individuals had very long ROH, stretching over nearly entire scaffolds (Figure 2). The longest ROH was 61.2 Mb and was found in a native female born in 2010, which corresponds to the year with the highest average inbreeding previously observed in the population (Hasselgren et al., 2018 Figure S3). There was no significant difference between length of ROH between native foxes and immigrant F2 + F3 F I G U R E 2 Runs of homozygosity (ROH) (coloured blocks) reflecting inbreeding due to shared ancestors less than 10 generations back (>8 Mb), mapped along the 23 largest scaffolds in native individuals and immigrant descendants of arctic foxes. *Immigrant F3; **founder [Colour figure can be viewed at wileyonlinelibrary.com] (95% CI [-0.09-1.29], df = 20, p = .09) or between immigrant F1 and immigrant F2 + F3 (95% CI [-0.26-1.37], df = 20, p = .19).

| Inbreeding coefficients (F ROH ) and individual ancestry
We estimated individual genomic inbreeding coefficients (F ROH ) as the proportion of the autosomal genome comprised of ROH. Overall, there was a wide range in F ROH coefficients; from 0.11 to 0.52 when considering inbreeding due to both recently and historically shared ancestors (<850 generations, ROH >100 kb) and from 0.004 to 0.35 considering inbreeding from only the most recent ancestors (<10 generations, ROH >8 Mb). In native foxes, more than half of the individual inbreeding (56%) was due to close relatedness between ancestors from less than 10 generations back in time, whereas in immigrant F1 the same fraction consisted of only 21% (Figure 3b). Immigrant F1 had 81% lower F ROH than native foxes when only considering inbreeding from recent ancestors (<10 generations, ROH >8 Mb; 95% CI [−0.21 to −0.06], df = 20, p = .002), 66% lower when considering inbreeding from ancestors up to 45 generations back (>2 Mb; 95% CI [−0.25 to −0.07], df = 20, p = .002) and 51% lower when considering inbreeding from ancestors deep back in history (<850 generations, ROH >100 kb; 95% CI [−0.24 to −0.07], df = 20, p = .002). Immigrant F2 and F3 did not differ from F1 (p > .07) or natives (p > .14) for any of the time frames.

| DISCUSS ION
We have linked inbreeding coefficients and genome-wide heterozygosity to individual ancestry and juvenile survival in 23 individuals of a threatened carnivore from a population that has survived at an extremely low size for over 100 years. We found extensive inbreeding in many individuals, with ROH stretching over nearly complete scaffolds (corresponding to chromosomal length), reflecting high inbreeding due to recent common ancestors within the past 10 generations. The substantially higher heterozygosity and lower inbreeding in immigrant F1 compared to native individuals demonstrate an effect of gene flow (Figure 3). Nevertheless, some long homozygous stretches (>8 Mb) were also found in the immigrant F1 genomes.
An explanation could be that parts of the immigrant ancestry is descended from Blåfjellet (Hasselgren et al., 2018), a population that is geographically close (approx. 120 km apart) and was up until recently connected to the study population.
We found that immigrant F2 and F3 had significantly lower heterozygosity than immigrant F1. However, we did not find any difference in F ROH between immigrant F2, F3 and the other groups but this is probably influenced by the low power due to small sample size and large individual variation in F ROH . One F2 and one F3 individual were however highly inbred, with 13%-17% of the genome contained within very long ROH (Figure 2; Table S1). Previous results from the same population revealed a consistent reduction in pedigree inbreeding (F PED ) from outcrossing with immigrants (Hasselgren et al., 2018). These results however imply a different pattern, as genomic inbreeding (F ROH ) in fact increased to similar levels as before immigration within two to three generations, at least in two out of four individuals. This in combination with no detected fitness elevation beyond immigrant F1 generation (Lotsander et al., 2021)  Similarly, in Isle Royale wolves, F ROH dropped drastically as a migrant male entered the population, but had increased to preimmigration levels within two generations (Robinson et al., 2019). Furthermore, in an isolated population of grey wolves in South Iberia, individuals recently hybridized with dogs showed long ROH and low genetic diversity, due to continued inbreeding between dog-wolf hybrids (Gómez-Sánchez et al., 2018). Our study thus contributes to a growing body of data implying that genetic rescue can be very short lived in wild carnivore populations (Gómez-Sánchez et al., 2018;Robinson et al., 2019;Saremi et al., 2019). In our case, immigration was comprised of only three individuals (of which two were full siblings) entering a small and fluctuating population. During recent years, close inbreeding within immigrant lineages has been recorded, potentially resulting in a genetic sweep (Lotsander et al., 2021). Additional gene flow into the population is needed to maintain beneficial effects of outcrossing.
We found that the inbreeding history corresponded well with the documented demographic history of the population. In native individuals, more than half (55%) of the F ROH originated from ancestors close in time (<10 generations back). This reflects the very recent re-establishment of the subpopulation and confirms the very rapid increase seen in F PED (Norén et al., 2016). The large fraction of inbreeding traced to recent ancestors, in combination with the consistent underestimation of F PED observed ( Figure S5), imply close relatedness between founders (whose kinship to each other is unknown). However, we also found that a substantial part of F ROH (27%) originated from common ancestors up to 45 generations back. This reflects the time period when human persecution caused the Scandinavian arctic fox to drastically decline in numbers (Lönnberg, 1927).
A key finding in this study is the negative relationship between inbreeding (F ROH ) and survival (Figure 4; Figure S4)

| Future perspectives
To better understand the genomic basis of inbreeding depression, it is essential to quantify mutational load in wild populations as well as the effect of deleterious mutations on protein coding genes and determine whether small populations are more likely to purge or accumulate mutational load. While some studies have addressed these questions (Dussex et al., in press;Kardos et al., 2018;Robinson et al., 2019;Xue et al., 2015), data documenting these processes are still scarce. Furthermore, from an applied conservation perspective, understanding the genomic basis of inbreeding depression and genetic rescue could potentially be an important tool for captive breeding programmes or choosing the most genetically suitable individuals for translocation into populations exposed to genomic erosion.

DATA AVA I L A B I L I T Y S TAT E M E N T
The de novo assembly has been deposited to Genbank under BioProject PRJNA704825 with the accession number JAFLEK000000000. Resequencing data (fastq files) are publicly available. ENA (PRJEB43377). Data set and r-code is available at Dryad (doi: https://doi.org/10.5061/dryad.gmsbc c2mn).