<figure><img src="image_1.png" /></figure>

HAL
open science

# Genomic footprints of recovery in the European bison

Tom Druet, Kamil Olénski, Laurence Flori, Amandine R Bertrand, Wanda Olech, Malgorzata Tokarska, Stanislaw Kaminski, Mathieu Gautier

<figure><img src="image_2.png" /></figure>

► To cite this version:

Tom Druet, Kamil Olénski, Laurence Flori, Amandine R Bertrand, Wanda Olech, et al. Ge-
nomic footprints of recovery in the European bison. Journal of Heredity, 2020, 111 (2), pp.194-203.
10.1093/jhered/esaa002. hal-02622499

HAL Id: hal-02622499

https://hal.inrae.fr/hal-02622499

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Genomic footprints of recovery in the European bison

Tom Druet<fnref n="1" />*, Kamil Olénski<fnref n="2" />, Laurence Flori<fnref n="3" />, Amandine R. Bertrand<fnref n="1" />, Wanda Olech<fnref n="4" />, Malgorzata

Tokarska<fnref n="5" />, Stanislaw Kaminski<fnref n="2" />, Mathieu Gautier<fnref n="6" />

Corresponding author:

Tom Druet, Unit of Animal Genomics, GIGA (B34 +1), Quartier Hôpital, Avenue de l'Hôpital, 11, B-4000 Liège, Belgium

Tel: +3243669172; Fax:+3243664151; e-mail: tom.druet@uliege.be

<footnote n="1">¹ Unit of Animal Genomics, GIGA-R & Faculty of Veterinary Medicine, University of Liège, Liège, Belgium</footnote>

<footnote n="2">² Department of Animal Genetics, University of Warmia and Mazury in Olsztyn, Olsztyn, Poland</footnote>

<footnote n="3">³ SELMET, INRA, CIRAD, Montpellier Supagro, Univ. Montpellier, Montpellier, France</footnote>

<footnote n="4">⁴ Department of Animal Genetics and Breeding, Warsaw University of Life Sciences – SGGW, Warsaw, Poland</footnote>

<footnote n="5">⁵ Mammal Research Institute, Polish Academy of Sciences, Białowieża, Poland</footnote>

<footnote n="6">⁶ INRA, UMR CBGP (INRA - IRD - Cirad - Montpellier SupAgro), Montferrier-sur-Lez, France</footnote>
**Abstract**

After extinction in the wild in the beginning of the twentieth century, the European bison has been successfully recovered in two distinct genetic lines from only twelve and seven captive founders. We here aimed at characterizing the levels of realized inbreeding in these two restored lines to provide empirical insights into the genomic footprints left by population recovery from a small number of founders. To that end, we genotyped 183 European bison born over the last 40 years with the Illumina BovineHD beadchip that contained 22,602 informative autosomal SNPs after data filtering. We then identified homozygous-by-descent (HBD) segments and classified them into different age-related classes relying on a model-based approach. As expected, we observed that the strong and recent founder effect experienced by the two lines resulted in very high levels of recent inbreeding and in the presence of long HBD tracks (up to 120 Mb). These long HBD tracks were associated with ancestors living approximately from 4 to 32 generations in the past suggesting that inbreeding accumulated over multiple generations after the bottleneck. The contribution to inbreeding of the most recent groups of ancestors were however found to be decreasing in both lines. In addition, comparison of Lowland individuals born at different time periods showed that the levels of inbreeding tended to stabilize, HBD segments being shorter in animals born more recently which indicates efficient control of inbreeding. Monitoring HBD segment lengths over generations may thus be viewed as a valuable genomic diagnostic tool for populations in conservation or recovery programs.

**Keywords:** European bison, inbreeding, recovery, ROH, homozygous-by-descent, RZooRoH
Introduction

The last wild European bison (Bison bonasus) were killed in 1919 in the Białowieża Forest in Poland, for the B. bonasus bonasus sub-species, and in 1927 in the Western Caucasus, for the B. bonasus caucasicus sub-species (Pucek, 1991). A restoration program was established in 1924 by the International Society of European Bison Protection. At that time, the population consisted of 54 captive animals living in European zoos and private breeding centers (Olech, 2009; Pucek, Belousova, Krasiński, & Olech, 2004; Tokarska, Pertoldi, Kowalczyk, & Perzanowski, 2011). This program led to the creation of two distinct genetic lines, namely the Lowland and the Lowland-Caucasian lines, still currently managed as separated populations. The reintroduction in the wild then started in 1952 when the first two bison were released for a trial period (Pucek et al., 2004). The Lowland line now includes close to 4,000 individuals, with approximately 3,100 individuals living in free-ranging herds, the largest herd being in the Białowieża forest where the recovery program started (Raczyński, 2018). Pedigree records showed that the Lowland line actually derives from only seven B. bonasus bonasus founders, two of which (the male Plebejer and the female Planta) contributing to approximately 80 % of the genome of contemporary individuals (Pucek et al., 2004; Tokarska et al., 2011). On the other hand, the Lowland-Caucasian line, that now consists of ca. 1,900 free-living (in the territory of Poland, Russia, Romania, Slovakia, Germany and Ukraine) and ca. 1,300 captive individuals (Raczyński, 2018), derives from eleven B. bonasus bonasus (including the seven founders of the Lowland line) and one B. bonasus caucasicus (named Kaukasus) founders. Both the Lowland and the Lowland-Caucasian lines thus represent striking examples of successful recovery after extinction in the wild. Yet, the strong founder effect resulted in a dramatic loss of genetic diversity (e.g., Hartl & Pucek, 1994; Tokarska et al., 2009) and high levels of inbreeding. For instance, the average inbreeding coefficient estimated from pedigree records was found equal to 0.50 in Lowland line individuals born after 1995 (Tokarska et al., 2011), and close to 0.30 in contemporary Lowland-Caucasian individuals (Olech, 2003).
At the genome level, individuals belonging to a population that experienced a strong and recent bottleneck have increased probability of carrying identical-by-descent (IBD) alleles (i.e., inherited twice from a common ancestor), including deleterious variants that may contribute to inbreeding depression and threaten the recovery process. Within individual genomes, IBD alleles are generally located in extended homozygous-by-descent (HBD) segments that appear as long stretches of homozygous genotypes called runs-of-homozygosity (ROH). The expected HBD segment lengths are inversely related to the number of generations to the common ancestor and their frequency to past effective population size and individual inbreeding coefficients. Therefore, the distribution of ROH lengths in individual genomes is informative on the history of their underlying population of origin (e.g., Ceballos, Joshi, Clark, Ramsay, & Wilson, 2018). Assessing the distribution of ROH within individual genomes has thus become popular to characterize genomic footprints of past demographic events in a wide range of model species including humans (Kirin et al., 2010; McQuillan et al., 2008; Pemberton et al., 2012), domesticated and intensively selected species such as livestock species and dogs (Bosse et al., 2012; Druet & Gautier, 2017; Ferenčaković et al., 2013) or wild populations (Gómez-Sánchez et al., 2018; Kardos et al., 2018; Kardos, Qvarnström, & Ellegren, 2017; Xue et al., 2015). Focusing on the Lowland and Lowland-Caucasian restored European bison lines, we here aimed at providing empirical insights into the dynamics of recovery after a strong founder effect. To that end, we characterized the distribution of HBD segments over recent times and the contribution of past ancestors to the levels of realized inbreeding in contemporary individuals. Capitalizing on the close phylogenetic relationship between the European bison and the cattle species (Gautier et al., 2016), we chose to rely on the available cattle genomic resources, regarding both genotyping assay (Wojciechowska et al., 2017) and genome assembly that provides a reasonable proxy to the European bison physical map (Wang et al., 2017). We built a large data set consisting of 183 European bison (154 Lowland and 29 Lowland-Caucasian individuals) born over the past 40 years that were genotyped with the commercial Illumina BovineHD beadchip that contained 22,602 autosomal and informative SNPs. We further relied on our recently developed hidden
Markov model (HMM), named ZooRoH (Druet & Gautier, 2017), to identify HBD segments and to evaluate the contribution of past generations of ancestors to inbreeding. Indeed, this model has been shown to remain robust to marker ascertainment bias (by accounting for population allele frequencies) and accurate with sparse genotyping data, i.e., with low and locally variable marker density (e.g., Druet & Gautier, 2017; Solé et al., 2017). We validated this approach by comparing the overall observed genetic diversity in HBD and non-HBD segments with Whole Genome Sequences (WGS) at a 10x coverage available for two of our genotyped individuals (Gautier et al., 2016).

