The type of bottleneck matters: Insights into the deleterious variation landscape of small managed populations

Abstract Predictions about the consequences of a small population size on genetic and deleterious variation are fundamental to population genetics. As small populations are more affected by genetic drift, purifying selection acting against deleterious alleles is predicted to be less efficient, therefore increasing the risk of inbreeding depression. However, the extent to which small populations are subjected to genetic drift depends on the nature and time frame in which the bottleneck occurs. Domesticated species are an excellent model to investigate the consequences of population bottlenecks on genetic and deleterious variation in small populations. This is because their history is dominated by known bottlenecks associated with domestication, breed formation and intense selective breeding. Here, we use whole‐genome sequencing data from 97 chickens representing 39 traditional fancy breeds to directly examine the consequences of two types of bottlenecks for deleterious variation: the severe domestication bottleneck and the recent population decline accompanying breed formation. We find that recently bottlenecked populations have a higher proportion of deleterious variants relative to populations that have been kept at small population sizes since domestication. We also observe that long tracts of homozygous genotypes (runs of homozygosity) are proportionally more enriched in deleterious variants than the rest of the genome. This enrichment is particularly evident in recently bottlenecked populations, suggesting that homozygosity of these variants is likely to occur due to genetic drift and recent inbreeding. Our results indicate that the timing and nature of population bottlenecks can substantially shape the deleterious variation landscape in small populations.


