African indigenous sheep originated in the Near East. They arrived, in the first instance, in North Africa via the Isthmus of Suez by the seventh millennium before present (BP) (Marshall, 2000). These sheep were of thin-tail type and their dispersion southwards into East Africa followed possibly the Nile river valley and the Red Sea coastline (Blench and MacDonald, 2006; Gifford-Gonzalez and Hanotte, 2011). The second wave brought fat-tail sheep into North and Northeast Africa via two entry points, the Isthmus of Suez and the Horn of Africa across the straits of Bab-el-Mandeb, respectively. The fat-rump sheep are a recent introduction and represent the third wave of arrival and dispersal of the species into eastern Africa (Epstein, 1971; Ryder, 1983; Marshall, 2000).

Sheep fulfill important socio-cultural and economic roles in the Horn of Africa. In Ethiopia they provide a wide range of products, including meat, milk, skin, hair, and manure, and are a form of savings and investment (Assefa et al., 2015). Ethiopia hosts many indigenous breeds of sheep, with currently 14 recognized populations/breeds, which are defined based on their geographic location and/or the ethnic communities managing them (Gizaw, 2008). Based on structure analysis, Edea et al. (2017) showed that the five Ethiopian indigenous sheep populations they analyzed clustered together based on their geographic distribution and tail phenotypes.

Fat depots act as an energy reserve that allows sheep to survive extreme environments and conditions such as prolonged droughts, cold, and food scarcity (Atti et al., 2004; Nejati-Javaremi et al., 2007; Moradi et al., 2012). Based on the combination of tail type and length, Ethiopian indigenous sheep can be classified as short fat-tail, long fat-tail, thin-tail, and fat-rump sheep. The short fat-tail inhabit sub-alpine mountainous regions, the long fat-tail predominate in mid- to high-altitude environments and the fat-rump sheep occur in semi-arid and arid environments (Gizaw et al., 2007). These populations are considered to be adapted to their production environments and they represent an important model species for investigating and enhancing our knowledge on the genome profiles of environmental adaptation, tail morphology, and fat localization.

Different approaches, that contrast groups of fat- and thin-tail sheep, have been used to identify candidate regions and genes associated with tail formation and morphotypes. Moradi et al. (2012) identified three regions on chromosomes 5, 7 and X associated with tail fat deposition in Iranian breeds. Using two fat-tail (Laticauda and Cyprus fat-tail) and 13 Italian thin-tail breeds, Moioli et al. (2015) identified BMP2 and VRTN as the likely candidate genes explaining fat-tail phenotype in the studied populations/breeds. Zhu et al. (2016) detected several copy number variations intersecting genes (PPARA, RXRA, and KLF11) associated with fat deposition in three Chinese native sheep (Large-tail Han, Altay, and Tibetan). Several candidate genes with possible links to fat-tail development, i.e., HOXA11, BMP2, PPP1CC, SP3, SP9, WDR92, PROKR1, and ETAA1, were identified using genome scans that contrasted fat- and thin-tail Chinese sheep (Yuan et al., 2017). Whole genome sequencing of extremely short-tail Chinese sheep revealed the T gene as the best possible candidate, among other nine genes, influencing tail size, following its association with vertebral development (Zhi et al., 2018). There is, so far, no information on the genetic basis of variation in tail fat distribution and size in African indigenous sheep.

In this study, using the Ovine 50 K SNP BeadChip genotypes, we investigated the (i) genetic relationships and structure within and between Ethiopian indigenous sheep of different fat-tail morphotypes alongside other sheep populations and breeds from the Caribbean, European, Middle East, China and Africa, and (ii) candidate genome regions and genes associated with tail morphology, fat deposition and possible eco-climatic adaptation in African indigenous sheep. For the latter, 11 Ethiopian indigenous sheep of different fat-tail morphotypes and two populations of thin-tail sheep from Sudan were analyzed.