## Material and Methods

**Genotyping Data.** A total of 189 European bison (154 Lowland individuals and 35 Lowland-Caucasian individuals) were genotyped with the Illumina BovineHD beadchip (Illumina, San Diego, CA) in three separate experiments. First, 147 Lowland individuals sampled between 1986 and 2011 in the Białowieża National Park (Poland) and five Lowland-Caucasian individuals were genotyped at the University of Warmia and Mazury (Olsztyn, Poland) (Kaminski, Olech, Olenski, Nowak, & Rusc, 2012). Second, seven Lowland individuals were sampled in 1991 in the Białowieża National Park (Gautier et al., 2016), and their genomic DNA was extracted, at that time, from blood samples (Jeanpierre, 1987) allowing their long-term preservation. They were further genotyped in 2015 at the Labogena laboratory for genetic analysis in animals (Jouy-en-Josas, France). Third, 30 Lowland-Caucasian individuals sampled mainly in the Caucasus Mountains in Russia by Taras Sipko between 2002 and 2012, were newly genotyped at the GIGA research center at the University of Liège (Belgium). At the University of Warmia and Mazury, genomic DNA was extracted from blood tissue using the BioSprint 96 DNA Blood kit (Qiagen) or from soft tissues using the DNeasy Blood and Tissue Kit (Qiagen) whereas at the University of Liège, it was extracted from tissue samples with the KingFisher extraction robot and using the NucleoMag Tissue Kit. In all three locations, SNP genotyping was conducted with the manual Illumina Infinium II Assay protocol and using standard procedures. We
only considered SNPs mapping to the 29 bovine autosomes on the Bovine UMD3.1 reference genome assembly (Liu et al., 2009) and discarded, using PLINK (Purcell et al., 2007), SNPs that i) were monomorphic, ii) displayed a call-rate <0.95, and iii) deviating from Hardy-Weinberg equilibrium ($p \le 0.001$) with a test statistic computed over the whole sample of individuals irrespective of their line of origin. Among the 23,901 remaining SNPs, 1,299 (5.43%) were also discarded because they showed strong genotype calling discrepancies (they were called heterozygous in a batch of 41 genotyped Lowland individuals while almost all the other individuals were called homozygous for the same allele) that suggested problems in the genotype cluster definition. Finally, six Lowland-Caucasian individuals with a call rate < 0.90 were excluded. The final analyzed data set consisted of genotyping data for 22,602 autosomal SNPs and 183 individuals (154 Lowland and 29 Lowland-Caucasian individuals).

**Population structure and genetic divergence across the two lines.** The structuring of genetic diversity across the 183 genotyped individuals was assessed by performing i) a multi-dimensional scaling (MDS) of the 183x183 matrix of genome-wide Allele Sharing Distances (ASD) using the --cluster and --mds-plot 4 options from PLINK v1.9 (Purcell et al., 2007); and ii) an unsupervised hierarchical clustering with K=2 clusters using default options with the ADMIXTURE v1.3.0 software (Alexander, Novembre, & Lange, 2009).

Pairwise $F_{ST}$ across the two European bison lines was estimated using the estimator by Weir and Cockerham (1984). We further estimated branch-specific differentiation of each line from their ancestral founding population using KimTree v2.0.1 (Clemente, Gautier, & Vitalis, 2018; Gautier & Vitalis, 2012) run with default options. KimTree is a Bayesian hierarchical model where the allele frequencies are modeled along each branch of a population tree (consisting only here of the two leaves underlying the two European bison lines) assuming a pure-drift model of divergence without migration using Kimura's time-dependent diffusion approximation for genetic drift (Kimura, 1964). The resulting estimated branch lengths correspond to divergence times $t$ measured on a diffusion
timescale (i.e., τ = t/N where t is the number of generations and N is the haploid effective population size along the whole period from the ancestral population). Note that the founding population of the two lines (root of the assumed tree) actually corresponds here to the 12 Lowland-Caucasian line founders (11 *B. bonasus bonasus* and 1 *B. bonasus caucasicus*) since these also include the seven founders of the Lowland line. Such a sub-sampling of founders may be interpreted under the KimTree model as an extra contribution to the divergence of the Lowland line from the modeled founding population thus leading to a (presumably very slight) overestimation of the Lowland line divergence time from its actual seven founders.

