domain: sdv info sde

text:
Tom Druet
email: tom.druet@uliege.be

Kamil Olénski

Laurence Flori

Amandine R Bertrand

Wanda Olech

Malgorzata Tokarska

Stanislaw Kaminski

Mathieu Gautier

Genomic footprints of recovery in the European bison

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

à 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.

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 [START_REF] Pucek | History of the European bison and problems of its protection and management[END_REF]. 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 [START_REF] Olech | The changes of founders' number and their contribution to the European bison population during 80 years of species' restitution[END_REF][START_REF] Pucek | European bison: Status survey and conservation action plan[END_REF][START_REF] Tokarska | Genetic status of the European bison Bison bonasus after extinction in the wild and subsequent recovery[END_REF]. 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 [START_REF] Pucek | European bison: Status survey and conservation action plan[END_REF]. 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 [START_REF] Raczyński | European Bison Pedigree Book[END_REF]. 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 [START_REF] Pucek | European bison: Status survey and conservation action plan[END_REF][START_REF] Tokarska | Genetic status of the European bison Bison bonasus after extinction in the wild and subsequent recovery[END_REF]. 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 [START_REF] Raczyński | European Bison Pedigree Book[END_REF], 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., [START_REF] Hartl | Genetic depletion in the European bison (Bison bonasus) and the significance of electrophoretic heterozygosity for conservation[END_REF][START_REF] Tokarska | Effectiveness of microsatellite and SNP markers for parentage and identity analysis in species with low genetic diversity: The case of European bison[END_REF] 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 [START_REF] Tokarska | Genetic status of the European bison Bison bonasus after extinction in the wild and subsequent recovery[END_REF], and close to 0.30 in contemporary Lowland-Caucasian individuals [START_REF] Olech | Wpływ inbredu osobniczego i inbredu matki na przezywalnosc cielat zubra (Bison bonasus)[END_REF]. 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., [START_REF] Ceballos | Runs of homozygosity: Windows into population history and trait architecture[END_REF]. 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 [START_REF] Kirin | Genomic runs of homozygosity record population history and consanguinity[END_REF][START_REF] Mcquillan | Runs of homozygosity in European populations[END_REF][START_REF] Pemberton | patterns of homozygosity in worldwide human populations[END_REF], domesticated and intensively selected species such as livestock species and dogs [START_REF] Bosse | Regions of homozygosity in the porcine genome: Consequence of demography and the recombination landscape[END_REF]Druet & Gautier, 2017;[START_REF] Ferenčaković | Estimates of autozygosity derived from runs of homozygosity: Empirical evidence from selected cattle populations[END_REF] or wild populations [START_REF] Gómez-Sánchez | On the path to extinction: Inbreeding and admixture in a declining grey wolf population[END_REF][START_REF] Kardos | Genomic consequences of intensive inbreeding in an isolated wolf population[END_REF][START_REF] Kardos | Inferring individual inbreeding and demographic history from segments of identity by descent in Ficedula flycatcher genome sequences[END_REF][START_REF] Xue | Mountain gorilla genomes reveal the impact of long-term population decline and inbreeding[END_REF]. 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 [START_REF] Gautier | Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences[END_REF], we chose to rely on the available cattle genomic resources, regarding both genotyping assay [START_REF] Wojciechowska | Panel of informative SNP markers for two genetic lines of European bison: Lowland and Lowland-Caucasian[END_REF] and genome assembly that provides a reasonable proxy to the European bison physical map [START_REF] Wang | The genome sequence of the wisent (Bison bonasus)[END_REF]. We built a large data set (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;[START_REF] Solé | Age-based partitioning of individual genomic inbreeding levels in Belgian Blue cattle[END_REF]. 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 [START_REF] Gautier | Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences[END_REF]. Second, seven Lowland individuals were sampled in 1991 in the Białowieża National Park [START_REF] Gautier | Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences[END_REF], and their genomic DNA was extracted, at that time, from blood samples [START_REF] Jeanpierre | A rapid method for the purification of DNA from blood[END_REF] [START_REF] Liu | Bos taurus genome assembly[END_REF] and discarded, using PLINK [START_REF] Purcell | PLINK: A tool set for whole-genome association and population-based linkage analyses[END_REF], SNPs that i) were 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 [START_REF] Purcell | PLINK: A tool set for whole-genome association and population-based linkage analyses[END_REF]; and ii) an unsupervised hierarchical clustering with K=2 clusters using default options with the ADMIXTURE v1.3.0 software [START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF].