In this study, we used Ovine 50 K SNP BeadChip generated genotype data to investigate autosomal genetic diversity in Ethiopian indigenous sheep. Including populations from other regions of the world and the African continent allowed us to assess this diversity in a global geographic context. Our findings showed that the Ethiopian indigenous sheep are genetically differentiated from the other populations including other African fat-tail sheep (Figures 2, 3). The finding that the Ethiopian fat-tail sheep are distinct from those found in North Africa, support the presence of at least two genetic groups of fat-tail sheep in the continent and two separate introduction events, one via the Northeast Africa and the Mediterranean Sea coastline, and the other via the Horn of Africa crossing through the strait of Bab-el-Mandeb, respectively. The distinct clustering of the thin-tail sheep suggests its independent introduction into the continent. The fact that the South African Ronderib and Namaqua sheep occur on the same PC planar axis with Ethiopian sheep (Figure 2) may suggest, a common genetic heritage between the two rather than with the North African breeds. The movement of sheep southwards remains speculative; some linguistic evidence suggests movement of bantu speaking populations from West Africa to South Africa through central Africa and following a western route rather than the more traditionally postulated eastern routes from East to South Africa (Newman, 1995). In such context a close clustering of the thin-tail West African sheep with some fat-tail southern African sheep breeds, such as the Namaqua from Namibia studied here is worth mentioning as it offers some possible insights. This however, will require further investigation beyond the scope of this study.

Our results agree with previous findings that were arrived at using microsatellite loci (Muigai, 2003) and 50 K SNP genotype data (Mwacharo et al., 2017). They are also in line with archaeological and anthropological evidences indicating the introduction first, of thin-tail sheep into the continent followed by fat-tail sheep, initially through the Sinai Peninsula and later the Horn of Africa region (Gifford-Gonzalez and Hanotte, 2011; Muigai and Hanotte, 2013).

Interestingly, the PCA results involving Ethiopian and Sudanese sheep separate the Ethiopian populations into three groups while ADMIXTURE revealed four genetic clusters in Ethiopian sheep irrespective of their geographic origins in the country. TreeMix revealed extensive gene flow between populations of different geographic origins and tail-types. These results suggest, most likely, current and historical intermixing of sheep following human socio-cultural and economic interactions. This appears to be a common feature in Ethiopia and most likely the Northeast and eastern Africa region as it was also observed in Ethiopian goats by Tarekegn et al. (2018). We propose here that the common D genetic background present in short fat-tail and fat-rump sheep may represent historical introgression of the thin-tail gene pool into short fat-tail and fat-rump genepool. This result calls for further investigation.

Our findings on the genetic relationships and differentiation between Ethiopian sheep populations agree with findings of previous studies, which were performed using either microsatellites (Gizaw, 2008) or 50 K SNP genotype data (Edea et al., 2017) and which indicated a grouping of Ethiopian indigenous sheep populations based on their tail phenotypes. However, uniquely in our study, the long fat-tail populations were further subdivided into two secondary groups representing sheep populations from western and southern Ethiopia (Figure 4). These two groups were also defined by different genetic backgrounds by ADMIXTURE (Figure 6) and they clustered separately in TreeMix (Figure 7). In addition, although they are defined by the same tail phenotype, the two populations of Ethiopian short fat-tail sheep did not cluster together. Geographic isolation coupled, most likely, with adaptation to different eco-climates, as well as ethnic, cultural and religious practices and differences, that can impede gene flow, may have shaped this genetic sub-structuring (Madrigal et al., 2001; Gizaw et al., 2007).