**Characterization of individual levels of inbreeding.** We relied on the ZooRoH model (Druet & Gautier, 2017) as implemented in the RZooRoH package (Bertrand, Kadri, Flori, Gautier, & Druet, 2019) to characterize individual inbreeding at both a global (genome-wide) and a local (locus-specific) scale. This model describes individual genomes as mosaics of HBD and non-HBD segments using a multiple HBD-classes Hidden Markov Model (HMM) to model the succession of observed individual SNP genotypes as a function of marker allele frequencies, the genotyping error rates and inter-marker genetic distances. Stretches of homozygous genotypes are indeed expected in HBD segments whereas genotypes are assumed to follow Hardy-Weinberg proportions in non-HBD classes. Each HBD class k is defined by its own rate R<sub>k</sub> that determines the expected length (equal to 1/R<sub>k</sub> Morgans) of the associated HBD segments and that is approximately equal to twice the number of generations to the common ancestor that transmitted the DNA segment (see Druet and Gautier (2017) and Bertrand et al. (2019) for a comprehensive description of the model). In other words, a given HBD class represents the contribution of ancestors living approximately 0.5*R<sub>k</sub>* generations in the past to the individual level of inbreeding, with small (large) values of R<sub>k</sub> corresponding to recent (old) ancestors and long (short) HBD segments. Here we considered for both lines a model with five HBD classes (with R<sub>k</sub> equal to 4, 8, 16, 32 and 64) and one non-HBD class whose rate was equal to 64
(see Druet and Gautier (2017)). HBD classes with higher rates were not fitted because their associated HBD segments were not captured with the available marker density. Thus, the model captured the contributions of ancestors living up to approximately 30 generations in the past. Allele frequencies were estimated for each line separately and a genotyping error rate of 0.1% was assumed. Finally, recombination distances were computed from physical distances (according to the UMD 3.1 Bos taurus assembly) assuming a cM to Mb ratio equal to 1 which provides a very good proxy given the marker density considered here, at least for cattle (e.g., Gautier et al., 2007; Kadri et al., 2016). As an output, the HMM provided the probabilities of being assigned to the different HBD classes at each marker position. The proportion of the genome associated with a specific HBD class (e.g., the proportion of the genome that was located within HBD segments of a certain length) was then obtained by averaging the corresponding HBD-class probabilities over all marker positions. As a result, the ZooRoH model allowed estimation of the contributions of each of the five age-based classes of ancestors (living approximately 2, 4, 8, 16 and 32 generations ago, respectively) to the overall level of individual inbreeding (i.e., the probability of a random locus to be in an HBD segment from any length). The latter was simply estimated by averaging over all the SNPs the sum of their five estimated HBD probabilities. Finally, to characterize the HBD segment length distribution, HBD segments were identified using the Viterbi algorithm that determines the sequence of HBD and non-HBD states with the highest likelihood.

To assess the efficiency of the ZooRoH model with the characteristics of the data set analyzed in the present study (marker spacing and frequencies), we performed a simulation study described in Supp. File 1 and similar to that we carried out in Druet and Gautier (2017).

**Validation of HBD segments identification with whole-genome sequencing data.**

To evaluate the accuracy of our characterization of HBD segments in European bison based on the BovineHD chip and the UMD3.1 Bos taurus genome assembly, we used the WGS data available for two of the genotyped individuals (Gautier et al., 2016). Indeed, the 183 genotyped European bison
included the two Lowland individuals (namely, BBO_3569 and BBO_3574 that were sampled in 1991 in the Białowieża forest) previously sequenced at an approximately 10x whole genome coverage and described in Gautier et al. (2016). In particular, Gautier et al. (2016) ran *mlrho* (Haubold, Pfaffelhuber, & Lynch, 2010) on each of these two individual WGS mapped onto the UMD3.1 bovine assembly (see Gautier et al., 2016 for details on data processing) to estimate the overall (scaled) mutation rate ($\theta$ = $4N_e\mu$). Relying on the output of these same analyses, i.e. the likelihoods $l_1^i(l_2^i)$ that each individual genomic position (considering only those positions covered by 3 to 30 reads, see Gautier et al. (2016)) is homozygous (heterozygous), we computed the probability that the individual is heterozygous at position *i* as
$$p_{het}^i = \frac{l_2^i \hat{\theta}}{l_1^i(1-\hat{\theta})+l_2^i \theta}$$

where $\hat{\theta} = 1.37 \times 10^{-3}$ and $\theta = 1.45 \times 10^{-3}$ for BBO_3569 and BBO_3574, respectively (Gautier et al., 2016). Genotypes were subsequently called heterozygous (or respectively homozygous) if their respective probability $p_{het}^i > 0.95$ (or respectively $p_{het}^i < 0.05$). Based on these (reference-free) genotype calls, we further computed the average heterozygosity levels as the proportion of called sites found as heterozygotes in 100kb bins centered on each SNP and compared it to its corresponding HBD probability (or HBD status) as estimated from the genotyping data (see above). Bins encompassing SNPs lying in true HBD segments were expected to present much lower heterozygosity levels. To obtain a genome-wide estimate of the expected heterozygosity in low and high heterozygosity bins, we computed the heterozygosity in all 100 kb non-overlapping windows over the whole individual genomes and fitted a mixture of two normal distributions with the mixtools R package (Benaglia, Chauveau, Hunter, & Young, 2009) to discriminate windows with low or high heterozygosity levels. Assuming the two normal distributions are associated with HBD (smallest estimated mean heterozygosity) and non-HBD (highest estimated mean heterozygosity) segments, this analysis provided an estimator of the genetic heterozygosity per bp in the two groups (HBD versus non-HBD segments) that accounted for noise introduced by errors in the sequencing data.
## Results