Material and Methods

Genotyping

Pairwise FST across the two European bison lines was estimated using the estimator by [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF]. We further estimated branch-specific differentiation of each line from their ancestral founding population using KimTree v2.0.1 [START_REF] Clemente | Inferring sex-specific demographic history from SNP data[END_REF][START_REF] Gautier | Inferring population histories using genome-wide allele frequency data[END_REF] 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 [START_REF] Kimura | Diffusion models in population genetics[END_REF].

The resulting estimated branch lengths correspond to divergence times τ measured on a diffusion Characterization of individual levels of inbreeding. We relied on the ZooRoH model (Druet & Gautier, 2017) as implemented in the RZooRoH package [START_REF] Bertrand | RZooRoH: An R package to characterize individual genomic autozygosity and identify homozygous-by-descent segments[END_REF] (see [START_REF] Druet | A model-based approach to characterize individual inbreeding at both global and local genomic scales[END_REF]). 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., [START_REF] Gautier | Genetic and haplotypic structure in 14 European and African cattle breeds[END_REF][START_REF] Kadri | Coding and noncoding variants in HFM1, MLH3, MSH4, MSH5, RNF212, and RNF212B affect recombination rate in cattle[END_REF]. 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 [START_REF] Druet | A model-based approach to characterize individual inbreeding at both global and local genomic scales[END_REF].

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 [START_REF] Gautier | Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences[END_REF]. Indeed, the 183 genotyped European bison

Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002

A c c e p t e d M a n u s c r i p t 9 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 [START_REF] Gautier | Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences[END_REF]. In particular, [START_REF] Gautier | Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences[END_REF] ran mlrho [START_REF] Haubold | MlRho-a program for estimating the population mutation and recombination rates from shotgun-sequenced diploid genomes[END_REF] on each of these two individual WGS mapped onto the UMD3.1 bovine assembly (see [START_REF] Gautier | Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences[END_REF] for details on data processing) to estimate the overall (scaled) mutation rate (θ = 4Neµ). Relying on the output of these same analyses, i.e. the likelihoods 𝑙 1 𝑖 (𝑙 2 𝑖 )

that each individual genomic position (considering only those positions covered by 3 to 30 reads, see [START_REF] Gautier | Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences[END_REF]) is homozygous (heterozygous), we computed the probability that the individual is heterozygous at position i as

