Genomic inbreeding and population structure in rams of Tunisian D’man sheep

Attachments:

Volume 61, Article 03

[Volume 61, Article 03]

456 kB

S. BEN JEMAA¹

S. KDIDI²

M. BOUSSAHA³

E. REBOURS³

M. H. YAHYAOUI²

¹ Laboratory of animal and forage production, National Institute of Agronomic Research of Tunisia, University of carthage, Rue Hédi Karray, Ariana, 2049, Tunisia.

² Livestock and Wildlife Laboratory, Arid Lands Institute, Route Djorf Km 22, 4119 Medenine, Tunisia

³ GABI, INRA, AgroParisTech, University of Paris Saclay, 78350 Jouy-en-Josas, France.

Abstract – D’man sheep is a minority breed in Tunisia that was established nearly two decades ago from a small Moroccan flock. Tunisian D’man is an isolated population which has experienced no gene flow since its introduction. Currently, this population is suffering from an increased incidence of stillbirths and recurrent abortions caused by inbreeding. The availability of dense SNP markers has facilitated the quantification of genomic inbreeding in farm animals. The aim of the present study was to estimate the kinship between the six D’man rams of a Tunisian breeding station and to assess their genomic structure using the Illumina OvineSNP50 BeadChip and comparisons with the French Lacaune and the Tunisian west thin tail breeds. We found moderate to high levels of inbreeding between the old rams ranging from 0.0069 to 0.1202 and low level of genetic diversity (expected heterozygosity ~ 0.26). Substantial level of west thin tail introgression was detected in the young D’man rams although these were supposed to be purebred. Analysis of runs of homozygosity (ROHs) showed that all individuals had at least one ROH > 10 Mb. ROH islands identified within the D’man rams harbored genes whose mutation leads to stillbirth, dystocia, embryo quality reduction and fleece phenotype variation. Our study demonstrates the usefulness of molecular markers in the management of inbreeding in small isolated populations and would be very helpful in the implementation of planned mating scheme based on sire genotypes.

Keywords: D’man sheep, Inbreeding, Single nucleotide polymorphism, Runs of homozygosity.

1. Introduction

The D’man sheep breed originated in the southern part of the Maghreb region (southwestern Algeria and southeastern Morocco). This breed was introduced for the first time in Tunisia in 1994 from a Moroccan founder flock composed of 200 ewes and 25 rams (Rekik et al. 2005). D’man population expanded widely in the country (particularly in the Tunisian oases) since its earliest introduction as it was appreciated by the breeders for its high prolificacy. A recent study showed that the high prolificacy of D’man is mainly due to the segregation of a major prolificacy locus, named FecL within the Tunisian D’man population (Ben Jemaa et al. 2019). This locus was originally shown to have a major effect on the prolificacy of the French Lacaune sheep and to be in complete linkage disequilibrium with an single nucleotide polymorphism (SNP) located in the intron 7 of B4GALNT2 gene (OAR11:36938224T>A) (Drouilhet et al. 2013). The estimated mutated allele effect of the FecL^L locus on litter size was +0.4 to +0.5 lambs per lambing (Martin et al. 2014). Currently, the Tunisian D’man breed can be considered as an island population as it is isolated and it has a small size. Indeed, for sanitary reasons, Tunisian D’man has experienced no gene flow from neighboring countries since its establishment nearly two decades ago thus leading to an increase in the level of inbreeding within this population. Consequently, a frequent onset of lethal genetic disorders in new born lambs was observed within many farms. The inbreeding problem is made worse by the absence of a National pedigree and performance recording system allowing the use of optimal mating strategies. However, the problem is less marked within the state flock of the breeding station of Chenchou which belongs to the Livestock and Pasture Office of Tunisia (OEP) which possesses a pedigree recording system allowing a planned mating scheme based on sire rotation across the seven family lines forming the flock. Chenchou breeding station is the leading supplier of rams for most of the small farmers of Southern Tunisia.

The development of high-throughput SNP genotyping platforms coupled with the gradual reduction of genotyping costs makes it possible to accurately estimate kinship between individuals, thus helping control inbreeding in isolated populations. Furthermore, the availability of dense SNP genotypes allows the identification of potential continuous stretches of homozygous SNP genotypes known as runs of homozygosity (ROH). Characterization of ROHs, particularly their size and frequency, provides information about relatedness within a population (Mastrangelo et al., 2018).