| INTRODUC TI ON
Deleterious mutations are expected to be held at low frequency by an equilibrium between the rate at which they arise by mutation and the efficacy of purifying selection at removing them from the population (mutation-selection balance) (Lynch, 2010;Ohta, 1992).
However, the number and frequency of deleterious genetic variants segregating in a population are affected by many evolutionary forces, including artificial selection, genetic drift and genetic hitchhiking, which is the change in allele frequency of a variant that is passed along together with another variant under positive selection (Charlesworth, 2006;Smith & Haigh, 1974).
In small populations, the mutation-selection balance is challenged by population contractions, which reduce the efficacy of purifying selection to remove harmful mutations (Ohta, 1973). As a result, the genetic load, defined as the reduction in mean fitness in a population caused by deleterious variation relative to a mutation-free population (Kimura, Maruyama, & Crow, 1963), is predicted to be larger. In the long-term, the high genetic load and the rapid increase in frequency of harmful mutations could impact population survival and genetic diversity, increasing the risk of inbreeding depression (Kimura et al., 1963).
Genetic drift, or the random fluctuations in the number and frequency of alleles, is mostly responsible for the deleterious genetic landscape in small populations. However, as studies in plant and animal species have suggested, the extent to which small populations are subjected to genetic drift considerably varies depending on the nature and time frame in which the bottleneck occurs (Liu, Zhou, Morrell, & Gaut, 2017;Mardsen et al., 2016;Zhang, Zhou, Bawa, Suren, & Holliday, 2016). For instance, a long-term population decline is expected to result in a lower proportion of amino acid changing variants, along with a reduction in the additive genetic load, due to purifying selection acting against deleterious variants (Mardsen et al., 2016). However, if populations have undergone recent and sudden declines, deleterious variation is predicted to be mainly shaped by genetic drift (Ohta, 1973).
Domesticated species are an excellent model to investigate the consequences of population bottlenecks on genetic and deleterious variation. This is because their demographic history is characterized by multiple population contractions associated with domestication, breed formation and intense selective breeding (Bosse, Megens, Derks, Cara, & Groenen, 2019;Makino et al., 2018;Mardsen et al., 2016;Moyers, Morrell, & McKay, 2017). Domestication involves the (partial or complete) isolation of a number of individuals from a wild progenitor population and entails drastic changes in the nature and strength of selective forces acting on the population, as well as its size (Larson & Fuller, 2014). The domestication bottleneck is usually followed by a long period of relatively weak and varying artificial selection, during which the reduced N e may either be stable or fluctuate depending on human-driven selection. Contrary to the long-term domestication process, breed formation is a more recent event that often entails intense selection over short time periods and is coupled with limited recombination and an additional reduction in N e (Moyers et al., 2017). In this study, traditional fancy breeds of chicken were used as a model species to investigate the consequences of two types of bottlenecks on deleterious variation: the severe domestication bottleneck occurred some thousands of years ago and the recent population decline accompanying breed formation in the last decades.
Since their development in the 16th and 18th century (Dana et al., 2011), traditional fancy breeds have persisted at small population sizes and comprised normal-sized (large fowl) and miniature (bantam) breeds. These traditional breeds experienced domestication only, which was based upon preferential breeding of birds exhibiting specific morphological features. The subsequent long period of weak and varying artificial selection resulted in the foundation of numerous breeds that are nowadays identified by an accurate phenotypic description (Tixier-Boichard, Bed'hom, & Rognon, 2011). In the last decades, hobby breeders have become interested in miniature forms of historical large breeds, which are called neo-bantams, and were initially created by crossing a large fowl with a bantam individual. Even though mating between neo-bantams has recently started to become very popular among hobby breeders, the selection purpose of obtaining an individual exhibiting all of the standard large fowl characteristics still remains (Bortoluzzi et al., 2018). The recent creation of neo-bantam breeds involved, on top of domestication, an additional population bottleneck. As we showed in our previous study, the reduced N e and parent-offspring mating pursued within a neo-bantam breed to consolidate favourable traits considerably increased the level of inbreeding (Bortoluzzi et al., 2018). Although we expect the recent and sudden bottleneck to have acted differently on the accumulation of deleterious variants relative to the domestication bottleneck experienced by historical breeds, its effect on genome-wide patterns of deleterious variation remains unclear.
Accurate predictions of deleterious variants are essential when assessing their contribution to phenotypic variation (Kono et al., 2016). To date, numerous approaches have been developed and applied to nonhuman species (Liu et al., 2017;Makino et al., 2018;Mardsen et al., 2016;Renaut & Rieseberg, 2015;Robinson et al., 2016;Zhang et al., 2016), of which the Sort Intolerat From Tolerant (SIFT) approach is among the most widely used. However, as shown in Kono et al., (2016) and Derks et al. (2018), additional filtering criteria should be applied to the set of deleterious mutations to improve the reliability of the prediction. These criteria should include orthologous genes to minimize the effect of off-site mapping of sequence reads, RNA expression of protein-coding variants and the use of different prediction approaches (Derks et al., 2018;Kono et al., 2016).
We here expanded the approach of Derks et al. (2018) to predict deleterious mutations in domestic chickens by addressing a potential source of bias not previously investigated. That is reference bias, which is the higher probability of calling a variant as reference. We corrected for that by polarizing all protein-coding variants by ancestral and derived state, rather than reference and nonreference, to not underestimate the inferred number of nonsynonymous and deleterious variants.
Whole-genome sequencing data from 97 chickens representing 39 traditional fancy breeds were here used to directly examine the impact of different population bottlenecks on patterns of deleterious variation in small populations. Overall, we find that the recent population bottleneck associated with the creation of neo-bantams has resulted in a higher proportion of deleterious variants relative to large fowl and bantam counterparts, as genetic drift has reduced the efficacy of purifying selection to eliminate harmful mutations. We also observe that most deleterious variants are found in long tracts of homozygous genotypes, suggesting that homozygosity of these variants is likely to occur due to genetic drift and recent inbreeding.
Our results indicate that the time frame and nature of the bottleneck can substantially shape the deleterious variation landscape in small populations.

| Samples and sequencing
DNA of 97 individuals from 39 traditional chicken breeds from the Netherlands was used for whole-genome sequencing on an Illumina HiSeq 3,000. Four samples from each of the known living Gallus species were also sequenced for this study (Table S1).
Based on their demographic and selection history, samples were

| Principal component analysis
Genetic relationships among the 95 sequenced individuals (Gallus species excluded) were investigated through a principal component analysis (PCA) performed on the filtered vcf files in PLINK v.1.9 (Purcell et al., 2007). First and second principal components were plotted using R v.3.2.0 (R Core Team, 2013).

| Heterozygosity analysis
Individual heterozygosity was estimated for each of the 95 chickens on a genome-wide scale by dividing the genome into nonoverlapping windows of 10 kb (Table S2 in Appendix S1). Within each window, heterozygous variants were called only if their depth of coverage met the minimum and maximum threshold, which were set at four and two times the average coverage, respectively ( Figure S2). The total number of heterozygous sites called in a 10-kb window was corrected for the total number of sites not called because of low coverage, following (Bosse et al., 2012).
Insufficiently covered bins (10-kb windows with less than 1,000 well-covered sites) were excluded from the genome-wide autosomal heterozygosity analysis.

| Runs of homozygosity (ROHs)
Runs of homozygosity were extracted from the genome of the 95 sequenced individuals implementing the method developed by (Bosse et al., 2012). Information on heterozygosity was used to identify autosomal ROHs, which are here defined as genomic regions showing lower heterozygosity than expected based on the genome-wide average. To identify ROHs, we considered 10 consecutive bins at a time (100,000 bps) in which we calculated the average window heterozygosity that was then compared to the average genome-wide heterozygosity. The 10 consecutive bins were retained as candidate ROHs only if their level of heterozygosity was below 0.25 the average genomic diversity. All 10 consecutive bins that did not meet these criteria were considered to contribute to the genome-wide heterozygosity level outside ROHs. Local assembly or alignment errors were avoided as much as possible by relaxing the threshold within candidate homozygous stretches allowing for maximally twice the average heterozygosity in a bin, only if the heterozygosity within the candidate ROH did not exceed 1/3 the average genomic diversity (for more information refer to [Bosse et al., 2012]). Insufficiently covered bins were not considered in the actual size of each ROH but were considered in the calculation of the actual ROH length (assuming that all bins were highly covered). ROHs with insufficient coverage (less than 2/3) were removed from our calculations. From the final list of individual ROHs (Table S2 in Appendix S1), we classified ROHs into three size classes, each of them corresponding to a specific demographic event, including past relatedness (short ROHs: <100 kb), background relatedness (medium: 0.1-3 Mb) and recent relatedness (long: ≥3 Mb).

| Population history estimation
SMC++ was used on unphased whole-genome sequencing data to estimate population history (Terhorst, Kamm, & Song, 2017). Only samples with an average genome coverage >10× and percentage of missing sites <10% were considered (Table S5). Population history was estimated for each breed separately setting the mutation rate to 1.9 × 10 -9 site -1 year -1 (Nam et al., 2010) and generation time at 1 year.

| Genetic load
Genetic load was calculated as the ratio of nonsynonymous deleterious to synonymous sites in each individual and averaged across individuals within each of the three management groups. Genetic load was separately estimated for heterozygous and homozygous derived alleles.

| Site-frequency spectrum (SFS)
The derived allele frequency (DAF) spectrum was calculated for synonymous, tolerated and deleterious variants, considering only bi-allelic SNPs. We then generated a histogram with 10 bins (with steps of 0.1 allele frequency) starting from a very low (0-0.10) to a very high (0.90-1.0) derived allele frequency.

| Enrichment of ROHs for deleterious variants
The distribution of putatively deleterious mutations inside and outside of ROHs was investigated following the method proposed by Szpiech et al. (2013). Homozygous derived variants were grouped into nondamaging or putatively neutral (e.g. synonymous and tolerated) and damaging (e.g. deleterious and LoF). The occurrence of damaging and nondamaging variants was investigated inside and outside each ROH size class. Coordinates of ROHs were used to calculate the fraction of the genome covered by any ROH and by each ROH size class as: where L g is the total length of the genome, L ROH is the total length of ROHs, i is the individual, and j is the ROH class jò S,M,L,A representing small, medium, long and any ROH, respectively (Table S3 in Appendix S1).

| RE SULTS
Patterns of deleterious variation were investigated using whole-genome sequencing (WGS) of 39 traditional chicken breeds (Table S1).
On average, 13.4× coverage was generated for each individual (Table   S1 in Appendix S1). The population-based variant calling approach identified 17 million SNPs and 1.2 million insertions/deletions (indels) (Table S2). Variants were distributed with an average density of 20 SNPs/100-kb, ranging from 0 to 82. The average transition to transversion (Ts/Tv) ratio was 2.58 (Table S2), which is in line with previous findings in commercial chicken populations (Derks et al., 2018). Samples origin was validated with the principal component analysis ( Figure S1).

| Population history is responsible for the current autozygosity landscape
Genome-wide autosomal heterozygosity and ROHs were used to investigate the extent and nature of genetic variation in our populations. Whole-genome heterozygosity ranged from 12.2 to 40.8 SNPs/10-kb (Table S2 in Appendix S1). On average, neo-bantams showed slightly lower heterozygosity than their original large fowl counterparts, though the level of heterozygosity was considerably higher than that observed in the bantam breeds ( Figure 2a).
The level of heterozygosity increased almost two-fold in all breeds when excluding ROHs, with neo-bantams showing slightly higher heterozygosity than both source populations (large fowls and bantams) ( Table S2 in Appendix S1; Figure 2b). The lower genome-wide heterozygosity observed in neo-bantam and bantam breeds is therefore explained by their higher average ROH size (Figure 2d).
On average, 25% of the genome in neo-bantams was covered by long ROHs, 0.6% by short, and 17% by medium ROHs (Table S3 in Appendix S1). Of the bantam breeds, the Eikenburger bantam was the most inbred, with up to 70% of the genome covered by ROHs (Table S3 in Appendix S1).
To further investigate how historical demographic changes have shaped the genomic patterns of homozygosity observed in our populations, we decided to infer past effective population size (N e ). According to our results, the chicken ancestral population size remained stable up to approximately 10,000 years, after which it dropped from an initial N e of 10 6 to 10 4 -10 3 ( Figure S3). Both management groups showed a constant flattening population size which has hardly recovered since the bottleneck.

| More rare than fixed deleterious variants
The role of genetic drift and purifying selection was investigated by annotating variants with respect to their effects on the amino acid sequence (Table S3). We also annotated alleles as ancestral and derived using three wild Gallus species as an outgroup. After filtering for RNAseq coverage and 1:1 orthologues, the final set of variants comprised 61,567 synonymous, 16,840 nonsynonymous tolerated, 3,833 nonsynonymous deleterious and 755 loss of function (LoF) mutations. Of the initial set of deleterious variants, 1,674 were classified as deleterious by both SIFT and GERP++. The efficacy of selection at removing deleterious variants from a population was investigated by looking at the distribution of the derived allele frequency (DAF) spectrum. The frequency spectrum showed more rare (DAF <0.1) derived alleles than nearly fixed or fixed deleterious alleles (DAF ≥0.9; Figure 4). We observed similar DAF spectra for large fowls ( Figure S4a) and neo-bantams ( Figure S4b), with neo-bantams showing slightly higher derived allele frequency than large fowl counterparts, even up to a DAF of 0.5.
Moreover, neo-bantams showed fewer nearly fixed or fixed deleterious variants than large fowl counterparts.

| The effects of genetic drift on deleterious variation
In traditional fancy breeds, the total number of derived alleles (heterozygous and homozygous derived) was lower for deleterious and LoF variants relative to putatively neutral ones (synonymous and nonsynonymous tolerated) ( Figure S5). Compared with large fowls, neo-bantams were slightly more enriched in the total number of deleterious and LoF homozygous derived mutations ( Figure   S5d). Despite these differences in the total number of homozygous derived genotypes, we decided to perform a likelihood ratio test (LRT) to formally test for individual differences between neobantams and large fowls. According to the likelihood ratio test, the number of homozygous derived genotypes was not significantly different between large fowl and neo-bantam counterparts for deleterious (p-value: 0.730) and LoF (p-value: 0.272) variants (Table 1). Significant were the differences for synonymous (pvalue: 3.442e-08) and nonsynonymous tolerated (p-value: 0.015) variants. We also investigated in each of the three management groups the total genetic burden resulting from the accumulation of deleterious mutations (genetic load; Figure 5). The deleterious to synonymous ratio of heterozygous variants was, on average, lower in large fowls compared with bantam and neo-bantams.
The same ratio when considering homozygous derived variants was, on average, slightly higher in large fowls than neo-bantams, which, on the other hand, showed extensive variation ( Figure 5).
Contrary, bantam breeds showed little variation with an average higher deleterious to synonymous ratio than that of large fowls and neo-bantams.

| Deleterious variation and demographic history
The study of ROHs offers a new basis for assessing the mechanisms by which demography and selection produce patterns of deleterious variation. The total number of putatively neutral homozygous derived sites was higher than that of damaging sites ( Figure S6). 3.524e-08; Figure S6a). The fraction of damaging and nondamaging homozygous derived genotypes in ROHs positively correlates with the total genomic ROH coverage (Pearson correlation: 0.956, p-value:

Moreover
<2.2e-16 for damaging; Pearson correlation: 0.981, p-value: <2.2e-16 for nondamaging; Figure 6). We also observed that each ROH size class (short, medium, long and any) is more enriched for deleterious F I G U R E 4 Derived allele frequency spectra of the 39 traditional chicken breeds. Derived allele frequency was inferred for synonymous, nonsynonymous tolerated (SIFT score ≥0.05) and nonsynonymous deleterious (SIFT score <0.05) mutations  The null model states that the proportion of homozygous derived genotypes that a large fowl carries is equal to that of the same genotypes carried by a neo-bantam, so that p lf = p nb. b Contrary to the null model, in the alternative model the proportion of homozygous derived genotypes carried by a large fowl individual is different from that of a neo-bantam (p lf ≠ p nb ). c Likelihood ratio test for differences in the number of homozygous derived genotypes per individual between large fowl and neo-bantams. The chisquare distribution with one degree of freedom of Λ was used to calculate p-values.
homozygous derived variants than for nondamaging homozygotes.