𝑝 ℎ𝑒𝑡 𝑖 = 𝑙 2 𝑖 𝜃 l1 𝑖 (1-𝜃 ̂)+𝑙 2 𝑖 𝜃 ̂
where 𝜃 ̂ = 1.37x10 -3 and 𝜃 ̂ = 1.45x10 -3 for BBO_3569 and BBO_3574, respectively [START_REF] Gautier | Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences[END_REF]. Genotypes were subsequently called heterozygous (or respectively homozygous) if their respective probability 𝑝 ℎ𝑒𝑡 𝑖 > 0.95 (or respectively 𝑝 ℎ𝑒𝑡 𝑖 < 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 [START_REF] Benaglia | mixtools: An R package for analyzing finite mixture models[END_REF] 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.

Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski

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;[START_REF] Solé | Age-based partitioning of individual genomic inbreeding levels in Belgian Blue cattle[END_REF], 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 [START_REF] Olech | Wpływ inbredu osobniczego i inbredu matki na przezywalnosc cielat zubra (Bison bonasus)[END_REF][START_REF] Tokarska | Genetic status of the European bison Bison bonasus after extinction in the wild and subsequent recovery[END_REF], 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 Validation of HBD segment detection with Whole-Genome Sequence data. WGS data available for two of the genotyped individuals, BBO_3574 and BBO_3569 [START_REF] Gautier | Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences[END_REF], 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 -4 and 1.95x10 -4 .

Conversely, the average heterozygosity was one order of magnitude higher in high heterozygosity windows and equal to 2.47x10 -3 (for BBO_3574) and 2.41x10 -3 (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

Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002

A c c e p t e d M a n u s c r i p t 13 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 -3 around SNPs belonging to the non-HBD class and 2.41x10 -4 around SNPs belonging to an HBD class. Similarly, for individual BBO_3574, the average heterozygosity was 1.89x10 -3 for the non-HBD and 2.37x10 -4 for HBD positions. Using HBD probabilities, we found even more pronounced differences with an average heterozygosity of 2.42x10 -4 and 2.23x10 -4 around SNPs with an HBD probability higher than 0.99 compared to 2.47x10 - 3 and 2.62x10 -3 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 -2 ) was estimated for the very recent HBD class with Rk = 8 whereas a significant increase (p = 1.6 x 10 -3 ) was estimated for the class with Rk = 16, associated

Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002

A c c e p t e d M a n u s c r i p t 14 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

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., [START_REF] Olech | Wpływ inbredu osobniczego i inbredu matki na przezywalnosc cielat zubra (Bison bonasus)[END_REF][START_REF] Tokarska | Genetic status of the European bison Bison bonasus after extinction in the wild and subsequent recovery[END_REF]. 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., [START_REF] Kardos | Measuring individual inbreeding in the age of genomics: Marker-based measures are better than pedigrees[END_REF].

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

Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020) [START_REF] Kaminski | Single nucleotide polymorphisms between two lines of European bison (Bison bonasus) detected by the use of Illumina Bovine 50 K BeadChip[END_REF][START_REF] Wojciechowska | Panel of informative SNP markers for two genetic lines of European bison: Lowland and Lowland-Caucasian[END_REF], 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 [START_REF] Solé | Age-based partitioning of individual genomic inbreeding levels in Belgian Blue cattle[END_REF] 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≤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

Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. ( 2020 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 [START_REF] Gautier | Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences[END_REF]. 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 -4 ) considered by [START_REF] Gómez-Sánchez | On the path to extinction: Inbreeding and admixture in a declining grey wolf population[END_REF] 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

Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002