In selection signature analysis, we contrasted groups of Ethiopian indigenous sheep that showed variation in the size of the fat-tail with thin-tail sheep. Our results identified several genes as potential candidates controlling tail morphotype and fat localization in the study populations. Several genes occurred within candidate regions that overlapped between at least two of the three approaches used to detect signatures of selection (hapFLK, F_ST, Rsb). The F_ST approach detects signatures arising from an increase or decrease in allele frequency differentiation between populations/breeds, hapFLK detects the same but based on increase/decrease in haplotype frequency differentiation between populations while accounting for hierarchical population structure (Kijas, 2014) while Rsb detects signatures associated with the patterns of linkage disequilibrium between loci across the genome (Oleksyk et al., 2010; de Simoni Gouveia et al., 2014). Since these methods are based on different algorithms and assumptions, if common signatures are detected by at least two of the methods it suggests good reliability of the results while reducing the likelihood of interpreting false positives. They also detect signatures spanning different time periods; the F_ST and hapFLK detect signatures arising from long term differential selection while Rsb detects ongoing signatures of selection including those that arise in the short to medium term (Oleksyk et al., 2010).

In the E1^*S comparison, three genes associated with growth traits were present on the candidate region on Oar2, i.e., histidine triad nucleotide binding protein 2 (HINT2), sperm associated antigen 8 (SPAG8) and natriuretic peptide receptor 2 (NPR2). Previous studies reported these genes to be associated with birth and carcass weights, and fat depth, respectively, in cattle (Casas et al., 2000; McClure et al., 2010) and sheep (Moradi et al., 2012; Wei et al., 2015). We also identified two genes on Oar5 (ANGPTL8, INSR), which might be responsible for fat accumulation in adipose tissues. Angiopoietin-like 8 (ANGPTL8), when induced by insulin receptor (INSR), inhibits lipolysis and controls post-prandial fat storage in white adipose tissue and directs fatty acids to adipose tissue for storage during the fed state (Mysore et al., 2017). The ADAMTS3 (ADAM metallopeptidase with thrombospondin type 1 motif 3) gene was present in the region identified on Oar6. This gene is expressed in cartilage, where collagen II is a major component, as well as in embryonic bone and tendon, suggesting that it could be a major procollagen processing enzyme in musculoskeletal tissues (Dubail and Apte, 2015). The homeobox B13 (HOXB13) and ALX homeobox 4 (ALX4) were identified on the candidate region on Oar11 and Oar15, respectively. Mutations in the former result in overgrowth of caudal spinal cord and tail vertebrae in mice (Economides et al., 2003), while the latter is involved in the development of limbs and skeleton (Fariello et al., 2014).

Our enrichment analysis for the E1^*S genes revealed a cluster of genes (BMP4, MED1) with functions that could possibly be related to tail formation. Bone Morphogenetic Protein 4 (BMP4) was revealed by Rsb and F_ST to be on a candidate region on Oar7 and it has been implicated in tail formation (Moioli et al., 2015). Peroxisome Proliferator Activated Receptor Gamma (PPARG) expression has been associated with back-fat thickness in sheep (Dervish et al., 2011). Ge et al. (2008) reported Mediator Complex Subunit 1 (MED1) is an essential protein for the optimal functioning of PPARG. Despite this association, our analysis did not reveal any signals spanning PPARG, but two of our methods (Rsb and F_ST) revealed a signature on Oar20 that spanned the PPARD gene, a paralogue to PPARG.

In the same comparison (E1^*S), we identified a cluster of genes (CDH8, ADRB3, THRA, TRPM8, PLAC8) that are associated with the GO biological process, response to cold. This is not surprising considering that three out of the four E1 populations are living at a high altitude and therefore in a relatively cold habitat. Indeed, Adreno receptor Beta 3 (ADRB3) plays a major role in energy metabolism and regulation of lipolysis and homeostasis (Wu et al., 2012). It is also associated with birth weight, growth rate, carcass composition and survival in various breeds of sheep (Horrell et al., 2009). The ion channel TRPM8 has been reported to play a major role in eliciting cold defense thermoregulation, metabolic and defense immune responses in humans (Kozyreva and Voronova, 2015).