**Structuring of the European bison genetic diversity.** The plot of the first MDS coordinates of the 183 European bison from the Lowland (n=154) and the Lowland-Caucasian (n=29) lines clearly separated individuals according to their line of origin (Fig. 1A). Although the genotyping data set originated from three different sources (see Material and Methods), we did not observe an association between the genotyping laboratory and the first two MDS coordinates. Unsupervised hierarchical clustering with ADMIXTURE (Alexander et al., 2009) with K = 2 clusters provided an alternative quantification of the contribution of the five founders exclusive to the Lowland-Caucasian line (Fig. 1B). As expected, ancestry for a first cluster (represented in red in Fig. 1B) presumably representative of the seven founders that contributed to both lines was almost exclusive in all the 154 individuals from the Lowland line (0.9978 on average, 153 and 146 individuals having an estimated ancestry above respectively 0.95 and 0.99 for this cluster) and displayed varying contribution in the 29 Lowland-Caucasian individuals (from 0.00 to 0.85). The alternative and complementary cluster (represented in orange in Fig. 1B) being hence likely representative of the contribution of the five additional founders, including Kaukasus, a Caucasian European bison.

At the population level, the two lines were clearly differentiated, with a pairwise FST equal to 0.0735. The estimated divergence time (on a diffusion time scale) of the two lines from their founding population revealed that the Lowland line actually experienced increased drift ($\hat{t}_L$ = 0.139) compared to the Lowland-Caucasian line ($\hat{t}_{LC}$ = 0.0250). Assuming a generation time of six years (Gautier et al., 2016) and approximately tS=15 generations since the founding population (for all individuals), such divergence times would correspond to diploid effective population sizes (integrated over the whole recovery period) of 54 and 300 for the Lowland and the Lowland-Caucasian lines, respectively.
Characterization of individual levels of inbreeding

A simulation study was first performed based on the marker informativeness of our data, by using
the same genetic map and marker allele frequencies spectra (see Supp. File 1 for details). As
expected from our previous studies (Druet & Gautier, 2017; Solé et al., 2017), we found that
informativeness was high enough to accurately estimate the inbreeding levels in both lines. The
partitioning of the inbreeding in different HBD classes or the estimation of HBD segment lengths
may be more approximate.

The distribution of the individual inbreeding coefficients estimated in both lines is given in Fig. 2A. In
agreement with previous work based on pedigree records (Olech, 2003; Tokarska et al., 2011), the
average individual inbreeding levels in both the Lowland and the Lowland-Caucasian lines were
found particularly high and equal to 0.40 and 0.31 respectively. Such differences were also
concordant with the higher level of drift estimated for the Lowland line (see above). In addition, we
observed that the Lowland-Caucasian line displayed more variation in individual inbreeding levels
(inbreeding coefficients varying from 0.12 to 0.44, sd=0.09), including individuals with relatively low
levels of inbreeding, as compared to the Lowland line (inbreeding coefficient varying from 0.28 to
0.56, sd=0.05).

As shown in Fig. 2B, individual inbreeding levels were mostly associated with recent common
ancestors in both lines, i.e., mostly explained by long HBD segments that belonged to HBD classes
with low R<sub>k</sub>. In the Lowland line, HBD segments were concentrated in classes with R<sub>k</sub> equal to 8 and
16 whereas in the Lowland-Caucasian line, HBD classes with R<sub>k</sub> equal to 16 and 32 had the largest
contributions. In other words, inbreeding had a recent origin, the most recent common ancestors
tracing back to approximately 4 to 8 generations in the Lowland line and to 8 to 16 generations in
the Lowland-Caucasian line. These slight differences in the recent ancestral contributions to levels of
inbreeding between the two lines were also consistent with the genome-wide increased level of
genetic drift for the Lowland line we reported above.
Accordingly, at the individual level, Lowland individual genomes displayed a higher proportion of HBD segments associated with more recent common ancestors (HBD classes with lower R<sub>k</sub> values indicated by red and orange colors in Fig. 2C) compared to Lowland-Caucasian individuals. We even identified Lowland individuals with a high proportion of their genome associated with HBD classes belonging to the R<sub>k</sub> = 4 HBD class. The latter included extremely long HBD segments (25 cM on average) as illustrated in Supplementary Figure 1 for chromosome 10 on a subset of 30 Lowland and 29 Lowland-Caucasian individuals. Overall, 11,475 HBD segments (74.5 on average per individual) were identified in the Lowland line with lengths ranging from 0.2 to 123.7 Mb (median of 8.8 Mb and a mean of 12.4 Mb). In the Lowland-Caucasian line, 2,463 HBD segments were identified (84.9 on average per individual) with lengths ranging from 0.1 to 90.7 Mb (median of 5.1 Mb and a mean of 7.9 Mb). The Lowland and Lowland-Caucasian individuals carried on average 33.5 (13.5) and 21.7 (6.4) HBD segments longer than 10 Mb (20 Mb), respectively. Finally, several individuals from the Lowland line carried segments longer than 50 Mb (1.1 per individual on average). Such long segments were also observed in the Lowland-Caucasian line but to a lesser extent (0.6 per individual).

**Validation of HBD segment detection with Whole-Genome Sequence data.** WGS data available for two of the genotyped individuals, BBO_3574 and BBO_3569 (Gautier et al., 2016), allowed partitioning of their genomes into low and high heterozygosity 100 kb windows. Average heterozygosity in low heterozygosity windows, that represented 58.0% and 59.4% of the windows for individuals BBO_3574 and BBO_3569, respectively, was equal to 1.93x10<sup>-4</sup> and 1.95x10<sup>-4</sup>. Conversely, the average heterozygosity was one order of magnitude higher in high heterozygosity windows and equal to 2.47x10<sup>-3</sup> (for BBO_3574) and 2.41x10<sup>-3</sup> (for BBO_3569). We then estimated average genetic heterozygosity in 100 kb windows surrounding the positions (± 50 kb) from the 22,602 genotyped SNPs and summarized results according to HBD classification (Fig. 3). We
observed that regions surrounding SNPs that were declared as non-HBD presented much higher
levels of heterozygosity than those surrounding SNPs lying within estimated HBD segments. For
individual BBO_3569, the average heterozygosity was 1.81x10⁻³ around SNPs belonging to the non-
HBD class and 2.41x10⁻⁴ around SNPs belonging to an HBD class. Similarly, for individual BBO_3574,
the average heterozygosity was 1.89x10⁻³ for the non-HBD and 2.37x10⁻⁴ for HBD positions. Using
HBD probabilities, we found even more pronounced differences with an average heterozygosity of
2.42x10⁻⁴ and 2.23x10⁻⁴ around SNPs with an HBD probability higher than 0.99 compared to 2.47x10⁻⁵
³ and 2.62x10⁻³ around SNPs with an HBD probability lower than 0.01 for individuals BBO_3569 and
BBO_3574, respectively. These results indicate that SNPs declared HBD using the genotyping data
set are indeed in regions of strongly reduced heterozygosity (as estimated with the WGS data).
Moreover, the observed heterozygosity around SNPs estimated as lying in HBD segments matches
the heterozygosity observed in low heterozygosity windows identified in the sequence data.
However, some windows declared as HBD displayed high levels of heterozygosity (Fig. 3). Such
discrepancies might either be due to a lack of information in these regions in the genotyping data set
or to other errors such as local rearrangements between Bos taurus and European bison genomes.