A c c e p t e d M a n u s c r i p t 17 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 [START_REF] Solé | Age-based partitioning of individual genomic inbreeding levels in Belgian Blue cattle[END_REF] 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 [START_REF] Gómez-Sánchez | On the path to extinction: Inbreeding and admixture in a declining grey wolf population[END_REF][START_REF] Kardos | Genomic consequences of intensive inbreeding in an isolated wolf population[END_REF][START_REF] Xue | Mountain gorilla genomes reveal the impact of long-term population decline and inbreeding[END_REF]. 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 [START_REF] Xue | Mountain gorilla genomes reveal the impact of long-term population decline and inbreeding[END_REF]. Based on ROH analysis of WGS data, 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 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 [START_REF] Raczyński | European Bison Pedigree Book[END_REF][START_REF] Tokarska | Genetic status of the European bison Bison bonasus after extinction in the wild and subsequent recovery[END_REF]. 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 [START_REF] Xue | Mountain gorilla genomes reveal the impact of long-term population decline and inbreeding[END_REF]. Purging of deleterious alleles is indeed sometimes more efficient in small populations [START_REF] Lynch | Genetics and analysis of quantitative traits[END_REF]

  Downloaded from https://academic.oup.com/jhered/advance-article-abstract/doi/10.1093/jhered/esaa002/5728663 by guest on 24 February 2020 Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002 A c c e p t e d M a n u s c r i p t 4

  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)[START_REF] Kaminski | Single nucleotide polymorphisms between two lines of European bison (Bison bonasus) detected by the use of Illumina Bovine 50 K BeadChip[END_REF].

  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 Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002 A c c e p t e d M a n u s c r i p t 6 only considered SNPs mapping to the 29 bovine autosomes on the Bovine UMD3.1 reference genome assembly

  to characterize individual inbreeding at both a global (genome-wide) and a local (locusspecific) 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 Rk that determines the expected length (equal to 1/Rk 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*Rk generations in the past to the individual level of inbreeding, with small (large) values of Rk corresponding to recent (old) ancestors and long (short) HBD segments. Here we considered for both lines a model with five HBD classes (with Rk equal to 4, 8, 16, 32 and 64) and one non-HBD class whose rate was equal to 64 Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002 A c c e p t e d M a n u s c r i p t 8

  , S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002 A c c e p t e d M a n u s c r i p t 10 Results Structuring of the European bison genetic diversity. The plot of the first MDS coordinates of the 183European 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[START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF] 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 (𝜏̂𝐿= 0.139) compared to the Lowland-Caucasian line (𝜏̂𝐿 𝐶 = 0.0250). Assuming a generation time of six years[START_REF] Gautier | Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences[END_REF] 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. Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002 A c c e p t e d M a n u s c r i p t 11

  , 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 Rk. In the Lowland line, HBD segments were concentrated in classes with Rk equal to 8 and 16 whereas in the Lowland-Caucasian line, HBD classes with Rk 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. Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002 A c c e p t e d M a n u s c r i p t 12 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 Rk 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 Rk = 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).

  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 -2 for class with Rk = 8) of the genome associated with recent HBD classes (with a rate < 16) but higher proportions (p = 7.6 x 10 -2 for class with Rk = 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.

  ). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002 A c c e p t e d M a n u s c r i p t 16

  [START_REF] Kardos | Genomic consequences of intensive inbreeding in an isolated wolf population[END_REF] 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[START_REF] Gómez-Sánchez | On the path to extinction: Inbreeding and admixture in a declining grey wolf population[END_REF]. 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 Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002A c c e p t e d M a n u s c r i p t 18 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.

  Downloaded from https://academic.oup.com/jhered/advance-article-abstract/doi/10.1093/jhered/esaa002/5728663 by guest on 24 February 2020 Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002 A c c e p t e d M a n u s c r i p t 19 various age-based related classes of HBD segments may represent a valuable strategy even when samples of different ages are not available.

AFigureFigure 1 .

 1 Figure

Figure 2 .

 2 Figure 2. Characterization of individual inbreeding levels. A) Distributions of the inbreeding

Figure 3 .

 3 Figure 3. Average heterozygosity estimated around 22,602 genotyped SNPs as a function of

Figure 4 .AFigure 1 AFigure 2 AFigure 3 AFigure 4

 41234 Figure 4. Changes of the inbreeding coefficients (upper panel) and of the proportions of the genome

  Downloaded from https://academic.oup.com/jhered/advance-article-abstract/doi/10.1093/jhered/esaa002/5728663 by guest on 24 February 2020

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 Comment citer ce document : Druet, T. (Auteur de correspondance), Olénski , K., Flori, L., Bertrand, A. R., Olech , W., Tokarska, M., Kaminski, S., Gautier, M. (2020). Genomic footprints of recovery in the European bison. Journal of Heredity (février), 29 p. , DOI : 10.1093/jhered/esaa002 A c c e p t e d M a n u s c r i p t 5 Markov model (HMM), named ZooRoH

  monomorphic, ii) displayed a call-rate <0.95, and iii) deviating from Hardy-Weinberg equilibrium (p ≤

	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).

Downloaded from https://academic.oup.com/jhered/advance-article-abstract/doi/10.1093/jhered/esaa002/5728663 by guest on 24 February 2020

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

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.

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).