| D ISCUSS I ON
In this study, we used whole-genome sequencing data from traditional fancy breeds of chicken to investigate the consequences of population bottlenecks on genome-wide patterns of deleterious variation in small populations. To do so, we combined individuals from multiple breeds with similar demographic history into one population to better estimate genetic and deleterious variation.
Such approach was most suitable given the low number of individuals per breed (between 1 and 4), which is the direct consequence of the threatened population size of most of these breeds (http:// edepot.wur.nl/424249). Even though (slightly) different breeds were grouped into the same population, we expect potential bias to be minimal, as already shown in Bortoluzzi et al. (2018) when using genome-wide SNP chip data.
In line with the small-population paradigm, we showed that the size of a population (N e ) is an important evolutionary factor in determining the level of genetic variability and the effectiveness of purifying selection at removing harmful mutations (Caughley, 1994;Charlesworth, 2009). In fact, as we showed, populations of small N e have a lower genetic diversity and high level of inbreeding, along with being more affected by genetic drift.
In traditional large fowls, the accumulation of deleterious alleles is characteristic of a population that since the domestication bottleneck has persisted for a long period of time at small size.
Therefore, in these populations deleterious mutations of especially small effects are expected to have risen in frequency and become fixed (Hedrick, 2001). As several studies have shown, domestication substantially decreases the effective population size and efficacy of purifying selection, which in turn reduces the genetic diversity and increases the mutational load (Moyers et al., 2017). These major effects have been observed in many species despite the multiple domestication centres and large population size of the ancestor (in the case of chicken, the red jungle fowl) (Kanginakudru, Metta, Jakati, & Nagaraju, 2008;Makino et al., 2018;Miao et al., 2013). For example, (Mardsen et al., 2016) recently observed that the dog genome harbours more amino acid changing variants than that of the wild wolf ancestor. The higher proportion of deleterious variants has also led to an increase in the additive genetic load in many dog breeds, which clearly indicates that the efficacy of purifying selection is lowered by strong population contractions accompanying domestication (Cruz, Vila, & Webster, 2008;Mardsen et al., 2016). Similar conclusions have been reached in other domestic animal species (Makino et al., 2018;Schubert et al., 2014) and plants (Liu et al., 2017;Lu et al., 2006). Even though demographic contractions associated with domestication have a major impact on the genome-wide genetic and deleterious variation, processes that co-occurred during domestication have recently questioned the role of domestication itself in increasing the mutation load. For example, the shift in mating system from outcrossing to predominantly selfing rice has been suggested to have substantially influenced the occurrence of deleterious mutation in domesticated rice (Liu et al., 2017).
It is, however, not to exclude that also the long period of weak and varying artificial selection for desirable traits following domestication could have further reduced N e (Moyers et al., 2017). is not large enough to discount the effects of genetic drift. As a result, weakly deleterious mutations accumulate in the genome due to genetic drift, as purifying selection does not have the time to purge these harmful mutations (Luikart, Allendorf, Cornuet, & Sherwin, 1998). The central role of genetic drift observed in our recently bottlenecked populations on deleterious variation is also observed in small populations under natural selection (Abascal et al., 2016;Hedrick, Kardos, Peterson, & Vucetich, 2016;Pollinger et al., 2010;Robinson et al., 2016). In their study, Robinson et al. (2016) showed that a severe bottleneck occurred ~30 generations ago has substantially reduced the already small effective size of the island fox from originally 64 individuals to fewer than a dozen. This recent population contraction has severely affected not only the genetic variation of the species, but also the genetic load. In a recent study on commercial chicken lines, Derks et al. mating between family members is intentionally pursued to select for specific traits, the proportion of homozygous segments in individual genomes is expected to substantially increase along with that of slightly deleterious mutations. Therefore, we expect the viability of the traditional breeds investigated in this study to strongly depend on future breeding preferences, which, if not genetically managed, are likely to limit the full exploitation of their genetic potential.

| CON CLUS IONS
In this study, we showed that the timing and nature of a population bottleneck can substantially shape the deleterious variation landscape in small populations. In particular, we showed that populations kept at small size for long period of time since the bottleneck have a reduced burden of deleterious alleles compared with recently bottlenecked populations. The reduced deleterious burden in these populations, which is also linked to a reduced number and total length of ROHs across the genome, is likely responsible for their genetic success.
According to our study, facilitating purging of deleterious mutations through inbreeding avoidance should be at the core of future breeding and conservation programmes in small populations (de Cara, Villanueva, Toro, & Fernandez, 2013). However, genomic information on deleterious variation can, and should, be incorporated and used in the development of conservation programmes that assure the long-term survival and enhance the genetic diversity of small populations. Fitness-related traits should also be considered to better measure potential fitness consequences at the individual and population level associated with recessive deleterious mutations.

ACK N OWLED G EM ENTS
We would like to acknowledge the poultry breeding clubs and hobby breeders for their collaboration during sampling. A special thanks to Miguel Correa-Marrero and Carlos de Lannoy for their useful suggestions and input on improving the runs of homozygosity pipeline used in the study. The research leading to some of these results has been conducted as part of the IMAGE project, which received funding from the European Union's Horizon 2020 Research and Innovation Programme under the Grant Agreement No. 677353. The research was also partly funded by the STW-Breed4Food Partnership, project number 14283: From sequence to phenotype: detecting deleterious variation by prediction of functionality.

DATA AVA I L A B I L I T Y S TAT E M E N T
Whole-genome sequencing data from this study have been submit-