Recent changes in the levels of inbreeding in the Lowland line

To provide insights into the recent changes in inbreeding levels, we regressed the birth years
available for 142 Lowland individuals on the relative contributions of each of the five HBD classes to
their overall level of inbreeding. To assess the significance of the estimated regression coefficients
and between group differences, we estimated the distribution of the corresponding test statistics
under the null hypothesis by performing 100,000 random permutations of the birth years. The
regression coefficients were non-significant for HBD classes with rates equal to 4, 32 or 64. However,
a significant decrease (p = 3.2 x 10⁻²) was estimated for the very recent HBD class with Rk = 8
whereas a significant increase (p = 1.6 x 10⁻³) was estimated for the class with Rk = 16, associated
with older ancestors (Fig. 4). As a result of these recent changes in the different contributions, the trend for inbreeding levels was non-significant (Fig. 4). Overall, these trends suggested that inbreeding in the current generations was associated with more distant ancestors (the common ancestors of individuals born more recently tracing back more generations in the past than individuals born a few decades ago) and an efficient control of inbreeding in the restored Lowland line. To further illustrate the changes over time, we compared two groups of individuals born at different time periods (e.g., separated by two generations). Compared to recent ancestors (n=39 individuals born before 1996), present day Lowland individuals (n=29 individuals born after 2005) showed lower proportions (p = 1.3 x 10⁻² for class with R<sub>k</sub> = 8) of the genome associated with recent HBD classes (with a rate < 16) but higher proportions (p = 7.6 x 10⁻² for class with R<sub>k</sub> = 16) for more ancient inbreeding (Supp. Fig 2). Also, the average levels of individual inbreeding increased only slightly from 0.40 for recent ancestors to 0.41 for present-day individuals in agreement with the results from the linear regression.

## Discussion

To provide empirical insights into the genomic footprints left by population recovery from a small number of founders, we characterized the levels of inbreeding in the two restored lines of the European Bison, an emblematic example of a successful coordinated conservation effort. Based on pedigree records, previous studies had already shown that modern European bison displayed high levels of inbreeding (e.g., Olech, 2003; Tokarska et al., 2011). However, such pedigree-based estimates correspond to expected rather than realized levels of inbreeding that can nowadays be more accurately characterized using genomic information (e.g., Kardos, Luikart, & Allendorf, 2015). We thus herein relied on genotypes for 183 European bison representative of the two lines obtained with a commercial high density SNP genotyping array developed for the closely related cattle (Bos taurus). As observed in previous studies that also used bovine SNP arrays to genotype European
bison (Kaminski et al., 2012; Wojciechowska et al., 2017), a high proportion of the interrogated SNPs were monomorphic and less than 5% of the SNPs (22,602 markers) were finally conserved after filtering. The resulting minor allele frequency (MAF) distribution of the analyzed SNPs was also clearly L-shaped (Supp. Fig. 3), only 66.8% of which showing a MAF > 0.05 over all the 183 individuals. The distribution of the spacing between consecutive markers of the selected map was not perfectly uniform (Supp. Fig. 4). These characteristics of our genotyping data set justified the use of a model-based approach, such as the ZooRoH model we recently developed, to provide genomic estimates of the individual levels of inbreeding via the identification of HBD segments. Model-based approaches indeed allow to explicitly account for the marker informativeness (i.e., allele frequencies), together with other information such as genetic distances, in a probabilistic framework to estimate HBD probabilities at each SNP position. We previously illustrated with cattle data (Solé et al., 2017) that even with less than 7,000 markers we could efficiently capture recent HBD segments associated with common ancestors tracing back to approximately 16 generations in the past (i.e., belonging to the HBD class with a rate equal to 32 that includes HBD segments with a 3 Mb expected length). For the present work, we first verified that the marker informativeness was high enough to estimate the inbreeding levels using a simulation study (Supp. File 1). Overall HBD probabilities were accurately estimated but in some situations, the partitioning in HBD classes was moderately shifted towards younger classes. This occurred with high inbreeding levels and for young HBD classes. Similarly, when considering an alternative modeling with a so-called 1R-model (Druet & Gautier, 2017) consisting of a single HBD class whose rate is estimated (rather than fixed), we found that it could be overestimated (which is equivalent to an underestimation of the length of HBD segments) when the genotyping error rates were high (e.g., 1%) and the inbreeding levels were smaller ($F \le 0.1$). Indeed, high genotyping error rates, as also errors in the genetic map (or in the genome assembly), would tend to cut HBD segments in smaller pieces (older HBD segments). Note however that we observed rather long HBD segments in both lines, indicating that such errors appeared to have a relatively small impact in our study. Conversely, with high inbreeding levels, as
we observed here, several consecutive HBD segments may incorrectly be considered as a single
longer HBD segment, if marker informativeness is not high enough, leading to an underestimation of
the age of the segments. Overall, it is important to emphasize that the estimation of the HBD
segments age remains approximate, and most particularly, when the marker informativeness is low.
The effect of these potential approximations is however likely to be of similar magnitude across the
compared groups of individuals (i.e., Lowland versus Lowland-Caucasian and older versus recent)
and they may thus have little impact on the observed trends supporting our main conclusions (see
also Supp. File 1 for more details). We further confirmed the effectiveness of our inference of HBD
segments using WGS data available for two genotyped Lowland European bison (Gautier et al.,
2016). The HBD segments identified in these two individuals based on the genotyping data displayed
a clear reduction in heterozygosity in the sequence data when compared to the rest of the genome
(about one order of magnitude). As a matter of comparison, it should be noted that the observed
average heterozygosity in these HBD segments was half as high as the threshold (5x10⁻⁴) considered
by Gómez-Sánchez et al. (2018) to call HBD segments in WGS data. Nevertheless, if this approach
based on WGS data was effective to validate that positions declared as HBD with genotyping data
were truly HBD, it provided less information regarding the accuracy of HBD segment length
estimation. Sequencing a larger number of individuals (e.g., a dozen) may be valuable to provide a
detailed and accurate picture of HBD segments distribution in the European genome. Yet, the
marker informativeness provided by the bovine SNP high density assay appeared sufficient to
capture the genomic footprints left by the recovery of this species that experienced a recent
bottleneck.