In the present study we aimed to estimate the genomic kinship matrix in the six D’man rams of the Chenchou station and to assess genome-wide inbreeding ROH distribution within these genitors using the Illumina OvineSNP50 BeadChip. Furthermore, we provided, for the first time, a detailed assessment of the genetic structure of the D’man breed using comparisons with the Tunisian West thin tail (WTT) sheep and the French Lacaune (LAC) breed. The choice of these two breeds is prompted by the geographic proximity (for D’man and WTT) and the common presence of the prolific allele FecLL (for D’man and LAC).

2. Materials et Methods

2.1. Blood sampling, DNA extraction and genotyping

Jugular vein blood samples (5 ml per animal) were obtained from all the six rams of Chenchou breeding station under the Tunisian Veterinary Authorities’ rules. Blood samples were collected in EDTA Vacutainer tubes. D’man DNA was extracted using salting out method as described by Miller et al. (1989). DNA samples belonging to the six D’man rams and to 17 WTT individuals (available from Kdidi et al. 2015) were then genotyped on the OvineSNP50 Genotyping BeadChip (Illumina, San Diego, CA, USA) using standard procedures (http://www.illumina.com) resulting in 54242 genotyped SNPs. We also included genotyping data of 25 Lacaune individuals available from the International Sheep Genomics Consortium (Kijas et al. 2012).

2.2. SNP quality control and marker selection

Among the 54242 genotyped SNPs, we used the 49034 ones that were in common between the LAC genotypes on one hand and the D’man and WTT genotypes on the other hand. We then used PLINK ver.1.09 (Purcell et al. 2007) for genotyping data quality control. Samples genotyped for less than 90 % of markers were excluded from the analysis. SNPs genotyped for less than 90 % of the animals and those with MAF less than 0.01 were also discarded. Using these criteria, we excluded one WTT individual due to a low genotyping rate (<90 %). Similarly, 253 SNPs were excluded because of low genotyping rates, and 251 SNPs were also excluded as they were monomorphic (MAF <0.01). An exact test for Hardy-Weinberg Equilibrium (HWE) was then carried within each population separately on the remaining SNPs, which led to 47038 SNPs retained for further analysis.

2.3. Kinship estimation and characterization of runs of homozygosity

Genomic kinship coefficient for each pair of D’man individuals was estimated using Loiselle estimator (LS) implemented in the R package RClone (Bailleul et al. 2016). Loiselle estimator was selected because it uses a correction for small sample sizes, thus making him more accurate than most of the other estimators.

Because kinship coefficients are better estimated if the markers used are not in linkage disequilibrium (LD) with each other, a subset of SNPs in approximate linkage equilibrium was selected by LD-based SNP pruning in PLINK.

ROHs were identified in sliding windows of 20 SNPs using PLINK. No more than five missing calls and one heterozygous SNP were allowed in each window. The minimum length of a ROH segment was set to 1 Mb. A ROH was declared if it contained at least 20 SNPs. The minimum required SNP density was one SNP per 200 kb and the maximum gap allowed between any two consecutive SNPs was 1000 kb. ROHs were classified into six categories (1 to 2, 2 to 3, 3 to 4, 4 to 8, 8 to 10, 10 to 20 Mb).

The --homozyg-group function implemented in PLINK was used to assess ROH islands shared among D’man individuals. In the present study, ROH islands were defined as the homozygous segments shared by at least 50% of individuals. Furthermore, gene content within the identified ROH islands was extracted using BioMart tool on Ensembl. Functional enrichment analysis was then performed on the list of genes overlapping ROH islands using Database for Annotation, Visualization and Integrated Discovery (DAVID) software by comparison to the sheep genome background supplied by DAVID.

2.4. Population structure and genetic diversity

We provided a fine-scale assessment of the genetic structure of the Dman, West thin tail and Lacaune individuals using two methods. First, we used Discriminant Analysis of Principal Components (DAPC) implemented in the R package adegenet (Jombart, 2008) without providing any a priori assignment of the individuals to their population of origin. In this analysis, the optimal number of genetic clusters that best describes the data was identified by running a clustering algorithm called k-means and comparing the different clustering solutions using the Bayesian Information Criterion (BIC). Individuals were then assigned to these clusters.

Second, we used ADMIXTURE 1.23 software (Alexander et al. 2009) to estimate proportions of ancestry from each contributing genetic cluster in the three populations. DISTRUCT software (Rosenberg, 2004) was then used to graphically display ancestry within each individual.

Finally, we used expected (Hs) and observed heterozygosities (Ho) to assess the genetic diversity of each of the three populations.

3. Results

3.1. Genomic relatedness and characterization of runs of homozygosity

Genomic pairwise kinship matrix between the six D’man rams was estimated using 3530 SNPs that are almost in linkage equilibrium with each other (Table 1). High positive kinship coefficients (> 0.11) were observed between ram DMAN1 on one hand and rams DMAN4 and DMAN6 on the other hand indicating that these individuals share more DNA identical by state than the rest of rams. Moderate kinship estimate was found between rams DMAN2, DMAN3 and DMAN5 (ranging from 0.011 to 0.0431). Importantly,negative kinship values were obtained between rams DMAN1, DMAN4 and DMAN6 on one hand and rams DMAN2, DMAN3 and DMAN5 on the other hand which indicates that these two groups are unrelated.

 Using our ROH definition (see Material and Methods), a total of 278 ROHs were identified in the six D’man rams. The number of identified ROHs per individual ranged between 36 (DMAN4) to 59 (DMAN 3). Total sum of ROHs per individual ranged between 71.28 Mb (DMAN4) to 224.22 Mb (DMAN3) (Table 2). A high positive correlation (0.98) was observed between the total sum of ROHs and the number of long ROHs (> 8 Mb). In all individuals, the highest number of ROHs was found within the [1-2 Mb] category. All individuals presented at least one ROH > 10 Mb. The number of ROHs per chromosome ranged between one and six (Figure 1). A high correlation was found between the number of ROHs and chromosome length (correlation = 0.90). Approximately 5.5%, on average, of the entire genome of the D’man individuals was covered by ROH segments. The highest coverage was shown by OAR24 (24.51%) and OAR11 (10.32%) (Figure 1).

Ten ROH islands were identified within the D’man rams located on OAR1 (2), OAR3 (2), OAR4 (2), OAR9 (1), OAR17 (1), OAR19 (1) and OAR25 (1). These ROH islands included 375 genes.

Two significant (enrichment score (ES) > 1.3) functional clusters were detected in the annotation clustering analysis (Table S1). The one with the highest enrichment score (ES=2.9) was related to Calcium-activated chloride channel protein and included chloride channel accessory (CLCA) genes. Proteins of the CLCA family may contribute to parasite expulsion through mucus hydration across the gut epithelium and smooth muscle contraction (Rowe et al. 2009).

The second functional cluster was related to arylesterase activity (ES=2.5) and included Paraoxonase 1 (PON1), Paraoxonase 2 (PON2) and Paraoxonase 3 (PON3) genes. PON1 enzymes are present in ovarian follicular fluid. Mutations in PON1 affect embryo quality among women undergoing in vitro fertilization (Mackness and Durrington 1995).

3.2. Genetic diversity and population structure

Mean observed heterozygosity was similar between D’man, Lacaune and WTT breeds (0.367, 0.368 and 0.375, respectively). In contrast, mean expected heterozygosity was clearly lower for D’man compared to Lacaune and WTT (0.399 Vs 0.492 and 0.494, respectively).

ADMIXTURE results for ancestral populations K=2 to K=4, for D’man, WTT and LAC, are shown in Figure 2. Assuming two ancestral populations (K = 2), D’man and WTT shared substantial level of ancestry (on average, 95 % of the D’man genome was shared with 75 % of WTT genome) while only 10 % genome was common between D’man and LAC. For K= 3, Lacaune was subdvided into two distinct genetic clusters while for K =4, half of D’man rams showed substantial signals of WTT genome introgression. These trends were confirmed by the DAPC (Figure 3). Indeed, Lacaune individuals were assigned to two different genetic clusters (1 and 3) while the three D’man individuals that showed signals of WTT introgression were assigned to the same genetic cluster as WTT (cluster 4). The three other individuals were assigned to cluster 2.

4. Discussion

The main aim of the present study was to assess the genomic inbreeding and the population structure of the D’man rams in the breeding state farm of Chenchou which is the leading supplier of rams for most of the farmers in Southern Tunisia. The importance of this study stems from the inability to introduce new individuals or use ram semen from Morocco to reduce inbreeding in the Tunisian flocks (owing to sanitary reasons). Therefore, internal genetic management would be the only way to tackle the undesirable consequences of the increasing inbreeding levels in Tunisian D’man. Despite the possibility to estimate kinship based on pedigree information, relatedness between individuals is more accurately estimated using molecular information. Indeed, variation of kinship due to Mendelian sampling is better captured using genomic information than the pedigree-based method. In the present study, we purposely used Loiselle estimator because it is less biased and more accurate than most of the other estimators when the sample is small (Wang 2017). The kinship matrix estimated from unlinked SNP markers showed that ram DMAN1 is closely related to rams DMAN4 and DMAN6 which implies that the former should not be mated with ewes from families 4 and 6. By contrast, rams DMAN2, DMAN3 and DMAN5 (which were younger than the three other rams) showed low kinship values (<0.043). Surprisingly, these three rams were shown to be unrelated with the first ones (i.e rams DMAN1, DMAN4 and DMAN6) although all the six rams are supposed to come from the same founder flock. This unexpected result was explained by the DAPC and admixture analysis showing significant levels of WTT introgression into the genome of these three rams (up to 75% according to ADMIXTURE). It is more likely that in order to tackle the growing rate of genetic problems arising from consanguinous matings, the managers of Chenchou station frequently practiced crosses with other local breeds. When we considered only non introgressed rams, expected heterozygosity dropped from 0.39 to only 0.26 thus becoming significantly lower than observed heterozygosity. This difference is mainly due to the high inbreeding rate between the three rams.

In the present study, we used runs of homozygosity to assess autozygosity for the six D’man rams. The mean number of identified ROH and mean total sum of ROHs per individual were 46 ± 8 and 139 ± 52 Mb, respectively. Using a similar ROH definition, Michailidou et al. (2018) found significantly lower values in the greek Chios (n=11.88), Boutsko (n=11.59) and Karagouniko (n=3.1) breeds. Total sum of ROHs per individual found by these authors did not exceeded 60 Mb (Michailidou et al. 2018). Similarly, Mastrangelo et al. (2017) found that mean number of ROH per individual in the Italian Barbaresca breed was 15.9 which is also by far lower than the average of D’man individuals. Among the six rams, DMAN3 and DMAN6 showed the highest number of long (> 8 Mb) ROHs (9 and 6 ROHs, respectively). Long autozygous ROH segments have long been implicated as a cause of the higher prevalence of genetic disorders stemming from the accumulation of strongly deleterious mutations as well as mildly deleterious variants (Szpiech et al. 2013). Under the assumption that the expected length of an autozygous segment follows an exponential distribution with mean equal to 1/2 g Morgans, where g is the number of generations since the common ancestor (Howrigan et al. 2011), we can hypothesize that the autozygous segments of these two rams originated from a common ancestor 6-7 generations ago (assuming that 1 cM ~ 1 Mb). Assuming that the generation interval is 2.5 years, this common ancestor lived about 15-17.5 years ago. Our results confirm the trend that ROHs occur more frequently on large chromosomes. The longest ROH was found on OAR03 (19 Mb), however, OAR24 and OAR11 showed the highest proportion in coverage by ROH segments. On the other hand, only DMAN6 showed a ROH on OAR24 while all but ram DMAN5 showed a ROH on OAR11. Among these, rams DMAN1, DMAN2 and DMAN6 showed a ROH beginning at 41.38, 46.66 and 43.37 Mb, respectively. It is worth noting that the major prolificacy locus FecL which was shown to have a significant effect on D’man litter size (Ben Jemaa et al. 2019) is located on the OAR11 (36.929322 Mb -36.992982 Mb). Both observations lead us to hypothesize that these ROHs resulted from a positive selection surrounding the FecL locus. Indeed, Ben Jemaa et al. (2019) found that the prolific allele, named FecL^L exhibited a high frequency in the D’man population (0.65) owing to a decades-old breeding strategy based on the selection of ewe lambs born from large litter size which are most likely born from either L+ or LL parents. Such positive selection could have generated long ROH surrounding the FecL locus.

Among the ten identified ROH islands in the D’man rams, the one located on OAR4 (8.363 Mb-13.125 Mb) is particularly interesting. This ROH island included PON1, PON2 and PON3 genes. It was hypothesized that PON genes are involved in follicle maturation as it was observed that an increase in the activity of these genes is positively correlated with follicle size in humans (Meijide et al. 2017). Interestingly, only non introgressed D’man rams (DMAN1, DMAN4 and DMAN6) showed patterns of autozygosity islands at this position. Therefore, we can hypothesize that this ROH island may have occurred due to selection for traits related to D’man prolificacy. Furthermore, by examining individually the list of genes located in the ten ROH islands, we identified two interesting genes located on OAR25 within the [4 Mb-8.609 Mb] interval. The first one is Interferon Regulatory Factor 2 Binding Protein gene (IRF2BP2). This gene is involved in fleece phenotype in sheep (Demars et al. 2017). The second interesting gene is the Beta-1,3-N-Acetylgalactosaminyltransferase 2 (B3GALNT2). Mutations in this gene were shown to be associated with stillbirth and dystocia in horse (Ducro et al. 2015).

5. Conclusion

To our knowledge, this is the first study to describe the population structure of the D’man breed based on medium-density SNP array data. We show that the youngest D’man rams of the breeding station of Chenchou display consistent introgression signals from the Tunisian WTT. On the other hand, the oldest rams have preserved the ancestral genomic structure of the breed but have high average level of inbreeding. The introduction of new rams from flocks located in the Northern part of the country can remedy this problem. However, a molecular verification of racial authenticity of the selected rams is essential. Finally, the ROH analysis suggests the possibility of presence of other genes influencing the D’man prolificacy.

Acknowledgements

This research was supported by the International Foundation for Science (IFS grant B/5478-2). The authors wish to thank Tunisian Livestock and Pasture Office (OEP) for their precious help in blood sample collection.

6. References

Alexander DH, Novembre J, Lange K (2009) Fast model-based estimation of ancestry in unrelated individuals. Genome Res 19 :1655–1664. doi:101101/gr094052109

Bailleul D, Stoeckel S, Arnaud‐Haond S (2016) RClone: a package to identify MultiLocus Clonal Lineages and handle clonal data sets in r. Methods in Ecology and Evolution 7 : 966–970. doi:101111/2041-210X12550

Ben Jemaa S, Ruesche J, Sarry J, Woloszyn F, Lassoued N, Fabre S. (2019) The high prolificacy of D’man sheep is associated with the segregation of the FecL ^L mutation in the B4GALNT2 gene. Reprod Domest Anim. doi:10.1111/rda.13391

Demars J, Cano M, Drouilhet L, Plisson-Petit F, Bardou P, Fabre S, Servin B, Sarry J, Woloszyn F, Mulsant P, Foulquier D, Carrière F, Aletru M, Rodde N, Cauet S, Bouchez O, Pirson M, Tosser-Klopp G, Allain D (2017) Genome-Wide Identification of the Mutation Underlying Fleece Variation and Discriminating Ancestral Hairy Species from Modern Woolly Sheep. Mol Biol Evol 34 : 1722–1729. doi:101093/molbev/msx114

Drouilhet L, Mansanet C, Sarry J, Tabet K, Bardou P, Woloszyn F, Lluch J, Harichaux G, Viguié C, Monniaux D, Bodin L, Mulsant P, Fabre S (2013) The Highly Prolific Phenotype of Lacaune Sheep Is Associated with an Ectopic Expression of the B4GALNT2 Gene within the Ovary. PLOS Genetics 9, e1003809. doi:101371/journalpgen1003809

Ducro BJ, Schurink A, Bastiaansen JWM, Boegheim IJM, van Steenbeek FG, Vos-Loohuis M, Nijman IJ, Monroe GR, Hellinga I, Dibbits BW, Back W, Leegwater PAJ (2015) A nonsense mutation in B3GALNT2 is concordant with hydrocephalus in Friesian horses. BMC Genomics 16 : 761. doi:101186/s12864-015-1936-z

Howrigan DP, Simonson MA, Keller M C (2011) Detecting autozygosity through runs of homozygosity: A comparison of three autozygosity detection algorithms. BMC Genomics 12 : 460. doi:101186/1471-2164-12-460

Jombart T (2008) adegenet: a R package for the multivariate analysis of genetic markers Bioinformatics 24 : 1403–1405. doi:101093/bioinformatics/btn129

Kdidi S, Calvo JH, González-Calvo L, Sassi M, Khorchani T, Yahyaoui M H (2015) Genetic relationship and admixture in four Tunisian sheep breeds revealed by microsatellite markers. Small Ruminant Research 131 :64–69. doi:101016/jsmallrumres201508012

Kijas JW, Lenstra JA, Hayes B, Boitard S, Porto Neto LR, San Cristobal M, Servin B, McCulloch R, Whan V, Gietzen K, Paiva S, Barendse W, Ciani E, Raadsma H, McEwan J, Dalrymple B, International Sheep Genomics Consortium Members (2012) Genome-wide analysis of the world’s sheep breeds reveals high levels of historic mixture and strong recent selection. PLoS Biol 10, e1001258. doi:101371/journalpbio1001258

Mackness M I, and Durrington P N (1995) HDL, its enzymes and its potential to influence lipid peroxidation Atherosclerosis 115 : 243–253

Martin P, Raoul J, Bodin L (2014) Effects of the FecL major gene in the Lacaune meat sheep population Genet Sel Evol 46 : 48. doi:101186/1297-9686-46-48

Mastrangelo S, Ciani E, Sardina M T, Sottile G, Pilla F, Portolano B, and Bi.Ov. Ita Consortium (2018). Runs of homozygosity reveal genome-wide autozygosity in Italian sheep breeds. Anim Genet 49 : 71–81. doi:10.1111/age.12634

Mastrangelo S., Portolano B., Di Gerlando R., Ciampolini R., Tolone M., Sardina M. T., and International Sheep Genomics Consortium (2017). Genome-wide analysis in endangered populations: a case study in Barbaresca sheep. Animal 11 : 1107–1116. doi:10.1017/S1751731116002780

Meijide S, Pérez-Ruiz I, Hernández ML, Navarro R, Ferrando M, Larreategui Z, Ruiz-Sanz J-I, Ruiz-Larrea MB (2017) Paraoxonase activities in human follicular fluid: role in follicular maturation Reproductive BioMedicine Online 35 : 351–362. doi:101016/jrbmo201706008

Michailidou S, Tsangaris G, Fthenakis G C, Tzora A, Skoufos I, Karkabounas S C, Banos G, Argiriou A, and Arsenos G (2018) Genomic diversity and population structure of three autochthonous Greek sheep breeds assessed with genome-wide DNA arrays Mol Genet Genomics 293 : 753–768 doi:101007/s00438-018-1421-x

Miller S A, Dykes D D, and Polesky H F (1988) A simple salting out procedure for extracting DNA from human nucleated cells. Nucleic Acids Res 16 :1215

Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira M A R, Bender D, Maller J, Sklar P, de Bakker P I W, Daly M J, and Sham P C (2007) PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses. Am J Hum Genet 81 : 559–575

Rekik M, Ben Salem I, Ben Hamouda M, Diallo H, Ammar H, Aloulou R (2005) Numerical and weight productivities of crossbred D’man by local Queue Fine de l’Ouest ewes Revue Elev Méd vét Pays trop 58 : 81–88

Rosenberg N A (2004) distruct: a program for the graphical display of population structure. Molecular Ecology Notes 4 : 137–138 doi:101046/j1471-8286200300566x

Rowe A, Gondro C, Emery D, Sangster N (2009) Sequential microarray to identify timing of molecular responses to Haemonchus contortus infection in sheep. Veterinary Parasitology 161 : 76–87 doi:101016/jvetpar200812023

Szpiech ZA, Xu J, Pemberton TJ, Peng W, Zöllner S, Rosenberg NA, Li JZ (2013) Long Runs of Homozygosity Are Enriched for Deleterious Variation. Am J Hum Genet 93 : 90–102 doi:101016/jajhg201305003

Wang J (2017) Estimating pairwise relatedness in a small sample of individuals. Heredity (Edinb) 119 : 302–313. doi:101038/hdy201752