Several other genes occurring in the E1^*S candidate regions and which are associated with the GO term embryonic skeletal system development (GO:0048706) included HOXC6, SULF2, WNT11, and HOXB9. WNT11 was identified by ZF_ST on Oar15 while HOXC6 and HOXB9 were revealed by hapFLK on Oar3 and Oar13, respectively. The WNT gene family and the T gene have been implicated in vertebral development in laboratory mice (Greco et al., 1996), and with the short-tail phenotype in sheep (Zhi et al., 2018). In addition, the roles of the WNT gene family in lipid metabolic processes in fat-tail sheep have also been reported (Kang et al., 2017). The HOX genes represent transcriptional regulatory proteins that control axial patterning in bilaterians (Garcia-Fernàndez, 2005), where the inactivation of one of the HOX genes often causes transformations in the identity of vertebral elements (Mallo et al., 2010). The HOX genes are able to control morphologies along the anteroposterior axis (Lewis, 1978). Furthermore, HOXC11, HOXC12, and HOXC13 developmental genes were found to be expressed in the tail region indicating their possible associations with tail size and fat development in fat-tail sheep (Kang et al., 2017).

The candidate regions revealed by the E2^*S comparison, spanned 65 candidate genes. Three genes of the BPI fold Containing Family B (BPIFB3, BPIFB4, and BPIFB6) were present in a candidate region on Oar13. These, along with other paralogs (BPIFB1, BPIFA3, BPIFB2, BPIFA1), formed a cluster of functional genes related to the GO term, lipid binding functional process (Supplementary Table 6). In contrast to the E1^*S comparison, the cluster of genes identified in the E2^*S comparison were associated with the GO terms, Magnesium ion binding, response to gamma radiation and cellular response to heat. This suggests, most likely, the propensity of this group of sheep to adapt to the eco-climatic conditions prevailing in their home-tract. This is consistent with the humid highland and moist lowland conditions of the geographic area where the populations representing the E2 group (Bonga, Gesses, Kido) were sampled. High fecundity and prolificacy is a common reproductive trait preferred by farmers in the Bonga sheep (field observations by the last author). This may explain the occurrence of the CIB4 and PRKAA1 in a candidate region in the E2^*S comparison. The CIB4 gene was suggested to be linked, in some way, to high fecundity in the small Tail Han sheep (Yu et al., 2010) and PRKAA1 is involved in ewe's follicular development (Foroughinia et al., 2017).

The third comparison (E3^*S) resulted in 36 genes that occurred in candidate regions that were revealed by at least two methods used to detect selection signatures. Fatty acid binding proteins FABP3 and FABP1 found on candidate regions on Oar2 and Oar3, respectively are the genes that relate most closely to fat deposition. SREBF1 along with PPARG are the main transcription factors controlling lipogenesis in adipose tissue and skeletal muscle (Ferré and Foufelle, 2010), and are mainly regulated by fatty acid-binding proteins (FABP) (Lapsys et al., 2000). Recently, Bahnamiri et al. (2018) evaluated the effects of negative and positive energy balances on the expression pattern of these genes in fat-tail and thin-tail lambs. They observed differential transcriptional regulation of lipogenesis and lipolysis during periods of negative and positive energy balances in the two groups of lambs. In general, the cluster of genes identified in this comparison was significantly enriched for GO terms relating to skin development, wound healing and regulation of actin cytoskeleton reorganization (Supplementary Table 6).

The overlapped genes between all comparisons are shown in Supplementary Figure 7. The commonest genes between the three comparisons are TSPAN8, RXFP2, and RIN2. The TSPAN8 (Tetraspanin 8) occurred in the candidate region on Oar3; it is among the genes that are reported to be associated with insulin release and sensitivity, and obesity in humans (Grarup et al., 2008), while the relaxin family peptide receptor 2 (RXFP2) has been associated with horn morphology (Johnston et al., 2011; Wiedemar and Drögemüller, 2015).

