Background {#Sec1}
==========

The Djallonke sheep is recognised for its natural ability to withstand a harsh, hot and humid tropical climate, where it is faced with the challenges of persistent drought, diseases and feed scarcity \[[@CR1], [@CR2]\]. Adaptation is probably a consequence of natural selection over several millennia \[[@CR3]--[@CR5]\]. Genomic regions adjacent to loci under adaptive selection over time are usually characterised by low heterozygosity \[[@CR6]\]. The most important livestock diseases are trypanosomiasis and haemonchosis \[[@CR7]--[@CR10]\]. The natural ability of the Djallonke to survive and remain productive under trypanosome challenge with very low mortality and without the aid of trypanocidal drugs is referred to as trypanotolerance \[[@CR11]\]. Trypanosomiasis in sub Saharan Africa is estimated to cause annual losses of more than 4.5 billion dollars (US\$) through direct and indirect production costs \[[@CR12], [@CR13]\]. Development of trypanotolerance is considered to be the most economical and sustainable option for combating African trypanosomiasis \[[@CR9], [@CR14], [@CR15]\]. The potential of this trypanotolerant trait in mitigating the disease in Africa has recently been reviewed \[[@CR16]\]. However, because Djallonke sheep have a relatively small mature body weight (between 20 kg to 30 kg \[[@CR17]\];) farmers often cross-breed them with the larger, but more disease susceptible, Sahelian breed.

In spite of the importance of the Djallonke and Sahelian sheep to the region, genetic studies are scarce. There are no records of whole genome variant characterisation in either the Djallonke or Sahelian breeds, besides our preliminary report at the International society of animal genetics conference \[[@CR18]\]. The objectives of this study are: i) to identify and document variants in each breed and ii) to use these to investigate putative candidate genetic regions in both breeds.

