Deleterious alleles in the context of domestication, inbreeding, and selection

1.1. Inbreeding depression

Small populations are more likely to suffer from inbreeding, which can negatively impact health and reproduction (Lynch et al., 1995). The decline in fitness observed in inbred progeny, relative to outbred progeny, is known as “inbreeding depression” (Keller & Waller, 2002). Inbreeding depression has largely been attributed to the accumulation of recessive harmful mutations in the genome: Inbreeding increases the probability of these mutations to become homozygous (i.e., Agrawal & Whitlock, 2012; Charlesworth & Willis, 2009; Ohta, 1973). Apart from increased homozygosity of recessives, the overdominance hypothesis assumes a heterozygote advantage, and therefore, overall loss of diversity due to inbreeding reduces the advantage of the heterozygotes (Crow, 1970a,b). In this paper, we focus on the detrimental effect of mostly recessive deleterious alleles, and their combined effect referred to as the mutational load.

With advancing sequencing technologies, the mutational load in the genome of an individual can be estimated from sequence data with increasing accuracy. We should bear in mind that the number of deleterious mutations does not necessarily need to differ between an inbred and an outbred individual for differences in fitness to exist; the key concept here is that while harmful mutations generally have a small fitness effect in heterozygous state, in homozygous state, they are expressed, causing, for instance, heritable diseases (Figure 1). Since ROH are formed because both haplotypes are identical by descent (IBD), the probability of a recessive deleterious allele with a low‐frequency p to become homozygous is higher in such regions than outside IBD regions (frequency p ²: (Szpiech et al., 2013).

1.2. Identification of deleterious alleles

Alleles with a putative effect on the phenotype, both beneficial and deleterious, are thought to be younger than variants with no effect with the same allele frequency (Maruyama, 1974). Understanding the factors that cause harmful mutations to increase in frequency in a genome will facilitate prediction of genetic load in current populations, but also aid to avoid a high genetic load in the future (Garcia‐Dorado, 2012). A popular approach to identify harmful mutations that are (almost) lethal is to screen populations for the absence of specific variants in homozygous state. If populations contain heterozygous carriers of such variants, homozygotes would be expected based on allele frequencies and carrier–carrier matings (Derks et al., 2017; Pausch et al., 2015; VanRaden, Olson, Null, & Hutchison, 2011). This method relies on the biological implications of the lethality of a variant—resulting in an absence of the allele in homozygous state. Depending on the frequency of the lethal allele in the population, large sample size is often required to identify such alleles (Derks et al., 2017). Therefore, most of the studies implementing this depletion of homozygotes approach use large panels of genotyped animals. An alternative approach is to predict deleteriousness from the functionality of a mutation. Next‐generation sequencing has opened up exciting possibilities to actually pinpoint potentially harmful mutations in individual genomes (Henn, Botigue, Bustamante, Clark, & Gravel, 2015; Li et al., 2010). The deleteriousness of a variant can be predicted based on its effect on gene functioning (such as protein changing, stop‐gain, stop‐lost), for example by assessing the degree of conservation of an amino acid residue across species (Cooper et al., 2005; Kumar, Henikoff, & Ng, 2009; Wang, Li, & Hakonarson, 2010). The implementation of multiple genome annotations, such as specific gene function or regulatory elements, has proven successful for in silico predictions of the effect of disease causing mutations in humans (Kircher, Witten, Jain, B. J. O'Roak, & Shendure, 2014). These techniques have recently been applied to estimate genetic load in human populations and domesticated species (Charlier et al., 2016). However, we need to remain cautious with such assessments of deleteriousness, since they rely on predictions that mostly have not been validated with experiments. High‐impact variation may very well not be lethal in a domestic setting and could even be perceived beneficial if they serve a particular breeding goal. Nevertheless, such advances in predicting the effect of variants have high potential for bridging the gap between sequence information and fitness effects.

1.3. Application to domestic species

A combination of both approaches to detect deleterious variants has proven successful in cattle (Pausch et al., 2015). Recent findings based on predicted deleterious variants from sequence data corroborate the classical theory by Mukai et al. (1972) that deleterious alleles are generally at low frequency (e.g., Mezmouk & Ross‐Ibarra, 2014). Multiple factors contribute to inbreeding depression, and epistatic effects across loci should be considered as well. Nevertheless, quantifying harmful mutations in single genomes is an important first step toward the “genomic characterization of genetic load.” An increase in studies assess mutational load in sequence data from domestic species. Some general trends about the effects of domestication on the burden of harmful mutations are now emerging (Makino et al., 2018; Zhou, Massonnet, Sanjak, Cantu, & Gaut, 2017). Here, we review recent studies on mutational load in domestic species and use re‐sequence data from pig and chicken to demonstrate how specific hypotheses about the burden of deleterious variants can be tested. Specifically, we will discuss two major drivers of mutational load in genomes: bottlenecks and selection.

2.1. Domestication increases load

Several studies report the effect of domestication on mutational load. In line with predictions, the bottlenecks associated with domestication have not only reduced genetic diversity of domestic species, but also increased the mutational load in them as well. Domestication bottlenecks have indeed been suggested to have substantially increased the mutational load in dogs, referred to as the “legacy of domestication” (Cruz, Vila, & Webster, 2008). In dog genomes, the ratio of amino acid changing heterozygosity to silent heterozygosity (variants presumed to have no effect) was higher than in their wild ancestors, gray wolves (Marsden et al., 2016). These results indicate that the ability of purifying selection to remove weakly deleterious variants is lowered by the bottlenecks. A similar phenomenon has been reported for domestic horse genomes, where an excess of deleterious mutations is thought to be caused by domestication and inbreeding (Schubert et al., 2014). Also in crops, deleterious mutations seem to have accumulated in domestic lineages (Kono et al., 2016; Lu et al., 2006).

2.2. Type of bottleneck is important

The time frame in which the bottleneck occurs has a strong effect on the potential of natural selection to eliminate harmful mutations. If the bottleneck is severe and sudden, (local) genomic recombination is scarce, then drift can exert a maximal effect. If, however, the same population would decline over a long period of time, purging might enable deleterious variants to be removed from the population, resulting in a lower load (Charlier et al., 2016; Hedrick & Garcia‐Dorado, 2016). This begs the question of whether it was domestication itself that increased mutational load, or rather one of the processes that co‐occurred with it. The strong reduction in effective population size, going from wild to domesticated, could have driven the deleterious alleles to higher frequency (Liu, Zhou, Morrell, & Gaut, 2017). Artificial selection could have reduced the effective population size even further.

Most domestic animal species are thought to have experienced a relatively strong population bottleneck, although the picture can become complex for species that went through multiple domestications in different regions. In livestock, domestication is no longer seen as a single, discrete event. Rather, substantial and continuous gene flow from wild populations has occurred during the process of domestication (Frantz et al., 2015; Scheu et al., 2015). In rice, Liu and colleagues suggest that it was not domestication itself, but the shift in mating system from outcrossing to predominantly selfing that had a substantial influence on mutational load (Liu et al., 2017). This suggests that the process of domestication as well as the management regime under which current lines have been formed could have influenced the occurrence of deleterious mutations in domestics. The intensity of the domestication bottleneck is thought to have influenced the difference in mutational load between annual and perennial crops (Zhou et al., 2017). A recent analysis on multiple domesticated species concludes that in domestic plant and animal genomes, an elevated proportion of deleterious genetic variation is present, with European pigs as an exception (Makino et al., 2018). This creates opportunity to investigate further how domestication has elevated mutational load.

2.3. Case study in pig and chicken: application to sequence data

We used genotype and re‐sequence data from pigs and chicken to investigate the mutational load in wild and domestic populations, as described in supplementary materials. These species represent highly different domestication histories and selection regimes. We compared the ratio of predicted deleterious heterozygosity with silent heterozygosity within individual genomes to estimate mutational load in domestic and (semi) wild pigs and chickens, using the method described by Renaut and Rieseberg (2015), applied to sunflowers.

Pigs were domesticated at least twice, independently, giving rise to the current Asian and European‐based domestic clades (Kijas & Andersson, 2001; Larson et al., 2005). European and Asian wild boar form an excellent model for their wild ancestors, enabling direct comparisons between the wild and domestic form (Bosse, Megens, Frantz, et al., 2014; Bosse, Megens, Madsen, et al., 2014). The use of pigs from two different geographic regions (Asia and Europe) and subjected to different domestication events (wild and domestic, including local and commercial populations) enables the study of the impact of demography and domestication on the distribution of deleterious mutations. In chicken, however, the wide variety of domestic breeds are thought to stem from red jungle fowl (Fumihito et al., 1996) and involve more complex demographic history, including multiple regional centers of domestication across Asia (Miao et al., 2013).

2.4. Consequences of increased load

Whether the high estimated mutational load in domestic genomes will lead to lowered viability is questionable. Many commercial lines are bred especially for reproduction‐related traits or through the production of hybrids, which will dilute the effect of mutational load. Current ongoing work has shown that life history traits such as lifespan or propagule size are a key to understanding levels of genetic diversity and mutational load (Romiguier et al., 2014). So far, the results on mutational load estimated from the ratio of synonymous to nonsynonymous mutations are fairly evident in indicating that purifying selection is higher in the wild (Chen, Glemin, & Lascoux, 2017). The sources of genetic variation that create genetic diversity where selection can act upon are the same for both domestic and wild animals. However, in most circumstances, domestic animals have been so strongly selected for the sake of production, that their genetic diversity is extremely low. Moreover, artificial selection for specific traits sometimes results in so much inbreeding that it is threatening the survival of the breed. Clear examples can be seen in dogs such as the Norwegian Lundehund suffering from intestinal problems, but cancer, eye, and heart diseases are also common (Kettunen, Daverdin, Helfjord, & Berg, 2017; Schoenebeck & Ostrander, 2014). Combined, these effects of inbreeding in domestics have led to animals and plants that may have lost their ability to face environmental challenges, should they come from climate change or new diseases. Finally, the circumstances in which domestic animals are kept do not reflect the (harsh) environment they face in the wild. The loss of adaptive potential, disease resilience, and fitness reduction is of major concern for monocultures such as bananas and in salmon (Araki, Berejikian, Ford, & Blouin, 2008; Garcia de Leaniz et al., 2007). Together with the fact that the traits selected in domestic animals are usually deleterious in the wild make the detrimental effects of the potentially harmful mutations hard to observe (Hedrick & Garcia‐Dorado, 2016).

3.1. Increased load due to selection

In humans, deleterious alleles are thought to have increased in frequency due to linkage to sites that have been under positive selection. (Chun & Fay, 2011). Indeed, in domestic species, we find similar evidence of selected regions to be over‐represented with predicted deleterious alleles. In dogs, an increased load was found within regions under selection, inferred from the genome (Marsden et al., 2016). A different approach is to assess the mutational load in genes known to affect phenotypic traits of interest. In plants, such genes were shown to contain proportionally more deleterious variants (Kono et al., 2016; Lu et al., 2006; Mezmouk & Ross‐Ibarra, 2014). The linkage of deleterious variants to genes under balancing selection could also lead to a local overrepresentation of deleterious variants. A recent study in cattle (Kadri et al., 2014) showed that a long known antagonism between fertility and milk production is actually due to such a phenomenon. While a major QTL for fertility with effects in milk production is under balancing selection in nature, the change to directional selection in selection for milk production has led to an extremely reduced fertility in cows. Balancing selection in livestock has probably an underestimated role yet to be explored.

3.2. Case study in pig and chicken: load in runs of homozygosity

Together, inbreeding and strong directional selection can increase the proportion of homozygous segments in individual genomes. Following the rationale of Szpiech et al. (2013) that ROH are enriched for deleterious variants (in homozygous state), we tested this hypothesis in pig and chicken lines. In three white layer lines, the number of predicted deleterious homozygous variants was higher in ROH‐regions compared to the rest of the genome. As a result, ROH contain proportionally more homozygous deleterious alleles than the rest of the genome (Figure 4), in line with findings in humans (Szpiech et al., 2013). Also, in commercial European pigs, ROH‐regions contained proportionally more deleterious homozygous variants. Interestingly, however, European wild boar did not display this pattern (Figure 5c).

3.3. Disentangling selection and drift

A higher proportion of deleterious alleles in ROH were also observed in cattle (Zhang, Guldbrandtsen, Bosse, Lund, & Sahana, 2015). In humans, the longest class of ROH was most enriched for deleterious variants in humans, whereas in cattle, shorter ROH contained proportionally more deleterious alleles (Zhang et al., 2015). A possible explanation for these different patterns is that in humans, inbreeding is mostly responsible for the formation of long ROH and hard selective sweeps are rare (Hernandez et al., 2011) By contrast, in cattle, artificial selection intensity is high and sweeps are abundant (Kim et al., 2013). This indicates that the different mechanisms that result in ROH, namely selection and drift, result in different patterns of detrimental alleles. In dog genomes, when regions under selection were excluded, an increase in deleterious alleles compared to wolves could still be observed, suggesting that drift rather than the hitch‐hiking effect associated with selection is causing the increased load in dogs (Marsden et al., 2016). Also, the role of recombination in generating genomic diversity patterns and ROH should not be neglected; recombination rate locally influences effective population size in a genome, which in turn affects the efficiency of natural selection (Begun & Aquadro, 1992; Bosse et al., 2012). Finally, the number of generations since the last common ancestor and the speed of inbreeding influence the length of ROH stemming from inbreeding, and the associated mutational load (Hedrick & Garcia‐Dorado, 2016). Therefore, minimizing inbreeding in modern breeding practices by specifically avoiding long IBD segments in optimal contributions may avoid an increase in mutational load (de Cara, Villanueva, Toro, & Fernandez, 2013a,b).

Distinguishing regions under selection and IBD segments due to inbreeding may shed new light on the assumed detrimental effects of domestication. Specifically targeting genes known to be associated with favorable traits, or regions of reduced heterozygosity on a population scale, can aid in distinguishing ROHs stemming from inbreeding from ROHs stemming from selection.

4.1. Final remarks

Overall, domestication seems to have elevated load in most domestic species. However, domestication does not lead to genomes with high load per se, indicating that fitness does not need to decline. Accumulation of harmful mutations in the genome is facilitated by population bottlenecks and (artificial) selection, but the details of population demography and selection history are species‐ and breed specific. Some mechanisms seem to apply to multiple domestic species, such as an increase in load due to (domestication) bottleneck. Also, the ancestral load is an important factor influencing load in different breeds, since the source population determines which variants are present in the first place. Artificial selection may further increase load by the facilitation of hitch‐hiking of harmful mutations and further restrictions of effective population size. The rare example of lower load in European domestic pigs compared to European wild boars provides a more positive view on mutational load in domestics. In this case, either the process of domestication itself or the artificial selection and hybridization leading to current breeds may have lowered mutational load. We are only beginning to understand how different selection regimes have influenced the mutational load in domestic lineages, and how estimated mutational load and fitness effects are connected. Consistent or strong trends about the influence of domestication on mutational load are not yet clear. However, some common patterns are emerging. We have entered an exciting era in which population genetics theory can be tested with genomic data to shed new light on the interplay of demography, selection, mutational load and fitness. Integrating knowledge from humans, model species such as Drosophila, domesticated species such as livestock and crops, and endangered wild species in combination with in silico predictions will aid in understanding the genetic load across many different demographic and selection scenarios.

5.1. Sampling

Re‐sequence data from a total of 76 pigs were retrieved from the European Nucleotide Archive (ENA) accession code ERP001813 and included 43 domestic pigs and 33 wild boars from Asia and Europe. Pigs were classified into five groups in accordance with their geographic origin and domestication status: Asian domestic pigs (AD; N = 16), North Asian wild boars (AWN; N = 5), South Asian wild boars (AWS; N = 5); European commercial pigs (ED; N = 27) and European wild boars (EW; N = 23). For chicken, re‐sequence data from eight African village chicken coming from different ecotypes in Kenya were obtained from (Ngeno, 2015). In addition, re‐sequence data from one commercial pig line (N = 9) and three chicken lines (N = 51;66;43) were obtained from (Derks et al., 2017; Derks, Megens et al. (submitted) and Derks et al., 2018) respectively.

¹European Nucleotide Archive (ENA) https://www.ebi.ac.uk/ena/data/view/PRJEB1683; ²(Bosse et al., 2012) ³(Yang et al., 2017) ⁴(Derks et al., 2017) ⁵(Ngeno, 2015) ⁶(Ngeno et al., 2015) ⁷(Derks et al., 2018), ⁸(Derks, Megens et al. (submitted)).

5.2. Alignment and variant calling

All samples were sequenced with Illumina paired‐end short read sequencting technology (Illumina inc). Reads were trimmed to a phred quality >20 and minimum length of both pairs of 40 bp. For pigs, short read alignment was done against the Sus scrofa genome, build 10.2 (Groenen et al., 2012) using BWA‐mem. For chicken, reads were aligned to the chicken reference genome Galgal 4 (International Chicken Genome Sequencing 2004). Variants were called and filtered using GATK (Genome Analysis Toolkit 2.6‐4) on an individual basis to avoid biases resulting from unequal sampling across populations. INDEL realignment, base quality score recalibration, and UnifiedGenotyper module analyses were performed according to GATK Best Practices. These variants were further filtered for read depth between half and twice the average genome‐wide depth of that individual with a minimum depth of 4x. We only considered genic variants that had genotype quality scores >Q20.

5.3. Identification of potentially deleterious variants

To obtain the functional consequences of the variants for pigs, we used Variant Effect Predictor (VEP, (McLaren et al., 2016)) on the ENSEMBL Sus Scrofa 78gene set, 1‐to‐1 orthologues with cattle. For chicken, we used Variant Effect Predictor(VEP) v78 on the ENSEMBL Gallus gallus 78gene set. Variants within a gene with a missense annotation and variants labeled as synonymous variant were retained for downstream analyses. The deleterious effect of the missense mutations were predicted using SIFT (Kumar et al., 2009), and SIFT scores <0.01 were used as deleterious mutations in calculating the ratio deleterious/synonymous.

5.4. Deleterious variants and ROH

Runs of homozygosity (ROH) were estimated from genotype data from the same individuals, obtained from (Yang et al., 2017) and (Derks et al., 2017) for pigs and from (Derks et al., 2018) for chicken. SNP genotyping was performed on the species‐specific Illumina 60K iSelect Beadchip for both species (Ramos et al., 2009). A ROH was defined as a genomic region of at least 1 Mb with at least 20 SNPs supporting the homozygous state, using the –homozyg option in PLINK v.1.9. This length reflects roughly 1 cM and consanguinity 50 generations ago (Howrigan, Simonson, & Keller, 2011), which is sufficient inbreeding for the purpose of this study, yet feasible with the density of our SNP chips. Number of deleterious and tolerated nonreference homozygotes overlapping and nonoverlapping ROH (Del‐ROH and Del‐non‐ROH, respectively) were counted separately in each individual. Differences between the fraction of deleterious and tolerated variants in ROH per individual were tested for significance using the paired t test as implemented in R.