Indeed, we herein confirmed that both European bison lines displayed high levels of realized
inbreeding associated with recent ancestors, consistent with the pedigree-based expectations. These
levels were also in agreement with the respective history of both lines, with higher levels of
inbreeding in the Lowland line (equal to 0.40 on average) that was recovered from fewer founders
than the Lowland-Caucasian line (equal to 0.31 on average). These levels of inbreeding observed
several generations after restoration thus remained particularly high, even higher than values
estimated with a similar approach and with 50K marker panels for intensively selected cattle
populations (Solé et al., 2017) such as the Belgian Blue cattle breed (0.09) and for some sheep
populations (Druet & Gautier, 2017), like the Rambouillet (0.12) selected breed or the Soay Island
population (0.22) that experienced a strong founder event. The levels of inbreeding we reported in
the restored European bison lines were actually more comparable to those observed in highly inbred
dog breeds (Druet & Gautier, 2017) and other species that experienced a strong population decline
(Gómez-Sánchez et al., 2018; Kardos et al., 2018; Xue et al., 2015). For instance, one third of the
genome from the mountain and eastern lowland gorillas was reported to be in HBD tracks based on
whole-genome sequence (WGS) data analysis (Xue et al., 2015). Based on ROH analysis of WGS data,
Kardos et al. (2018) estimated that 27 % of the genome of an isolated Scandinavian wolf population
was associated with long HBD segments originating from recent ancestors present less than 10
generations ago. Similarly, the estimated inbreeding coefficient from one of the last wolves from the
Sierra Morena mountain range in Spain was as high as 0.42 (Gómez-Sánchez et al., 2018). As in these
extreme examples from highly inbred populations, long homozygous segments were found to be
extremely frequent in the genome of the contemporary European bison. On average, individuals
from the Lowland line carried more than 30 HBD segments longer than 10 Mb and one segment 50
Mb or longer. Nevertheless, the HBD classes with rates of 8 and 16, and that included also shorter
HBD segments (12.5 Mb down to 6.25 Mb on average), had the major contribution to inbreeding in
the restored Lowland line. The main contributing groups of ancestors traced back to four to eight
generations in the past in the Lowland line and eight to sixteen generations in the past in the
Lowland-Caucasian line. Although these estimated ages need to be interpreted with caution, these
would correspond to ancestors living a few generations after the bottleneck that occurred
approximately 15 generations ago (assuming a 6 year generation interval and an average birth date
in 2010 for all the individuals). Ideally, the Lowland-Caucasian individuals that were sampled
between 2002 and 2012 should be compared to Lowland individuals sampled during the same
period. Interestingly, the observed differences between the Lowland-Caucasian and the Lowland individuals (lower inbreeding levels and shorted HBD segments) were still true when comparisons were done using the present day Lowland individuals only.

The recovery process left distinct footprints in the genome of the European bison. First, the contributions of the different HBD classes to the overall inbreeding varied as a function of their birth date. Indeed, while the levels of inbreeding were only slightly higher in individuals born after 2005 compared to those born before 1996, the HBD segments tended to be shorter in more recent individuals. Such a pattern is actually expected in an expanding population after a bottleneck since the probability to mate with a close relative decreases as the number of parents increases. As a result, the expected time to the ancestor underlying HBD segments increases. However, the individual inbreeding levels don't decrease because if we trace back more generations in the past we finally find a common ancestor due to the extremely small founding population. In general, an observed decrease through time of the overall HBD segment lengths may be regarded as positive as this would suggest a limited additional inbreeding and an optimal mixing of the founder haplotypes.

The stabilization of the levels of inbreeding and the reduction of HBD segment lengths that we observed in the European bison lines indicate that the restoration plan has been successful to control inbreeding although genotyping (at high density) more individuals representing evenly the different time periods would be valuable to provide a refined picture. Monitoring HBD segment lengths over generations may hence be viewed as a valuable genomic diagnostic tool for conservation or recovery programs. The recovery process can also be evaluated using contemporary individuals only. Indeed, we observed in the current population that the contribution from the most recent group of ancestors was smaller than that of more ancient groups. Ancestors living a few generations after the bottleneck had the highest contribution. After that peak, the contributions to inbreeding decreased indicating that the population is expanding and that inbreeding is mostly associated to the period surrounding the bottleneck. Hence, genotyping contemporary individuals and using model-based approaches such as the ZooRoH model to partition their genomes into
various age-based related classes of HBD segments may represent a valuable strategy even when
samples of different ages are not available.

Beyond the distribution of HBD segments in the genome, a functional characterization of the genomic footprints left by the recovery process would be needed to better understand why the recovery was successful and why the modern European bison lines show modest inbreeding depression effects (Raczyński, 2018; Tokarska et al., 2011). Severe bottlenecks can indeed result in the accumulation of moderately deleterious variants but in extreme cases the most severe deleterious variants can be purged out of the population as was observed in mountain gorillas (Xue et al., 2015). Purging of deleterious alleles is indeed sometimes more efficient in small populations (Lynch & Walsh, 1998). Comparing the amount of deleterious variants in current and historical populations, as was done for instance in Eastern Gorillas (van der Valk, Díez-del-Molino, Marques-Bonet, Guschanski, & Dalén, 2019), may be informative to study the evolution of the genetic load in the European bison.

Funding

This work was supported by the Fonds de la Recherche Scientifique - FNRS (F.R.S.-FNRS) under Grants J.0134.16 and J.0154.18. The genotyping of 152 bison samples was funded by the University of Warmia and Mazury (Olsztyn, Poland) under the grant n° 0105-804.

Acknowledgements

The authors wish to thank Anna Madeyska-Lewandowska (Private Veterinary School, Warsaw, Poland), Hubert Levéziel and Cécile Grohs who provided samples for the seven Lowland individuals genotyped at INRA (AFROSEQ program supported by the INRA Animal Genetics Department) as part
of a collaborative research project on the characterization of casein variants. We would like to thank
all the anonymous reviewers for valuable comments on our work. Tom Druet is Senior Research
Associate from the F.R.S.-FNRS. We used the supercomputing facilities of the "Consortium
d'Equipements en Calcul Intensif en Fédération Wallonie-Bruxelles" (CECI), funded by the F.R.S.-
F.N.R.S.