Twelve genes (MELK, RNF38, GNE, CLTA, CCIN, RECK, HINT2, SPAG8, NPR2, FAM221B, MSMP, RGP1) were common between E1^*S and E2^*S comparisons. On Oar2, three genes were identified within the overlapping candidate region, i.e., CLTA which is associated with prion protein deposition in sheep (Filali et al., 2014), GNE which is important for the metabolism of sialated oligosaccharides in bovine milk (Wickramasinghe et al., 2011) and RECK which encodes an inhibitor of angiogenesis, invasion and metastasis, DNA methylation, and increased mRNA in cell lines in humans (Su, 2012). Other genes (i.e., HINT2, SPAG8, and NPR2) are associated with fat deposition in sheep as herein discussed for each of the three comparisons.

Furthermore, one gene (DIS3L2) was in a candidate region that overlapped between the E1^*S and E3^*S comparisons. DIS3 like 3'-5' exoribonuclease 2 (DIS3L2) has also been identified, among genes involved in cancer, cellular function and maintenance, and neurological disease, in a candidate region under selection in cattle (Gautier et al., 2009). In sheep, using F_ST, iHS, and Rsb, de Simoni Gouveia et al. (2017) indicated that DIS3L2 is among genes associated with height variation. In addition, DIS3L2 has reportedly been associated with the Perlman syndrome, which is characterized by overweight in humans (Astuti et al., 2012).

Finally, seventeen genes (PKD2L2, FAM13B, WNT8A, NME5, BRD8, KIF20A, CDC23, GFRA3, CTNNA1, LRRTM2, SIL1, SPATA24, DNAJC18, SMIM33, TMEM173, FRY, ATP10A) were in candidate regions that overlapped between the E2^*S and E3^*S comparisons. Among these, DnaJ heat shock protein family (HSP40) member C18 (DNAJC18) and spermatogenesis associated 24 (SPATA24) on Oar5 were reported among genes involved in heat stress tolerance and male reproductive function, respectively, in East African Shorthorn Zebu cattle (Bahbahani et al., 2015).

Overall, our results revealed four distinct autosomal genomic backgrounds (A, B, C, D) in Ethiopian indigenous sheep. The genotypes of most of the individuals analyzed were made up of at least two genetic backgrounds which could be accounted for by some level of current or historical admixture between populations. Selection signature analysis identified several putative candidate regions spanning genes relating to skeletal structure and morphology, fat deposition and possibly adaptation to environmental selection pressures. Our results indicate that Ethiopian indigenous sheep could be a valuable animal genetic resource that can be used to understand genetic mechanisms associated with body fat metabolism and distribution. This is especially important because fat deposits are a crucial component of adaptive physiology and excessive fat deposition in adipose tissue can result in obesity and overweight, and energy metabolism disorders in humans.

Genotypic data of 160 animals representing eleven Ethiopian and two Sudanese sheep populations are deposited and available at (https://www.animalgenome.org/repository/pub/NOTT2018.0423/).

The animals used in this study are owned by farmers. Prior to sampling, the objectives of the study were explained to them in their local languages so that they could make an informed decision regarding giving consent to sample their animals. Government veterinary, animal welfare and health regulations were observed during sampling. The procedures involving sample collection followed the recommendation of directive 2010/63/EU. Skin tissues importation and/or exportation was permitted by the Ethiopian Ministry of Livestock and Fisheries under Certificate No: 14-160-401-16.

AbA, JM, and OH conceived and designed the study. AbA analyzed the data and together with JM wrote the manuscript. JM and OH revised the manuscript. HB provided support in data analysis. SM, FP, and EC contributed to genotyping and genotype data of non-Ethiopian breeds (Najdi, Omani, and Libyan Barbary) and provided critical inputs on data analysis and in writing the manuscript. FA, MA, and MOA supported the sampling and genotyping of Najdi, Omani and Libyan sheep. AK and AyA lead and coordinated the sampling of Ethiopian sheep HM lead and coordinate the sampling of Sudanese sheep. All authors contributed to the interpretation of the results based on their knowledge on local indigenous sheep genetic resources of their respective countries. All the authors read and approved the final manuscript.