Methods {#Sec2}
=======

Animals {#Sec3}
-------

Four ewes and one ram of the Djallonke breed (DJ) were selected from the National Open Nucleus Breeding Station (ONBS) dedicated for Djallonke Sheep in Ejura in the Ashanti region (longitude 01^o^ 28′W and latitude 06^o^ 41′ N). Five Sahelian (SA) ewes were selected from the National Sheep Breeding Station in Pong Tamale (longitude 00^o^ 54′W and latitude 09^o^ 38′ N). All sheep were reproductively mature (13--24 months old) and were chosen in consultation with the management of the breeding stations to represent unrelated animals that were true to breed type (phenotypically similar to the breed ideal). The management relied on stock records for determination of relatedness among sheep on two breeding stations. Approximately 9 ml of blood was collected via the jugular vein into disodium EDTA vacutainers. All sampled sheep were monitored on farm for at least 24 h post sampling and no adverse effect was recorded. No sheep were sacrificed during this study. The samples were transported at 0^o^ -- 4^o^ C to the laboratory, centrifuged at 800 x g for 3 min at room temperature (15-25 °C) with the rotor bucket brake off. The buffy coat was used immediately for genomic DNA extraction or was stored at -20 °C. Genomic DNA was extracted from each of the buffy coat samples using the Zymo Quick-gDNA™ MiniPrep DNA purification Kit (according to the manufacturer's protocol). DNA quality and concentration were assessed using agarose gel electrophoresis (1% in 1xTAE) and by Nanodrop spectrophotometry.

Library construction and sequencing {#Sec4}
-----------------------------------

For each individual, 100 ng of DNA was sheared using the Covaris S2 System, to generate a broad range of DNA fragments with sizes from 100 to 1000 bp. The DNA fragments were ligated to T-overhang adaptors with the NEB Next Ultra kit (New England Biosciences). Each animal had a unique barcode (Ion Xpress Barcodes, Life Technologies). Fragments of approximately 300-330 bp were size-selected using the E-gel system (Invitrogen), and recovered fragments were further purified using AMPure XP SPRI beads (Beckman). Equimolar amounts of each library were combined and amplified using an Ion Chef system (ThermoFisher Scientific) via emulsion PCR, then sequenced on an Ion Proton™ system (ThermoFisher Scientific) using a PI chip. Genomic DNA from the ten individuals was separately sequenced, and the sequencing reads were then pooled by breed. After filtering and trimming, an average of 10 and 13% of the reads were excluded due to low quality, and 26 and 28% excluded due to polyclonality for the Djallonke and Sahelian samples, respectively. Coverage analysis was performed on a total (post QC) of 73 Gbp of sequenced data, comprising 404,755,012 pooled reads (average read length 185.4 nucleotides) and 57.6 Gbp comprising 303,136,043 pooled reads (average read length 176.6 nucleotides) for Djallonke and Sahelian sheep, respectively. The genome coverage depth obtained was calculated as being 27.90x and 22.01x for the Djallonke and the Sahelian respectively, and covered 97% of sheep reference assembly v3.1. All the variants were submitted to the European variant archive of the European Bioinformatics Institute with the accession number PRJEB15642.

Mapping and pre-processing of reads {#Sec5}
-----------------------------------

Base calling, de-multiplexing, quality control (QC) and alignment pre-processing \[[@CR19]\] was completed using Torrent Suite 4.6 on a Torrent Server (ThermoFisher Scientific). Briefly, polyclonal and uniformly low-quality reads were removed, and the remaining reads were trimmed from the 3′ end only. Mapping was also performed within Torrent Suite 4.6, using the Torrent mapping alignment program (TMAP). Individual libraries were mapped to the sheep reference genome Oar v3.1 (University of California, Santa Cruz (UCSC)). For each of the sheep breeds, all individual BAM files were merged and sorted using SAMtools v0.1.19-44,428 cd \[[@CR20]\], and coverage analysis was performed for both the individual and combined datasets through automated plugins in TorrentSuite 4.6. Duplicate reads were removed using Picard Tools v1.122.

Variant calling pipeline {#Sec6}
------------------------

Genome Analysis Tool Kit version 3.2.2 (GATK) RealignerTargetCreator and IndelRealigner were used to produce realignments of the pooled BAM files for each breed. GATK HaplotypeCaller was used in GVCF mode to call intermediate genome-wide variants separately for the pooled DJ genomes and pooled SA genomes, producing two pooled genomic variant call format (gvcf) files. GATK GenotypeGVCFs was then used to perform a joint genotyping of the two pooled gvcf files with minimum standard confidence thresholds for both calling and emitting variants set at 30 to produce a composite pooled variant call format (vcf) file (Pooled-Sheep VCF). This analysis was selected to ensure good quality variant calling and reduce false discovery rates. Finally, VCFtools (v0.1.15) \[[@CR21]\] was used to extract individual Djallonke and Sahelian samples from the composite joint genotyped vcf into separate vcf files, which were used for downstream analyses.

Genetic relationship matrix {#Sec7}
---------------------------

To determine genetic relatedness, a genomic relationship matrix (GRM) was computed on a composite vcf file that contained all the 10 vcf files generated from both the Djallonke and the Sahelian samples, using the Genome-wide Complex Trait Analysis (GCTA) software \[[@CR22], [@CR23]\]. Furthermore, Principal Components analysis (PCA) was used to compare the autosomal genomes of all individual samples to determine population substructure for each breed \[[@CR24]\], and assess the genetic relatedness using the SNPRelate implemented in R CRAN (<http://cran.r-project.org>) \[[@CR25]\].

Detection of regions of reduced heterozygosity {#Sec8}
----------------------------------------------

HomSI (Homozygosity Stretch Identifier) was used to identify regions of reduced heterozygosity in both genomes \[[@CR26]\]. The Djallonke genome was designated as the index case and compared against the Sahelian as the unaffected case for input settings for the HomSI analysis. Analysis of runs of homozygosity in Djallonke and Sahelian using HomSI, with the stringent settings of 5 Mb window size and 10 kb sliding size, allowed the capturing of a wide spectrum of different lengths of homozygosity throughout the genome \[[@CR26]--[@CR29]\]. Integrative Genomics Viewer (IGV 2.3.46, [www.BroadInstitute.org](http://www.broadinstitute.org)) was used to view vcf and BAM file tracks aligned to the sheep reference genome Oar v3.1, selecting regions based on genomic coordinates of regions of contrasting reduced heterozygosity identified by HomSI in order to identify candidate genes. Prominent regions were investigated using IGV for the specific genes of interest \[[@CR30], [@CR31]\] and were used to identify candidate genes within the region. For every prominent region of low heterozygosity in the HomSI output, the co-localised candidate gene or genes were inferred from the available Ensembl (Release 85 and 86) annotated sheep reference assembly (version 3.1) or from the conserved synteny for other mammalian genomes from the Ensembl genome database \[[@CR32]--[@CR34]\]. The approach presented here for investigating genetic evidence of trypanotolerance and resistance to haemonchosis was by directly linking regions with contrasting heterozygosity in this dataset to the reported candidate genes in the database for Animal Quantitative Trait Loci (Animal QTLdb) \[[@CR35]\] and other previously published genetic association studies for the two traits. There have been several previous genomic investigations of resistance to nematode infection including *H. contortus* in multiple sheep breeds and these results were compared with our Djallonke and Sahelian sheep results \[[@CR36]--[@CR40]\]. In contrast, there has been no previous genomic investigation of trypanosomiasis in any sheep breed; therefore, comparison was made with the reported trypanotolerance associated candidate genes in Ndama cattle \[[@CR41]--[@CR45]\].

Annotation and functional analysis of genomic variants {#Sec9}
------------------------------------------------------

As there were unequal numbers of males and females used between the two breeds, for the purpose of a balanced comparison, autosomal chromosomes were extracted from each of the pooled vcf files using VCF tools v0.1.15. Known SNPS were annotated in the vcf files with SnpSift v4.2 (annotate command) using the Ensembl Release-85 Variation reference vcf for *Ovis aries* as the database. ([ftp.ensembl.org/pub/release-85/variation/vcf/ovis_aries/](ftp://ftp.ensembl.org/pub/release-85/variation/vcf/ovis_aries/)) \[[@CR32]\]**.** SnpEff v4**.**2 (Cingolani et al., 2012) was used for functional annotation of identified autosomal SNPs in both the Djallonke and Sahelian genomes based on the Ensembl sheep genome assembly Oar_v3.1.82. Pairwise comparison of Genomic SNPs and INDELS for Sahelian and Djallonke sheep was computed using the BEDTools suite v. 2.26.0 \[[@CR46]\].

Results {#Sec10}
=======

Genetic relationship matrix and principal component analysis {#Sec11}
------------------------------------------------------------

The GRM computed for these datasets supports the assumption that the ten individual animals sampled were unrelated. Additional file [1](#MOESM1){ref-type="media"}shows the GRM output for this analysis. The PCA computed for the 10 individual datasets from the two sheep breeds showed distinct clustering for the two breeds (Fig. [1](#Fig1){ref-type="fig"}). Fig. 1PCA of Autosomal Genomes for 5 Djallonke (DJ) and 5 Sahelian (SA) sheep showing distinct clustering of the two breeds

Comparison of genomes {#Sec12}
---------------------

On average, there was one variant every 191 base pairs (bp) in the Djallonke sheep and 1 variant every 193 bp in the Sahelian sheep (Table [2](#Tab2){ref-type="table"}). Transition to transversion ratios were similar in the Djallonke (2.47) and Sahelian (2.48) sheep (Table [1](#Tab1){ref-type="table"}). The estimated missense to silent SNP ratios were also similar between the Djallonke and Sahelian sheep (0.69 and 0.68 respectively). Similarly, the two datasets had equal insertion to deletion ratios (0.38). Similar proportions of SNPs, insertions and deletions were observed for Djallonke (86%: 4%: 10%) and Sahelian (87%: 4%: 9%) sheep. Approximately 84% of variants in both Djallonke and Sahelian sheep were present in the Ensembl Variation database (release 85) hence approximately 16% are unidentified variants (Table [1](#Tab1){ref-type="table"}). Analysis of only SNPs, however, revealed that approximately 94% were present in the database (Table [1](#Tab1){ref-type="table"}), indicating that approximately 6% of SNPs were previously unreported in sheep. Table 1Summary of SnpEff variant analysis of Djallonke and Sahelian sheep genomesDjallonke all variantsSahelian all variantsDjallonke SNP onlySahelian SNP onlyNumber of variants12,821,83612,654,76111,100,61910,978,689Number of known variants10,764,740 (83.956%)10,671,465 (84.328%)10,460,303 (94.232%)10,372,765 (94.481%)Number of multi-allelic entries9235919955465534Number of effects14,652,59014,454,27812,641,43412,496,971Variant rate1 per 191 bases1 per 193 bases1 per 220 bases1 per 223 basesTs/Tv ratio2.47362.48192.47362.4819

The distribution of variants by chromosome {#Sec13}
------------------------------------------

The distribution of the variants (i.e. the sum of SNPs, insertions and deletions) was similar across all the autosomes, with chromosomes 11 and 26 having the lowest and highest frequencies for variants in both breed genomes (Table [2](#Tab2){ref-type="table"}), respectively. A total of 12,821,836 and 12,654,761 variants were identified in Djallonke and Sahelian sheep, respectively, with 12,556,638 (96.30%) shared between the two genomes. In total, 324,760 (2.49%) variants were specific to the Djallonke breed, whereas 158,085 (1.21%) variants were specific to the Sahelian breed. Therefore, the total number of variants identified for the two breeds is 13,039,483. Analysis with BEDTools intersect indicated that 242,572 SNPs and 82,609 indels were specific to the Djallonke, whereas 120,652 SNPs and 37,762 indels were unique to the Sahelian breed. Table 2Comparison of the Djallonke and the Sahelian GenomesChromosomeLength (bases)DjallonkeSahelianVariantsRatio of nucleotides to variantsVariantsRatio of nucleotides to variants1275,612,8951,415,8831941,397,7991972248,993,8461,252,1691981,236,4712013224,283,2301,126,6471991,112,2392014119,255,633614,660194606,6521965107,901,688547,117197540,3041996117,031,472638,104183629,7071857100,079,507521,899191514,989194890,695,168463,424195457,602198994,726,778508,149186501,4521881086,447,213458,020188451,9261911162,248,096309,854200305,9052031279,100,223418,404189412,8991911383,079,144420,072197414,5772001462,722,625313,626199309,4452021580,923,592431,999187426,2461891671,719,816397,765180392,4981821772,286,588390,027185384,8751871868,604,602361, 551189356,6431921960,464,314314,693192310,7771942051,176,841272,206188268,3691902150,073,674279,334179275,2261812250,832,532279,610181275,5331842362,330,649342,228182337,6401842442,034,648233,669179230,6261822545,367,442259,240175256,1301772644,077,779251,486175248,131177Total2,452,069,99512,821,83619112,654,761193

Distribution of autosomal SNPs by genomic region {#Sec14}
------------------------------------------------

Most SNPs in both sheep breeds were intergenic or intronic (Table [3](#Tab3){ref-type="table"}), with approximately 1% located in the remaining genic regions (i.e. untranslated regions (UTR), exons, and splice sites). Although the Djallonke sheep had a higher number of SNPs than the Sahelian sheep, the ratios of the SNPs in these three regions (intergenic, intronic and "other" genic regions including exons) were similar for Djallonke (68.78%: 30.04%: 1.18%) and Sahelian sheep (68.81%: 30.02%: 1.17%). A comparison of the "other" genic category revealed similar proportions of synonymous, non-synonymous, splice site, UTR and miscellaneous variants for each breed (Table [3](#Tab3){ref-type="table"}). Table 3Comparison of Autosomal SNPs by type and functional categoryVariant typesDjallonkeSahelianFunctional class Missense33,984 (40.7%)33,395 (40.5%) Nonsense344 (0.41%)338 (0.41%) Silent49,171 (58.9%)48,732 (59.1%) Intergenic region8,861,5168,747,769 Upstream gene611,224601,689 Downstream gene605,975596,609 Intronic4,413,0524,350,066 Intragenic199181 3 prime UTR29,54828,988 5 prime UTR681672 5 prime UTR Truncation22 5 prime UTR54125307 Splice acceptor15071454 Splice donor16401553 Splice region14,54414,262 Exon Synonymous49,14148,702 Non synonymous33,90233,313 Non coding exon22,86522,315 Non coding transcript5252 Exon loss1212 Initiator codon78 Start lost9792 Stop gained541528 Stop lost4443 Stop retained2828 Transcript ablation11 In frame deletion349330 In frame insertion674646 Disruptive in frame deletion476447 Disruptive in frame insertion607592 Frameshift16,99716,661

Regions of homozygosity {#Sec15}
-----------------------

Approximately 2.5 Gbp of autosomal chromosomal DNA was resolved into about 50,000 detection windows for each breed. HomSI analysis, identified regions having reduced heterozygosity (Fig. [2](#Fig2){ref-type="fig"}; blue) of various sizes (1 to \> 100 kb) across genic and intergenic regions within both breeds. Seventy of these reduced heterozygosity regions co-localised with known candidate genes. There were also several genic and intergenic regions that showed reduced heterozygosity but did not contain any known putative candidate gene. In addition, there were regions for which the Djallonke show complete fixation for one allele (blue) and the Sahelian showing complete fixation for the alternative allele (white) e.g. *TRHDE* (Fig. [2](#Fig2){ref-type="fig"}). Fig. 2Comparison of the HomSI analysis of three genes associated with adaptive traits for Djallonke and Sahelian sheep genomes

Putative regions for tolerance to trypanosomiasis {#Sec16}
-------------------------------------------------

Eight regions of reduced heterozygosity were observed to be co-localised with previously reported trypanotolerance associated genes (Table [4](#Tab4){ref-type="table"}) \[[@CR42], [@CR47], [@CR48]\]. Six of the eight reported genes (*CTSS, ARHGAP15, INHBA, STX7, RAB35, CD19*) were co-localised with regions of reduced heterozygosity in the Djallonke genome only, and the other two (*SCAMP1, TICAM1*) were co-localised with regions of reduced heterozygosity in both Djallonke and Sahelian genomes (Fig. [3](#Fig3){ref-type="fig"}). The Djallonke sheep show longer runs of reduced heterozygosity (blue) than the Sahelian sheep at the *INHBA* and *RAB35* gene regions, but both breeds show reduced heterozygosity across approximately 2-kb (16,922 kb--16,924 kb) of the TICAM1 gene (Fig. [3](#Fig3){ref-type="fig"}). The Sahelian breed shows increased heterozygosity (orange) between 9279 kb to 9286 kb in the *SCAMP1* region. In contrast, the Djallonke shows reduced heterozygosity within the same region (blue). Table 4Trypanotolerance candidate genes co-localised with regions of reduced heterozygosity (RORH) in Djallonke and Sahelian sheep and the reported orthologs in the cattle genomeSheep BreedGenomic region (bp) co-localised with RORH (Oar_v3.1)Candidate genesOrthologous loci in Cattle UMD\
v 3.1Djallonke1:99,623,159-99,650,475*CTSS*3:20,024,302-20,047,228Djallonke2:165,154,368-165,495,812*ARHGAP15*2:53,065,587-53,732,838Djallonke4:79,355,697-79,366,930*INHBA*4:79,986,254-79,997,754Both7:9,272,393-9,359,760*SCAMP1*10:9,369,310-9,520,700Both5:16,922,626-16,924,761*TICAM1*7:20,547,964-20,550,264Djallonke8:57,581,762-57,631,158*STX7*9:71,381,757-71,455,585Djallonke17:62,105,682-62,115,444*RAB35*17:64,724,244-64,742,928Djallonke24:25,905,178-25,915,071*CD19*25:164,039-26,169,956 Fig. 3Comparison of the HomSI analysis of four Trypanotolerance associated genes for Djallonke and Sahelian sheep genomes

Putative regions for resistance to Haemonchosis {#Sec17}
-----------------------------------------------

There were 25 regions of reduced heterozygosity that co-localised with previously reported candidate gene for Haemonchosis resistance (Table [5](#Tab5){ref-type="table"}). Twenty-one of the regions had reduced heterozygosity in the Djallonke sheep. The remaining four regions *(MHCII-DRB1, PIK3CD, MUC15, IL17RB*) had reduced heterozygosity in both Djallonke and Sahelian sheep. Figure [4](#Fig4){ref-type="fig"} shows three regions associated with resistance to *Haemonchus contortus* infection: the *IFNG* gene, the *CHIA* gene, and the *SUGT1* gene. In each case, only the Djallonke sheep displayed reduced heterozygosity. Table 5Regions of reduced heterozygosity (RORH) in Djallonke and Sahelian sheep co-localised with reported candidate genes for resistance to nematodesSheep breedGenomic region (bp) co-localised with RORH (Oar_v3.1)Candidate GeneTrait InferenceReferenceBoth20:25,400,738-25,402,966*MHCII-DRB1*Gastrointestinal nematodesSchwaiger et al. (1995)Djallonke1:27,524,283-27,601,025*LRP8H. contortus*Benavides et al. (2016)Djallonke1:87,657,990-87,674,113*DENND2DH. contortus*McRae et al. (2014)Djallonke1:87,710,082-87,728,410*CHI3L2H. contortus*McRae et al. (2014)Djallonke1:87,788,905-87,811,255*CHIAH. contortus*McRae et al. (2014)Djallonke3:151,527,165-151,535,188*IFNG*Mixed intestinal parasitesColtman et al., 2001Djallonke3:123,851,175-125,982,479*ATP2B1H. contortus*Benavides et al. (2016)Djallonke8:62,006,022-62,039,859*IL20RAH. contortus* + othersPeriasamy et al., 2014Both12:41,923,922-41,973,979*PIK3CDH. contortus* + othersPeriasamy et al., 2014Djallonke12:62,163,057-62,283,591*LAMC1H. contortus*Benavides et al. (2016)Both15:55,404,807-55,417,235*MUC15H. contortus*Benavides et al. (2016)Djallonke17:52,168,248-52,191,479*ABCB9H. contortus*Yang et al., 2015Djallonke10:11,465,154-11,505,274*SUGT1H. contortus*Yang et al., 2015Djallonke14:48,062,306-48,071,138*PAK4H. contortus*Yang et al., 2015Djallonke14:14,206,784-14,215,316*FCER2H. contortus*Yang et al., 2015Djallonke4:44,668,137-45,205,142*RELNH. contortus*McRae et al., 2014Djallonke6:89,053,717-89,061,196*AREGH. contortus*Zhengyu et al., 2016Djallonke11:57,796,766-57,800,093*SOX9H. contortus*Benavides et al., 2016Djallonke6:70,189,729-70,234,612*KITH. contortus*Zhengyu et al., 2016Djallonke16:66,450,575-66,482,861*NSUN2H. contortus*McRae et al., 2014Both19:47,044,394-47,059,322*IL17RBH. contortus*Zhengyu et al., 2016Djallonke19:55,173,564-55,175,033*HRH1H. contortus*McRae et al., 2014Djallonke25:44,535,996-44,543,183*CXCL12H. contortus*Zhengyu et al., 2016Djallonke19:53,290,059-53,291,096*CXCR6H. contortus*Zhengyu et al., 2016Djallonke20:34,055,201-34,055,659*UBE2NH. contortus*Benavides et al., 2016 Fig. 4Comparison of the HomSI analysis of three Haemonchotolerance associated genes for Djallonke and Sahelian sheep genomes

Putative regions for adaptation to tropical conditions {#Sec18}
------------------------------------------------------

There were also genomic regions of reduced heterozygosity that were co-localised with genes known to be associated with immune responses and natural adaptation (Fig. [2](#Fig2){ref-type="fig"}). A total of 37 candidate genes fell within these reduced heterozygosity genomic regions in the Djallonke sheep, including 14 that were shared with the Sahelian sheep (Table [6](#Tab6){ref-type="table"}). Three of the gene regions associated with adaptive selection (*MSRB3* gene, *APC2* gene, and *TRHDE* gene) are shown in Fig. [2](#Fig2){ref-type="fig"}. Differences between these genes were observed with respect to the polymorphism patterns. The Djallonke sheep have reduced heterozygosity (blue) at the region of the *MSRB3* gene whilst the Sahelian was more polymorphic (Fig. [2](#Fig2){ref-type="fig"}). Over much of the *TRDHE* gene region, the two sheep breeds were fixed for alternative alleles. Both sheep breeds, however, showed reduced heterozygosity for the same allele between coordinates 41,267,000 and 41,272,000, encompassing 7 exons (exons 16 to 22) of the *APC2* gene. Table 6Candidate genes for tropical adaptation co-localised with regions of reduced heterozygosity in Djallonke and Sahelian sheepSheepChr.Genomic coordinatesCandidate geneTrait InferenceReferenceDJ185,955,810-86,011,841*GNAI3*Melanogenesis (Thermo-tolerance)Kim et al., 2016DJ1188,388,916-188,441,236*LMLN*MelanogenesisKim et al., 2016DJ142,584,598-42,656,778*IL12RB2*Immune functionsRoffler et al., 2016DJ1168,393,395-168,624,986*ALCAM*Immune functionsRoffler et al., 2016DJ1121,075,675-121,168,117*SYNJ1*Phosphatidylinositol dephosphorylationRoffler et al., 2016Both252,423,842-52,445,175*NPR2*Skeletal Morphology and body sizeKijas et al., 2012DJ49,433,282-9,465,962*KRIT1*Regulation of endothelial cell proliferation and migrationRoffler et al., 2016DJ485,316,865-85,381,180*TSPAN12*Regulation of signal transduction of cell surface receptorsRoffler et al., 2016Both541,256,802-41,272,546*APC*Immune functions (Tumour suppressor)Roffler et al., 2016Both636,514,210-36,556,824*ABCG2*Urea Metabolism (Homeostasis)Kijas et al., 2012Djallonke694,584,400-94,605,575*FGF5*Regulation of fibroblast growth factor receptorKijas et al., 2012Both763,450,344-63,456,226*BMP4*Body size and developmentKim et al., 2016DJ3204,447,104-204,461,390*OLR1*Internalization, degradation of oxidised low density lipoprotein by endothelial cellsRoffler et al., 2016Djallonke3108,235,641-108,685,027*TRHDE*Regulation of appetite and digestionKim et al., 2016Djallonke3154,219,234-154,397,986*MSRB3*Regulation of response to oxidative stressKijas et al., 2012Both335,907,955--36,031,445*ALK*Immune function (Protein phosphorylation)Kim et al., 2016Both399,472,045-99,509,159*IL1R1*Immune functionKim et al., 2016Djallonke954,817,997-54,825,977*IL7*Immune functionKim et al., 2016Both1036,838,524-36,858,872*ATP12A*Homeostasis (Potassium and Sodium)Kim et al., 2016Djallonke1040,800,056-40,821,770*PCDH9*Homophilic cell adhesionKim et al., 2016Both1136,083,204-36,098,540*PDK2*HomeostasisKijas et al., 2012Both1118,245,395-18,411,418*NF1*HomeostasisKijas et al., 2012Both1378,815,423-78,893,076*NFATC2*Immune functionKijas et al., 2012DJ13666,266-1,154,524*PLCB1*ThermotoleranceKim et al., 2016DJ1545,551,281-45,552,222*OR2AG1*Response to stimulusKijas et al., 2012Both1638,969,273-39,028,126*PRLR*ReproductionKijas et al., 2012DJ1718,131,831-18,226,233*ELF2*Regulation of transcriptionKim et al., 2016DJ1729,240,707-29,257,289*PGRMC2*ReproductionKim et al., 2016DJ1819,723,286-19,802,578*ABHD2*Wound healingKijas et al., 2012DJ184,690,980-4,728,935*ALDH1A3*Energy, digestive MetabolismKim et al., 2016Both1838,107,388-38,110,333*FOXG1*Regulation of transcriptionKijas et al., 2012Both1931,583,789-31,811,540*MITF*MelanogenesisKijas et al., 2012DJ197,255,507-7,331,066*GLB1*Cellular metabolismKijas et al., 2012DJ1933,852,131-34,140,194*SUCLG2*Cellular metabolismKim et al., 2016Both2026,649,266-26,651,191*HSPA1A*HomeostasisBoth2149,011,232-49,012,130*IFITM21*Immune functionsRoffler et al., 2016DJ2142,526,284-42,531,851*BATF2*Immune functionsRoffler et al., 2016

Discussion {#Sec19}
==========

A comparison of the genomes of Djallonke and Sahelian sheep in this study has shown that Djallonke sheep have a genetic variant every 191 base pairs while the Sahelian sheep have a genetic variant every 193 base pairs. Approximately 16% of the variants had not been previously reported in sheep. The two breeds also had similar ratios of transitions to transversions (2.5), missense to silent mutations 90.7) and insertions to deletions (0.4). These breeds also had similar proportions of SNP to indels. The distribution of variants across the autosomal chromosomes was also similar; in both breeds chromosome 11 had the lowest frequency of variants while chromosome 26 had the highest frequency.

The transition to transversion ratio obtained for Djallonke (2.47) and Sahelian (2.48) are similar to the expected values observed for other mammalian genomes: 2.26 for cattle whole genome \[[@CR49]\], 2.13 for human intergenic SNPs \[[@CR50]\], and 2.81 for human exonic SNPs \[[@CR51]\]. These comparable ratios support the reliability of the sequenced datasets in this study, and they are therefore expected to contain low numbers of false positives (Type 1 errors) caused by random sequencing errors. This is further underscored by the high sequencing coverage statistics obtained for both genomes (i.e. \> 97% and \> 20x), which is suitable for "high-confidence" variant calling \[[@CR52]\]. The advantage that sequencing has over medium or high-density SNP genotype datasets, is that it provides higher resolution and power for the detection of selection signatures over relatively short distances \[[@CR53], [@CR54]\]. For instance, the Illumina Ovine 50KSNP BeadChip and Illumina Bovine HD 800KSNP BeadChip only provide a SNP density of approximately 1 SNP for every 5 million and 3 million bases, respectively. Furthermore, the use of markers on breeds that were not included in the training set for the marker development introduces further possible ascertainment bias into the analysis.

The proportion of SNPs in the intergenic, intronic and the remaining genic regions including exons for the two genomes are similar to the proportions recorded in Korean cattle breeds \[[@CR49]\]. The exonic regions, although containing the least number of SNPs, represent the most important subset of SNP, because they are more likely to be associated with changes in protein sequence, structure and function than intronic and intergenic SNPs \[[@CR55]\]. In particular, population-specific, rare exonic SNPs have been shown to be the most consequential determinants of fitness traits in humans \[[@CR56]\]. Fixed non-synonymous SNPs, which are described as SNPs for which only one allele (of a given locus) is present in a population, are of major interest in identifying breed or population specific traits \[[@CR49]\].

The high number of novel variants identified: 2,057,096 (16.03%) in the Djallonke breed and 1,983,296 (15.67%) in the Sahelian sheep confirms that these breeds are an important genetic resource for world sheep diversity. More than 0.5 million SNPs in each of the two sheep breeds are probably novel. There were also high numbers of breed specific variants; 242,572 SNP and 82,609 indels in the Djallonke and 120,652 SNP and 37,762 indels in the Sahelian breed. These breed specific variants could facilitate the sustainable management of these breeds and aid in confronting future emerging livestock diseases as well other global challenges, such as the uncertain consequences of climate change \[[@CR57]\]. Recent reports indicate that most of the indigenous African livestock breeds are endangered \[[@CR58]\] and might become extinct.

The HomSI scan permitted the identification of regions of reduced heterozygosity in greater detail than other sliding window algorithms such as the "Integrated haplotype homozygosity score (iHS)" \[[@CR59], [@CR60]\] and "the composite of likelihood ratio (CLR)" statistics \[[@CR53], [@CR61]\]. Furthermore, iHS detects only "ongoing sweeps" and CLR detects only "completed sweeps" in a target genome. Additionally, selection sweeps identified using HomSI are of higher resolution in comparison to the other methods (with sliding windows of 10,000 versus 50,000 base pairs). Runs of reduced heterozygosity were identified in reported candidate gene regions and may have resulted from selective sweeps. There were 70 regions of identified to have relatively reduced heterozygosity that co-localised with previously reported candidate genes for tolerance to trypanosomiasis, resistance to haemonchosis or adaptation to tropical conditions.

Five of the eight candidate trypanotolerance genes previously reported in a peripheral blood mononuclear cell gene expression study in experimentally infected trypanotolerant Ndama cattle \[[@CR48]\] fell within the regions of reduced heterozygosity identified in this study (*STX7, SCAMP1, RAB35, CD19, CTSS*). We identified putative selection signatures that co-localised with four of these five genes in Djallonke sheep, but the fifth candidate gene (*SCAMP1*) had similar values of heterozygosity in both Djallonke and Sahelian sheep. It is possible that Sahelian sheep may have also undergone some selection for trypanotolerance. Interestingly a sixth candidate gene, the INHBA gene, also fell within a region of reduced heterozygosity in the Djallonke. The INHBA candidate loci is the most significantly associated trypanotolerant loci reported in the Animal QTLdb to date \[[@CR35]\]. The *INHBA* gene was identified through fine mapping analysis of four a priori identified trypanotolerant associated loci in 360 Ndama cattle under natural infection conditions \[[@CR42]\]. The *INHBA* gene has been shown to regulate the differentiation of hematopoietic cells in mammals \[[@CR62]--[@CR65]\]. This is consistent with the hypothesised mechanisms of trypanotolerance, because the trait is strongly associated with the host's capacity to control anaemia \[[@CR4], [@CR42], [@CR66]\]. The last two trypanotolerance candidate genes, *ARHGAP15* and *TICAM1,* fell within regions of reduced heterozygosity. These two genes were previously identified in a combined transcriptomic and selective sweep analysis of infected trypanotolerant Ndama and Boran cattle \[[@CR43]\]. These genes co-localised within regions of reduced heterozygosity in the Djallonke dataset, whereas only *TICAM1,* but not *ARHGAP15,* was co-localised with a region of reduced heterozygosity in the Sahelian dataset.

Previous studies on trypanotolerance have used a lower density of molecular markers \[[@CR41], [@CR42], [@CR47], [@CR48]\], and the confidence limits of the reported candidate loci are quite large \[[@CR16]\]. Comparison of the Djallonke and Sahelian sheep revealed several putative selective sweeps of varying sizes (down to 2 Kilo-bases resolution). Although trypanotolerance is a complex quantitative trait and controlled by many genes, it is highly unlikely that all of the variants captured in these regions are causative variants. It is more likely that some variants are in linkage disequilibrium with the causal variants and hitch-hiked over time \[[@CR67]\].

Gene ontology (GO) revealed that the 25 haemonchosis associated regions contain genes involved in multiple biological processes such as immune response and chemotaxis (*MHCII-DRB1, IL20RA, IL17RB, FCER2, HRH1*), response to pain and tissue homeostasis (*RELN, SOX9*), and protein coding, binding, methylation and phosphorylation (*ATP2B1, SOX9, MUC15, UBE2N, LRP8, RELN, NSUN2, LAMC1, ABCB9, PIK3CD, SUGTI, PAK4*) \[[@CR32], [@CR33]\]. Other functions of the identified candidate genes include calcium binding and transport (*LRP8, LAMC1*) and carbohydrate metabolism (*CHI3L2, CHIA*) \[[@CR32], [@CR33]\].

Six of the genes falling within regions of reduced heterozygosity identified in this study (*LRP8*, *ATP2B1*, *LAMC1*, *SOX9*, *MUC15, UBE2N)* were also associated with resistance to *H. contortus* infection in a recent GWAS study using a backcross population of Red Maasai and Dorper sheep under natural infection conditions \[[@CR38]\]. A further six genes (*CHI3L2*, *CHIA*, *DENND2D*, *RELN*, *NSUN2*, and *HRH1)* were among the previously reported top 1% of candidate genes for resistance and susceptibility to gastrointestinal nematodes in divergent populations of Romney and Perendale sheep \[[@CR37]\]. Two of the genes (*IL20RA*, *PIK3CD*) were associated with resistance to experimental challenge with *H. contortus* \[[@CR36]\]. In a gene expression study of deliberately infected Chinese Hu sheep, four genes (*ABCB9*, *SUGT1*, *PAK4*, *FCER2)* were found to contribute to the key immunological responses \[[@CR39]\]. More recently, five of the genes (*AREG*, *KIT*, *IL17RB*, *CXCL12*, *CLCR6)* were also found to be up regulated in *H. contortus* resistant Canarian hair sheep \[[@CR40]\]. Three of the 23 genes (*IL20RA, PIK3CD, RELN*) have also been associated with resistance to other gastrointestinal nematodes such as *Trichostrongyle species, Teladorsagia circumcincta* and other *Nematodirus species* \[[@CR36], [@CR38]\].

A total of 37 regions with reduced heterozygosity contained genes associated with adaptive responses. Some of the genes were involved in immune functions (e.g. *IL12RB2, ALCAM*, *APC2, IL1R, 1IL7*), homeostasis (e.g. *HSPA1A, ATP12A, PDK2, NF1, ABCG2*), melanogenesis/ thermotolerance (*GNAI3, LMLN, PLB1, MITF*) and cellular and digestive metabolism (*GLB1, SUCLG2, TRHDE, OLR1*) \[[@CR31], [@CR68], [@CR69]\]. These genes are plausible candidates for resistance to disease, heat tolerance, or the ability to exist on low quality diets in the harsh, hot and humid climatic conditions faced by these sheep breeds.

Twelve of the 37 low heterozygosity regions (*NPR2*, *ABCG2*, *FGF5*, *MSRB3*, *PDK2*, *NF1*, *NFATC2*, *OR2AG1*, *PRLR*, *ABHD2*, *MITF*, *GLB1*) were also reported in the top 0.1% of candidate genes identified in a previous genome-wide study for signatures of recent selection in 74 different sheep breeds selected from various regions of the world \[[@CR68]\]. Fourteen of the 37 regions (*GNAI3*, *LMLN*, *BMP4*, *TRHDE*, *ALK*, *IL1R1*, *IL7*, *ATP12A*, *PCDH9*, *PLCB1*, *ELF2*, *PGRMC2*, *ALDH1A3*, and *SUCLG2*) were also among the genes recently reported as candidate adaptive genes in indigenous Egyptian sheep and goat breeds \[[@CR31]\].

More recently, in a study of natural local environmental adaptation, ten regions (*IL12RBB2*, *ALCAM*, *SYNJ1*, *KRIT1*, *TSPAN12*, *APC*, *OLR1*, *IFTM21*, and *BATF2*) were among those reported as being important for adaptation in Dall sheep (*Ovis dalli dalli*) \[[@CR69]\]. This approach combined targeted resequencing of a priori identified candidate adaptive genes of immunity and metabolism in domestic sheep (*O. aries*) and bighorn sheep (*Ovis canadensis*) to develop a panel of SNP markers \[[@CR69]\]. As with Djallonke sheep, Dall sheep have undergone many centuries of natural selection with limited human intervention. In contrast to the tropical climatic conditions for Djallonke and Sahelian sheep, the Dall sheep breed evolved under Arctic and sub-Arctic climatic challenges, and hence the common swept regions have direct bearing to only immune functions and not climatic adaptation.

The many shared adaptive signatures of selection between the Djallonke and Sahelian sheep in this study can be attributed to common selection pressures due to their shared environment over several centuries. Historical admixture has been reported in Djallonke sheep populations in different regions of sub-Saharan Africa \[[@CR70], [@CR71]\]. The high number of shared variants (96%) also supports the possibility of migration between the breeds. Introgression from a breed with a high frequency of a homozygous region may reduce heterozygosity in the recipient breed.

Conclusions {#Sec20}
===========

A whole genome analysis of the Djallonke and Sahelian sheep breeds identified over 1 million novel genomic variants. This large number of novel variants suggests that the two sheep breeds represent unique genetic resources, and hence are important for world sheep diversity. The considerable number of breed-specific SNPs identified in Djallonke and Sahelian sheep could aid the sustainable management of each breed. The results also appear to support previous reports of genetic regions associated to trypanotolerance, resistance to *H. contortus* infection and adaptation to a harsh tropical climate. The genomic evidence of trypanotolerance, inferred from conserved orthologues of trypanotolerant Ndama cattle, suggests evidence of similar adaptive selection response for a common disease in two different ruminant species. However, a more comprehensive genetic study in a larger dataset coupled with clinical parasitology will be required to a make any definitive statement.

Supplementary information
=========================

 {#Sec21}

**Additional file 1:** Shows the details of Genomic relationship matrix output for this analysis for all the individual Djallonke and Sahelian sheep sampled for this study.

Animal QTLdb

:   Animal Quantitative Trait Loci

CLR

:   composite of likelihood ratio

DJ

:   Djallonke

GATK

:   Genome Analysis Tool Kit

GCTA

:   Genome-wide Complex Trait Analysis

GRM

:   Genomic relationship matrix

HomSI

:   Homozygosity Stretch Identifier

iHS

:   Integrated haplotype homozygosity score

PCA

:   Principal component analysis

SA

:   Sahelian

SNP

:   Single nucleotide polymorphism

UTR

:   Untranslated regions

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information
=========================

**Supplementary information** accompanies this paper at 10.1186/s12864-019-6198-8.

The authors thank Dr. Abdulrahman Iddriss of University of Development Studies in Ghana for invaluable support during in Ghana. A special thanks to the farm staff and the managers of Ejura and Pong Tamale sheep breeding station in Ghana for facilitating the sampling of blood. We greatly thank CHIRI Bioscience of Curtin University for the use of its facilities.

DMG and KM supervised the project and provided a substantive contribution to all aspects of this work presented, it's conceptualisation, experimental design, supervision of experiments, performance, analysis and interpretation of the data. MJS provided technical advice on the project concept. MY did all the blood sampling and DNA extractions. MY drafted the manuscript under the direct supervision and advice of DMG, KAM and MJS. RJNA performed the initial pre-processing and mapping of sequenced data with Torrent Suite 4.6 and reviewed the associated method section of manuscript. MY, EM, and KAM performed all the GATK pipeline analyses. MY, EM and KAM did the HomSI analysis and the snpEff annotation. MY did the GO analyses. All the genomes were sequenced in RJNA's laboratory. All Authors have read and approved the final manuscripts.

MY was funded by a Curtin University International Strategic Research Scholarship (CIRS). This work was part funded by the BBSRC Animal Health Research Grant BB/L004070/1. Neither the BBSRC nor Curtin University played any role in the design of the study, collection, analysis, interpretation of data, and in writing the manuscript.

The data has been generated from this study has been release to the European variant archive of the European Bioinformatics Institute with the accession number PRJEB15642. The data has also been shared with the International Sheep Genomics Consortium.

The study was performed according to the Australian Code of Practice for Care and Use of Animals for Scientific Purposes. The Curtin University Animal Ethics approval number is AEC_2014_35.

Not applicable.

The authors declare that they have no competing interests.