Data Accessibility

The genotyping data used in the present study were deposited in Dryad. The sequencing data was
downloaded from the Sequence Read Archive repository (accession number SRP070526).
References

Alexander, D. H., Novembre, J., & Lange, K. (2009). Fast model-based estimation of ancestry in unrelated individuals. *Genome Research*.

Benaglia, T., Chauveau, D., Hunter, D., & Young, D. (2009). *mixtools: An R package for analyzing finite mixture models*. *Journal of Statistical Software*, **32**(6), 1–29.

Bertrand, A., Kadri, N. K., Flori, L., Gautier, M., & Druet, T. (2019). *RZooRoH: An R package to characterize individual genomic autozygosity and identify homozygous-by-descent segments*. *Methods in Ecology and Evolution*.

Bosse, M., Megens, H.-J., Madsen, O., Paudel, Y., Frantz, L. A. F., Schook, L. B., ... Groenen, M. A. M. (2012). Regions of homozygosity in the porcine genome: Consequence of demography and the recombination landscape. *PLoS Genetics*, **8**(11), e1003100.
https://doi.org/10.1371/journal.pgen.1003100

Ceballos, F. C., Joshi, P. K., Clark, D. W., Ramsay, M., & Wilson, J. F. (2018). Runs of homozygosity: Windows into population history and trait architecture. *Nature Reviews Genetics*, **19**(4), 220.

Clemente, F., Gautier, M., & Vitalis, R. (2018). Inferring sex-specific demographic history from SNP data. *PLoS Genetics*, **14**(1), e1007191.

Druet, T., & Gautier, M. (2017). A model-based approach to characterize individual inbreeding at both global and local genomic scales. *Molecular Ecology*, **26**(20), 5820–5841.

Ferenčaković, M., Hamzić, E., Gredler, B., Solberg, T. R., Klemetsdal, G., Curik, I., & Sölkner, J. (2013). Estimates of autozygosity derived from runs of homozygosity: Empirical evidence from selected cattle populations. *Journal of Animal Breeding and Genetics = Zeitschrift Fur Tierzuchtung Und Zuchtungsbiologie*, **130**(4), 286–293. https://doi.org/10.1111/jbg.12012

Gautier, M., Faraut, T., Moazami-Goudarzi, K., Navratil, V., Foglio, M., Grohs, C., ... Lathrop, G. M. (2007). Genetic and haplotypic structure in 14 European and African cattle breeds. *Genetics*, **177**(2), 1059–1070.

Gautier, M., Moazami-Goudarzi, K., Levéziel, H., Parinello, H., Grohs, C., Rialle, S., ... Flori, L. (2016). Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences. *Molecular Biology and Evolution*, **33**(11), 2801–2814.

Gautier, M., & Vitalis, R. (2012). Inferring population histories using genome-wide allele frequency data. *Molecular Biology and Evolution*, **30**(3), 654–668.

Gómez-Sánchez, D., Olalde, I., Sastre, N., Enseñat, C., Carrasco, R., Marques-Bonet, T., ... Ramírez, O. (2018). On the path to extinction: Inbreeding and admixture in a declining grey wolf population. *Molecular Ecology*, **27**(18), 3599–3612.

Hartl, G. B., & Pucek, Z. (1994). Genetic depletion in the European bison (Bison bonasus) and the significance of electrophoretic heterozygosity for conservation. *Conservation Biology*, **8**(1), 167–174.
Habold, B., Pfaffelhuber, P., & Lynch, M. (2010). *MIrho*–a program for estimating the population mutation and recombination rates from shotgun-sequenced diploid genomes. *Molecular Ecology*, **19**, 277–284.

Jeanpierre, M. (1987). A rapid method for the purification of DNA from blood. *Nucleic Acids Research*, **15**(22), 9611–9611.

Kadri, N. K., Harland, C., Faux, P., Cambisano, N., Karim, L., Coppieters, W., ... Boichard, D. (2016). Coding and noncoding variants in HFM1, MLH3, MSH4, MSH5, RNF212, and RNF212B affect recombination rate in cattle. *Genome Research*, **26**(10), 1323–1332.

Kaminski, S., Olech, W., Olenski, K., Nowak, Z., & Rusc, A. (2012). Single nucleotide polymorphisms between two lines of European bison (Bison bonasus) detected by the use of Illumina Bovine 50 K BeadChip. *Conservation Genetics Resources*, **4**(2), 311–314.

Kardos, M., Åkesson, M., Fountain, T., Flagstad, Ø., Liberg, O., Olason, P., ... Ellegren, H. (2018). Genomic consequences of intensive inbreeding in an isolated wolf population. *Nature Ecology & Evolution*, **2**(1), 124.

Kardos, M., Luikart, G., & Allendorf, F. W. (2015). Measuring individual inbreeding in the age of genomics: Marker-based measures are better than pedigrees. *Heredity*, **115**(1), 63.

Kardos, M., Qvarnström, A., & Ellegren, H. (2017). Inferring individual inbreeding and demographic history from segments of identity by descent in Ficedula flycatcher genome sequences. *Genetics*, **205**(3), 1319–1334.

Kimura, M. (1964). Diffusion models in population genetics. *Journal of Applied Probability*, **1**(2), 177–232.

Kirin, M., McQuillan, R., Franklin, C. S., Campbell, H., McKeigue, P. M., & Wilson, J. F. (2010). Genomic runs of homozygosity record population history and consanguinity. *PloS One*, **5**(11), e13996.

Liu, Y., Qin, X., Song, X.-Z. H., Jiang, H., Shen, Y., Durbin, K. J., ... Ren, Y. (2009). Bos taurus genome assembly. *BMC Genomics*, **10**(1), 180.

Lynch, M., & Walsh, B. (1998). *Genetics and analysis of quantitative traits* (Vol. 1). Sinauer
Sunderland, MA.

McQuillan, R., Leutenegger, A.-L., Abdel-Rahman, R., Franklin, C. S., Pericic, M., Barac-Lauc, L., ...
Wilson, J. F. (2008). Runs of homozygosity in European populations. *American Journal of Human
Genetics*, **83**(3), 359–372. https://doi.org/10.1016/j.ajhg.2008.08.007

Olech, W. (2003). *Wpływ inbredu osobniczego i inbredu matki na przezywalnosc cielat zubra (Bison bonasus)*. Wydawn. SGGW.

Olech, W. (2009). The changes of founders' number and their contribution to the European bison population during 80 years of species' restitution. *European Bison Conservation Newsletter*, **2**, 54–60.
Pemberton, T. J., Absher, D., Feldman, M. W., Myers, R. M., Rosenberg, N. A., & Li, J. Z. (2012).
Genomic patterns of homozygosity in worldwide human populations. *The American Journal of
Human Genetics*, 91(2), 275–292.

Pucek, Z. (1991). *History of the European bison and problems of its protection and management*. 19–39. Swiat Press.

Pucek, Z., Belousova, I. P., Krasiński, Z. A., & Olech, W. (2004). *European bison: Status survey and conservation action plan*. IUCN Gland, Switzerland and Cambridge.

Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A. R., Bender, D., ... Sham, P. C. (2007).
PLINK: A tool set for whole-genome association and population-based linkage analyses. *American
Journal of Human Genetics*, 81(3), 559–575. https://doi.org/10.1086/519795

Raczyński, J. (2018). *European Bison Pedigree Book 2017*. Białowieża National Park.

Solé, M., Gori, A.-S., Faux, P., Bertrand, A., Farnir, F., Gautier, M., & Druet, T. (2017). Age-based partitioning of individual genomic inbreeding levels in Belgian Blue cattle. *Genetics Selection Evolution*, 49(1), 92.

Tokarska, M., Marshall, T., Kowalczyk, R., Wojcik, J., Pertoldi, C., Kristensen, T., ... Bendixen, C. (2009). Effectiveness of microsatellite and SNP markers for parentage and identity analysis in species with low genetic diversity: The case of European bison. *Heredity*, 103(4), 326.

Tokarska, M., Pertoldi, C., Kowalczyk, R., & Perzanowski, K. (2011). Genetic status of the European bison Bison bonasus after extinction in the wild and subsequent recovery. *Mammal Review*, 41(2), 151–162.

van der Valk, T., Díez-del-Molino, D., Marques-Bonet, T., Guschanski, K., & Dalén, L. (2019). Historical Genomes Reveal the Genomic Consequences of Recent Population Decline in Eastern Gorillas. *Current Biology*, 29(1), 165–170.

Wang, K., Wang, L., Lenstra, J. A., Jian, J., Yang, Y., Hu, Q., ... Du, Z. (2017). The genome sequence of the wisent (Bison bonasus). *GigaScience*, 6(4), 1–5.

Weir, B. S., & Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. *Evolution*, 38(6), 1358–1370.

Wojciechowska, M., Nowak, Z., Gurgul, A., Olech, W., Drobik, W., & Szmatoła, T. (2017). Panel of informative SNP markers for two genetic lines of European bison: Lowland and Lowland-Caucasian. *Animal Biodiversity and Conservation*, 40(1), 17–25.

Xue, Y., Prado-Martinez, J., Sudmant, P. H., Narasimhan, V., Ayub, Q., Szpak, M., ... Cooper, D. N. (2015). Mountain gorilla genomes reveal the impact of long-term population decline and inbreeding. *Science*, 348(6231), 242–245.
## Figure legends

**Figure 1. Population structure across the European bison samples.** A) Multidimensional scaling analysis of the matrix of pairwise individual Allele Sharing Distance averaged over 22,602 SNPs. All the 183 genotyped individuals are plotted along the two first coordinates. Individuals from the Lowland and the Lowland-Caucasian lines are plotted in red and orange, respectively. B) Unsupervised hierarchical clustering from the ADMIXTURE program assuming K=2 clusters. Individuals are sorted according to the first coordinate of the MDS analysis. As a result, individuals from the Lowland-Caucasian line are on the right of the dashed line. See online version for full colors.

**Figure 2. Characterization of individual inbreeding levels.** A) Distributions of the inbreeding coefficients for the Lowland (154 individuals) and the Lowland-Caucasian (29 individuals) lines. Each distribution is represented as a violin plot, the inner part of which being equivalent to a boxplot (with median and interquantile range). The distribution clouds are kernel density plots estimated using default options of the vioplot function from the vioplot R-package. B) Proportion of individual genomes associated with different HBD classes. Each bar represents the average genome-wide HBD probability of an HBD class (in red and orange for the Lowland and Lowland-Caucasian lines, respectively), i.e., the average proportion of the genome associated with that HBD class for individuals from the given line. Individual proportions of the genome associated with a specific HBD class are obtained by averaging the corresponding HBD-class probabilities over all marker positions. The plotted values are then obtained by averaging individual values over all individuals. HBD classes with lower rates $R_k$ correspond to more recent common ancestor (living approximately 0.5 * $R_k$
generations in the past). C) Partitioning of individual genomes in different HBD classes. Each bar represents an individual and the height of the different stacks represents the proportion of the genome associated with the HBD class of the corresponding color (the colors used for the different rates of HBD classes are reported in the legend). The total height represents the overall level of inbreeding. See online version for full colors.

**Figure 3.** Average heterozygosity estimated around 22,602 genotyped SNPs as a function of their HBD classification estimated in two Lowland individuals (BBO_3569 and BBO_3574). The heterozygosity was estimated from WGS data in 100 kb regions surrounding each SNP (± 50 kb). The SNP HBD classifications (within an HBD segment or not) were obtained by running the multiple HBD classes model on the two individuals (based on allele frequencies estimated on the whole Lowland sample). For BBO_3569 (resp. BBO_3574), 9707 (resp. 8812) SNPs were classified as HBD and 12895 (resp. 13790) were classified as non-HBD. Results are presented as boxplots of the average heterozygosity observed in each class (HBD / non_HBD) for both individuals separately.

**Figure 4.** Changes of the inbreeding coefficients (upper panel) and of the proportions of the genome in HBD classes with rates $R_k = 8$ (middle panel) and $R_k = 16$ (lower panel) as a function of the birth year. The inbreeding coefficients and the contributions of the different HBD classes were estimated with the multiple HBD classes model on the 154 Lowland individuals and using 22,602 SNPs. The linear regressions plots were realized with the basicTrendline R package using default options. See online version for full colors.
<figure><img src="image_4.png" /><figcaption>Figure 1</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 2</figcaption></figure>

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>Figure 3</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 4</figcaption></figure>