Whole Genome Resequencing Reveals Genetic Diversity and Selection Signatures of Ethiopian Indigenous Cattle Adapted to Local Environments

Abstract

Cattle are among the most important domesticated bovid species in the world, of which Ethiopia possesses large populations adapted to different agro-ecologies and production systems. Though several molecular population genetic studies have been done on Ethiopian indigenous cattle, genomic diversity and selection signatures pertinent to adaptation to the different local environments have yet to be comprehensively characterized. Here, the whole genome sequences of 151 samples from 14 Ethiopian indigenous cattle populations were analyzed to assess genomic diversity and differentiation as well as signatures of positive selection (using Hp, iHS, F_ST, and XP-CLR) in comparison to Sudanese zebu, Asian zebu, Ankole, and African and European taurine cattle. High genomic differentiation was observed between Ethiopian and non-Ethiopian cattle populations, while low genomic differentiation and inbreeding were present between and within Ethiopian cattle populations. Sixteen genome regions overlapping with 40 candidate genes were commonly identified by at least three genome scan methods. High frequencies of missense SNPs in ITPR2, CHADL, GNAS, STING1, and KIT genes with high haplotype differentiations were observed in Ethiopian cattle compared to non-Ethiopian cattle. The candidate genes were significantly associated with several biological functions and molecular pathways responsible for nutrient metabolism, skeletal development, immune response, reproduction, water balance, coat color pigmentation, and circulatory homeostasis. Our results provide new insights into the adaptation of the Ethiopian indigenous cattle to the country’s diverse environments.

1. Introduction

Our knowledge of the origin of cattle in Africa stems from archaeological evidence and DNA analysis. According to the literature, the taurine Bos taurus taurus and indicine Bos taurus indicus cattle arose from the wild aurochs, Bos primigenius [1,2,3]. Taurine cattle were domesticated in the Fertile Crescent, and indicine cattle in the Indus River Valley of the Indian subcontinent [4,5]. Following the past presence of aurochs on the African continent, an African center of cattle domestication has also been proposed, this is, however, highly disputed today [4,6]. The East and the Horn of Africa have witnessed a wave of the introduction and dispersion of indicine cattle into the African continent after the introduction of taurine cattle [6,7,8,9]. The geographical proximity of Ethiopia to the introduction routes of indicine cattle into the African continent has made the country a center of cattle diversification [10].

The African continent possesses 150 unique indigenous cattle populations classified into zebu, taurine, sanga (crosses of zebu and taurine), and zenga (crosses of sanga and zebu) cattle breed groups [11,12], among which zebu, sanga, and zenga dominated the East African region [9]. Ethiopian cattle populations fall into the groups of sanga, zenga, and zebu, with the latter further sub-classified into large and small East African zebu cattle [11,12]. The Sheko, previously classified as taurine cattle, is now re-considered as a taurine × zebu crossbred [10,13]. These diverse cattle genetic resources are adapted to different agro-ecologies where they contribute to major economic activities such as food production, crop cultivation, social welfare, and cultural practices in different cattle production systems [14].

The different agro-ecologies of Ethiopia are home to different locally adapted cattle populations [15]. Over time, in line with the evolution of cattle production systems (pastoralism and mixed crop-livestock farming), genetic introgression, through crossbreeding and migration, has increased the level of admixture that enriched the genetic diversity of the populations with opportunities for rapid adaptation [13,16,17]. Physical and social barriers, including ecological differences, breeding systems, religions, and cultural preferences, have shaped the genetic backgrounds of these admixed populations while restricting gene flow across populations. These factors have all contributed to the unique, locally adapted indigenous cattle genetic resources of the country [18].

Different molecular techniques have been applied to characterize the genetic diversity and the unique features of cattle and other livestock species adapted to different environments. For example, previous reports have used mitochondrial DNA, microsatellite DNA, and Y-chromosome markers to investigate the origin, distribution, and genetic diversity of Ethiopian cattle [8,9,19,20,21], but these neutral markers only capture a subset of variations among individuals and between populations. The development of low to high-density genotyping arrays for single nucleotide polymorphisms (SNPs) further improved the accuracy of detecting genetic diversity within and among populations. As a result, several reports have revealed the diversity and signatures of selection in cattle and other livestock species across their genomes [13,20,21,22]. Still, these genotyping arrays only represent a subset of the genomic diversity. Compared to high-density SNP arrays, whole genome sequences (WGS) allow us to scan all polymorphic loci and increase the accuracy of detecting population genetic structure and diversity as well as the signature of selection across the entire genomes [23,24].

In a previous study, we examined the genetic adaptation of indigenous Ethiopian cattle to high altitudes [25]. We analyzed three Ethiopian highland zebu populations living at very high altitudes (>3000 m above sea level) and three Ethiopian lowland zebu populations. We aimed to identify the specific signatures of selection to high altitude adaptations. In this study, we used WGS data to assess the genomic diversity, population genetic structure, and signatures of selection of 14 Ethiopian indigenous cattle populations distributed in different agro-ecologies of the country. We aim to comprehensively report all adaptations in the Ethiopian cattle populations (e.g., climatic and disease challenges). The WGS of Ethiopian indigenous cattle was analyzed in comparison with non-Ethiopian cattle, including Sudanese zebu, Asian zebu, Ankole, and African and European taurine cattle. These cattle belong to different lineages and breed groups adapted to different environments and production systems. The comparative analysis was designed to capture the signatures of positive selection that contributed to the diversity and adaptation of Ethiopian indigenous cattle populations to the environments of the country. Our study is the most comprehensive so far on the genomic diversity and adaptation of Ethiopian indigenous cattle.

2. Materials and Methods

2.1. Cattle Populations and Sample Collection

Fourteen Ethiopian indigenous cattle populations (ETZ) living at different altitudes and agro-ecologies were included in this study (Bale, Choke, Semien—cold, humid, high altitude; Arsi, Horro, Fogera—cold, humid, mid altitude; Sheko—hot, humid, mid altitude; Mursi, Goffa, Bagaria, Begait, Boran—hot, humid, low altitude; Afar and Ogaden—hot, arid, low altitude) (Figure 1 and

Table 1. Ethiopian cattle populations, regions of sampling locations, altitudes, GPS coordinates, and agro-ecologies.

No.	Cattle Population	Region (District)	Altitude (m.a.s.l)	GPS Coordinates		Climate/Agro-Ecology
				Latitude	Longitude
1	Afar	Afar (Werer/Asayta)	804	9.35	40.17	Hot, arid, low altitude
2	Arsi	Arsi (Arsi Bekoji)	2800	7.58	39.28	Cold, humid, mid altitude
3	Bagaria	Benishangul (AL mahal)	680	11.55	35.16	Hot, humid, low altitude
4	Bale	Bale (Bale Mountain)	3586	6.77	39.76	Cold, humid, high altitude
5	Begait	Gonder (Humera)	895	14.10	37.22	Hot, humid, low altitude
6	Boran	Borena (Borena plain)	1368	4.55	38.10	Hot, humid, low altitude
7	Choke	Gojam (Choke Mountain)	3410	10.61	37.84	Cold, humid, high altitude
8	Fogera	Bahir Dar (Fogera plain)	1735	11.30	37.28	Cold, humid, mid altitude
9	Goffa	Goffa (Arba Minch)	1100	6.21	37.07	Hot, humid, low altitude
10	Horro	Welega (Horro)	1722	9.07	37.03	Cold, humid, mid altitude
11	Mursi	South Omo (Mursi)	1405	5.47	36.34	Hot, humid, low altitude
12	Ogaden	Somali (Kebri-Beyah)	1200	9.59	41.86	Hot, arid, low altitude
13	Semien	Gonder (Semien Mountain)	3732	13.24	38.14	Cold, humid, high altitude
14	Sheko	Kaffa Shaka	2240	7.65	35.51	Hot, humid, mid altitude

). We compared the genomes of these cattle with non-Ethiopian cattle populations, including the Sudanese zebu (Kenana and Butana, SUZ), Ankole (ANK, sanga), Gir (Asian zebu, ASZ), Angus and Holstein (European taurine, EUT), and Muturu and N’dama (African taurine, AFT). Blood samples from unrelated individuals were collected after discussion with local government authorities and cattle owners. The genomic DNA was extracted from whole blood using the QIAGEN DNeasy Blood & Tissue Kit (https://www.qiagen.com/us/, accessed on 10 March 2018). The integrity of the DNA was checked in a 1% agarose gel. The concentration and purity of the DNA were determined by spectrophotometer reading at 260 nm and 280 nm, respectively (DeNovix Inc., Wilmington, DE, USA).

2.2. Whole Genome Sequences, Reads Mapping, and SNPs Calling

Ethiopian cattle sequences were retrieved from the NCBI SRA database (https://www.ncbi.nlm.nih.gov/sra, accessed on 9 December 2018) with the bio-project accession numbers PRJNA574857 (Afar, Arsi, Begait, Boran, Fogera, Goffa, Horro, Mursi, and Sheko), PRJNA698721 (Bagaria, Bale, and Semien), PRJNA841948 (Choke), and PRJNA312138 (Ogaden). Non-Ethiopian cattle genomes were retrieved from the bio-project accession numbers PRJNA312138 (Ankole, Kenana, and N’Dama), PRJNA176557 (Angus and Holstein), PRJNA386202 (Muturu), and PRJNA574857 (Butana) (Supplementary Table S1).

The genomes (150 bp of paired end reads) were sequenced using the Illumina platform (Illumina, San Diego, CA, USA) (Table S1). Illumina adapters and low-quality bases from each raw sequence read were trimmed out, and the pair reads were mapped to the ARS-UCD1.2 taurine cattle reference genome [26] using the Burrows–Wheeler Aligner (bwa v0.7.17) [27]. After indexing the alignments and removing duplicates, base quality score recalibration and haplotype caller analyses were performed using the GATK v3.8-1-0-gf15c1c3ef according to the GATK best practices pipeline [28]. For base quality recalibration, the known variants were masked by the known sites listed in the ARS1.2PlusY_BQSR_v3.vcf.gz file provided by the 1000 Bull Genomes project. SNP variants extracted from each sample were combined and recalibrated to the 99.9 truth sensitivity filtering method to remove low-quality SNPs. After quality check and filtering, 51.9 million biallelic autosomal SNPs were detected and used for further analyses.

2.3. Population Genetic Structure

Principal component (PCA) and population admixture analyses were performed to detect the genetic differentiation and structure of all ETZ and non-Ethiopian cattle populations. We pruned the SNPs with LD (r²) ≥ 0.5 in 50-SNPs sliding windows with 10-SNPs step size and removed SNPs with minor allele frequency (MAF) < 0.05 and missing call rates greater than 10% from the dataset. After rigorous filtering, a total of 5.1 million non-linked SNPs were used for admixture analysis. The ADMIXTURE v1.3.0 software [29] was employed to analyze the ancestry proportion of individual cattle genomes (K = 1 to 10). The ancestry proportions were plotted and visualized using the R package. Similarly, the PCA was done after filtering out SNPs with MAF < 0.05 and a genotype missing rate of less than 90%. A total of 25.6 million SNPs was retained after the filtering and the PCA were calculated for two datasets using PLINK v1.9 [30]. The first dataset included all cattle populations (ETZ and non-Ethiopian cattle populations), while the second dataset contained only the 14 ETZ populations, to detect the genetic structure of ETZ in comparison to the other cattle population and within ETZ populations, respectively. The PC1 and PC2 eigenvector values for the two datasets were plotted with the ggplot2 in the R package.

2.4. Genomic Diversity

2.4.1. Nucleotide Diversity, Population Genetic Differentiation, and Heterozygosity

Genomic diversity was measured using nucleotide diversity (π), population genetic differentiation (F_ST), and genomic heterozygosity. The π and F_ST values were calculated within 100 kb windows with 50 kb step size along the autosomes using the vcftools v0.1.15 [31]. Genetic differentiation within and among ETZ and non-Ethiopian cattle populations was estimated using pairwise F_ST [32]. The mean weighted-F_ST values of individual windows were taken for pairwise comparison between populations. The observed heterozygosity was calculated as the proportion of total heterozygous SNPs to the total number of sites counted in each genome, based on the ARS-UCD1.2 taurine cattle reference genome. The observed and expected heterozygosity were calculated using the –het option of the PLINK v1.9 [30] for each genome and then averaged for each population.

2.4.2. Runs of Homozygosity-Based Genomic Inbreeding Coefficient (F_roh)

The genomic inbreeding coefficient for each animal was calculated from the number and the size of haplotype autozygosity along the genomic regions called runs of homozygosity (ROHs). ROHs are uninterrupted stretches of homozygous genotypes common among individuals within a population that enable reliable estimation of the level of inbreeding [33,34,35,36]. ROHs were estimated from ped and map plink files using the slidingRUNS.run function in the detectRUNS R-package [35]. The number of ROHs, ROH length, and genomic inbreeding coefficient (F_roh) per individual genome were summarized using the summaryRuns function. The algorithm to detect ROHs relies on a scanning window approach, and the analysis was done for a stretched homozygous haplotype length of ≥0.5 and ≥1 Mb genomic size (

Table 2. Length and number of ROHs, genomic inbreeding coefficient (F_roh), nucleotide diversity, observed and expected heterozygosity in Ethiopian and non-Ethiopian cattle populations.

		ROHs ≥ 0.5 Mb		ROHs ≥ 1.0 Mb		Nucleotide Diversity
Breed	N	No.	Froh	No.	Froh		E(Het)	O(Het)
Afar	14	416	0.0132	43	0.0024	3.49 × 10−3	0.279	0.273
Arsi	10	731	0.0224	107	0.0059	3.44 × 10−3	0.279	0.259
Bagaria	10	553	0.0157	57	0.0029	3.59 × 10−3	0.280	0.286
Bale	10	270	0.008	17	0.0014	3.61× 10−3	0.280	0.291
Begait	9	468	0.0197	55	0.0048	3.45 × 10−3	0.279	0.264
Boran	10	460	0.0124	32	0.0017	3.52 × 10−3	0.279	0.260
Choke	10	216	0.0069	7	0.0005	3.60 × 10−3	0.280	0.289
Fogera	12	695	0.0247	128	0.008	3.15 × 10−3	0.278	0.241
Goffa	13	462	0.0137	66	0.0035	3.11 × 10−3	0.279	0.242
Horro	11	681	0.0184	84	0.0051	3.38 × 10−3	0.279	0.254
Mursi	10	616	0.0202	70	0.0043	3.32 × 10−3	0.279	0.252
Ogaden	9	663	0.0316	142	0.0119	2.88 × 10−3	0.275	0.193
Semien	10	167	0.0054	6	0.0008	3.64 × 10−3	0.280	0.294
Sheko	13	622	0.0204	87	0.005	2.67 × 10−3	0.278	0.217
Butana	10	628	0.0225	129	0.0079	3.30 × 10−3	0.279	0.264
Kenana	10	1008	0.0347	248	0.0145	3.06 × 10−3	0.277	0.214
Ankole	10	763	0.0229	86	0.0057	2.39 × 10−3	0.276	0.177
Muturu	10	3864	0.119	596	0.0319	1.26 × 10−3	0.278	0.075
N’Dama	10	2200	0.0931	654	0.0525	2.26 × 10−3	0.277	0.125
Gir	9	1087	0.0355	220	0.0126	2.41 × 10−3	0.275	0.178
Angus	10	2762	0.129	1112	0.085	1.39 × 10−3	0.279	0.121
Holstein	11	2325	0.1024	847	0.0627	1.43 × 10−3	0.279	0.130

Note: ROHs = Runs of homozygosity; No. = Number of ROHs of sequence length > 0.5 Mb and >1 Mb; F_roh = genomic inbreeding based on ROHs; E(HET) = Expected heterozygosity; and O(HET) = Observed heterozygosity.

). Therefore, the genomic inbreeding coefficient of an animal based on the ROH length was calculated as: $F_{roh}=\sum L_{ROH}/L_{auto}$ [36], where ∑L_ROH is the sum of the length of all ROHs detected in an individual genome, and L_auto is the total length of all autosomes.

2.4.3. Linkage Disequilibrium (LD) and Effective Population Size

The squared correlation coefficient between two SNPs (r²) was used to measure the LD value of all paired SNPs. The pairwise LD between SNPs was calculated using the PopLDdecay v3.40 directly from the vcf file [37]. The r² was detected within a minimum of 20 SNPs apart and a maximum distance of 1 Mb [38]. Subsequently, the LD decay between SNP pairs was graphically plotted based on the average r² values for each 20 kb using the R package script in the PopLDdecay [37]. The effective population sizes (Ne) and their trends across generations were calculated for each cattle group using the SNeP v1.1 [38]. The SNeP estimates historical effective population sizes using the following equation ${E(r}^2)={(1+4N\mathrm{e}c)}^{-1}$ [39], where Ne is the effective population size, and c is the recombination rate in cM with 1 cM equal to 1 Mb [38].

2.4.4. Phylogenetic Relationships among Cattle Samples

Pairwise neighboring joining genetic distances for the 205 cattle samples were calculated after filtering SNPs having a minimum MAF of 0.1, no missing genotype, and pairwise r² between SNPs less than 0.05. We converted the SNPs in vcf format into nucleotide sequences in FASTA format using the vcf2phylip v1.5 (https://github.com/edgardomortiz/vcf2phylip/tree/v1.5). The nucleotide sequences of each sample in the FASTA format were aligned using the MEGA v11 [40]. The phylogenetic tree was reconstructed after aligning a total of 3.6 million SNPs from all 29 autosomes, and the evolutionary relationship was inferred using the neighbor-joining method [41]. The evolutionary distance was computed using the maximum composite likelihood nucleotide substitution model [42]. All ambiguous positions were removed for each nucleotide sequence using the pairwise deletion, and a phylogeny was built with 1000 bootstrap replicates.

2.5. Detection of Selection Signatures and Their Functional Annotations

Analysis of selection signatures was performed within a population using pooled heterozygosity (Hp) and integrated haplotype homozygosity (iHS). Among populations, it was performed using population differentiation (F_ST) and cross-population composite likelihood ratio (XP-CLR) tests. We used a 100-kb window with a step size of a 50-kb overlapping window. The iHS selection scan was performed using the program rehh of the R package after phasing and imputing the genotype allele for each autosome using the Beagle 5.1 [43]. It is a standardized ratio measuring the length of extended haplotype homozygosity at a given SNP, comparing the ancestral allele to the derived allele [44]. The XP-CLR signature of selection was done by calculating allele frequency differences in two populations. It detects alleles differentially selected in the genomic region of the target population compared to neutral alleles in the reference population [45].

The top significant regions for each selection scan test were run on the Ensembl Biomart online annotation tool (http://www.ensembl.org/index.html, accessed on 14 June 2020), using the ARS-UCD1.2 taurine cattle reference genome [26], after combining overlapping significant 0.5% windows. The Database for Annotation, Visualization, and Integrated Discovery (DAVID, v6.8, https://david.ncifcrf.gov/home.jsp, accessed on 19 June 2020) [46] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway [47] were used to annotate the candidate genes and their functions as well as the function of their products (gene ontology (GO) analysis). To illustrate the difference between candidate gene regions, we drew line plots of nucleotide diversity and F_ST values. To detect missense SNPs at candidate genes, we used the Variant Effect Predictor (VEP 104.3) for the functional annotation of the SNPs with MAF > 0.05.

3. Results

3.1. Population Genetic Structure in Ethiopian Cattle

The population genetic structure and admixture in ETZ populations were assessed using the PCA and admixture analysis. In the first dataset (Figure 2a), PC1 and PC2 explained 63.1% of the total variation among all samples. PC1 (53.3%) separated the zebu cattle (ETZ, SUZ, and ASZ) from the taurine cattle (EUT and AFT) and the admixed ANK cattle. PC2 (9.7%) differentiated the African cattle from non-African cattle (ASZ and EUT). In the second dataset, PC1 and PC2 accounted for 26.6% of the total variation. It separated ETZ populations into three clusters (Figure 2b). Among them, Sheko and Mursi were separately clustered from the other ETZ populations. These two populations are located in the Southwest part of the country, where trypanosome is endemic. Similarly, Begait and Bagaria clustered together, and they were separated from the other ETZ populations. These two populations are adapted to the hot and humid lowland climate in the North and Northwest parts of Ethiopia that border Sudan.

The admixture analysis of the first dataset captured the maximum ancestry proportion at K = 4, which was supported by the lowest coefficient of error variance (Figure 3a). The admixture plot showed that ETZ and SUZ shared similar ancestral proportions as compared to EUT, AFT, and ASZ populations (Figure 3b). The ancestry was derived primarily from indicine cattle B. t. indicus (91.9%), including 78.8% specific to African zebu (ETZ and SUZ) and 13.1% shared with ASZ, while AFT and EUT shared only 6.1% and 2.0% of the ancestry, respectively (Figure 3b). AmongETZ cattle populations, Bale and Boran showed the highest zebu ancestry (99.7%), while Sheko had the highest AFT ancestry (19.5%), followed by Mursi. Sheko displayed a pattern of ancestries similar to Ankole.

A K = 5, the Ankole cattle (a sanga from Uganda) was separated from the Sudanese and Ethiopian cattle. Also, K = 5 and above separated the Mursi and Sheko from the other East African cattle populations. Among the ETZ cattle populations, Bale, Choke, Semien, and Boran showed similar admixture patterns (Figure 3b).

3.2. Genomic Diversity in Ethiopian Cattle Populations

3.2.1. Nucleotide Diversity, Heterozygosity, and Population Genetic Differentiation

Compared to the non-Ethiopian cattle populations, the highest average nucleotide diversities (π) were observed in ETZ populations, ranging from 3.64 × 10⁻³ in Semien to 2.67 × 10⁻³ in Sheko. SUZ showed similar π estimates with ETZ, but Muturu had the lowest π value (1.26 × 10⁻³) among all cattle populations (Figure 4 and Table 2). Similarly, higher observed heterozygosity was observed in ETZ compared to AFT, EUT, and ASZ (Figure 4 and Table 2), e.g., Semien (0.294), Bale (0.291), and Choke (0.289) followed by Bagaria (0.286) and Afar (0.273), while the lowest observed heterozygosities were detected in Muturu (0.075) and Angus (0.121). The high observed heterozygosities in ETZ populations were substantiated by their high π values with a correlation coefficient of 0.96 (p < 0.0001) (Figure 4).

The overall population genetic differentiation (F_ST) values between ETZ populations were lower than the values obtained between non-Ethiopian cattle populations (

Table 3. Genetic differentiation (F_ST) between Ethiopian and non-Ethiopian cattle populations.

Breed	ANG	AFR	ANK	ARS	BAG	BAL	BEG	BOR	BUT	CHO	FOG	GIR	GOF	HOL	HOR	KEN	MUR	MUT	NDA	OGD	SEM
ANG
AFR	0.308
ANK	0.248	0.092
ARS	0.306	0.010	0.075
BAG	0.323	0.019	0.063	0.019
BAL	0.312	0.012	0.086	0.006	0.017
BEG	0.282	0.019	0.091	0.018	0.013	0.021
BOR	0.289	0.010	0.087	0.008	0.019	0.008	0.021
BUT	0.329	0.049	0.119	0.046	0.045	0.053	0.038	0.049
CHO	0.303	0.009	0.079	0.003	0.014	0.003	0.016	0.008	0.048
FOG	0.293	0.009	0.072	0.002	0.017	0.009	0.018	0.013	0.045	0.004
GIR	0.439	0.009	0.205	0.086	0.087	0.089	0.084	0.083	0.112	0.087	0.087
GOF	0.289	0.009	0.070	0.006	0.024	0.013	0.025	0.013	0.049	0.009	0.009	0.095
HOL	0.122	0.009	0.250	0.308	0.325	0.314	0.284	0.291	0.331	0.304	0.296	0.442	0.292
HOR	0.292	0.009	0.068	0.001	0.020	0.009	0.020	0.012	0.047	0.004	0.002	0.092	0.007	0.294
KEN	0.297	0.009	0.084	0.022	0.022	0.031	0.015	0.031	0.046	0.026	0.021	0.095	0.023	0.299	0.022
MUR	0.277	0.009	0.056	0.009	0.029	0.018	0.029	0.018	0.056	0.014	0.011	0.112	0.008	0.279	0.008	0.029
MUT	0.320	0.009	0.268	0.332	0.346	0.335	0.299	0.308	0.364	0.327	0.317	0.479	0.311	0.309	0.315	0.336	0.297
NDA	0.264	0.009	0.189	0.245	0.258	0.252	0.229	0.238	0.234	0.244	0.235	0.387	0.231	0.262	0.231	0.092	0.212	0.234
OGA	0.312	0.009	0.073	0.003	0.029	0.020	0.026	0.015	0.054	0.016	0.012	0.095	0.013	0.317	0.012	0.023	0.007	0.348	0.252
SEM	0.307	0.009	0.084	0.006	0.013	0.005	0.016	0.010	0.047	0.001	0.005	0.086	0.012	0.308	0.007	0.027	0.018	0.329	0.247	0.020
SHE	0.251	0.009	0.048	0.027	0.047	0.034	0.044	0.039	0.074	0.029	0.029	0.139	0.023	0.252	0.023	0.044	0.021	0.267	0.060	0.034	0.032

Note: ANG = Angus, AFR = Afar, ANK = Ankole, ARS =Arsi, BAG = Bagaria, BAL = Bale, BEG = Begait, BOR = Boran, BUT = Butana, CHO = Choke, FOG = Fogera, GIR = Gir, GOF = Goffa, HOL= Holstein, HOR = Horro, KEN = Kenana, MUR = Mursi, MUT = Muturu, NDA = N’Dama, OGD = Ogaden, SEM = Semien, and SHE = Sheko.

). Nevertheless, high F_ST values were detected between Mursi and Begait (0.029), Mursi and Bagaria (0.029), Goffa and Begait (0.025), and Goffa and Bagaria (0.024), while the lowest F_ST values were observed between Arsi and Horro (0.001), Semien and Choke (0.001), and Arsi and Fogera (0.002). In addition, cattle populations from the cold and humid high-altitude environment (Bale, Choke, Semien) showed low genetic differentiation (Table 3). Higher genetic differentiations were observed between ETZ and non-Ethiopian taurine breeds with F_ST values ranging from 0.346 (Muturu and Bagaria) to 0.291 (Holstein and Boran). Between the non-Ethiopian cattle populations, the highest differentiation was observed between Angus and Gir (0.476) (Table 3).

3.2.2. ROHs-Based Genomic Inbreeding Coefficient (F_roh)

The F_roh was calculated from homozygous sequences across each genome at ROHs ≥ 0.5 Mb and ≥1 Mb (Figure 5 and Table 2). We detected fewer ROHs ≥ 1 Mb in ETZ genomes compared to the taurine genomes (Table 2). For ROHs ≥ 1 Mb, the mean F_roh in ETZ populations ranged between 0.0005 (Choke) and 0.0119 (Ogaden) with seven and 142 ROHs, respectively (Table 2), while a relatively large number of ROHs and higher F_roh values were detected in Angus (ROHs = 1112 and F_roh = 0.085) and Holstein (ROHs = 847 and F_roh = 0.0627). For ROHs ≥ 0.5 Mb, the number of ROHs and F_roh values were similarly low across ETZ populations, ranging from F_roh = 0.0054 (Semien) to 0.0316 (Ogaden). The F_roh at ROHs ≥ 0.5 Mb were also higher in non-Ethiopian cattle populations, e.g., F_roh = 0.129 in Angus and 0.112 in Muturu (Figure 5 and Table 2).

3.2.3. Linkage Disequilibrium (LD) and Effective Population Size (Ne)

The LD and Ne were calculated for each cattle group (ETZ, ASZ, SUZ, EUT, AFT, and ANK). The lowest LD value but the highest LD decay were observed in ETZ (Figure 6a). Relatively, the taurine cattle showed the highest LD (r² = 0.7), supporting the presence of long homozygous regions. The Ne of the pooled samples for each cattle group indicated that ETZ, SUZ, and ANK had similar Ne, higher to the one calculated for AFT and closer to the ones of ASZ and EUT, nearly 1000 generations ago (Figure 6b and Table S2). A relatively larger Ne was observed in ETZ (n = 649) compared to non-Ethiopian cattle populations. Around 50 generations ago, it ranged from 151 in AFT to 373 in ANK (Table S2). The relative large Ne in ETZ populations were supported by their high nucleotide diversity, small inbreeding coefficient, and high observed heterozygosity (Figure 4 and Figure 5, Table 2).

3.2.4. Phylogenetic Relationships between Ethiopian and Non-Ethiopian Cattle

A neighbor-joining phylogenetic tree was constructed to assess the genetic relationships among all 205 cattle samples (Figure 6c). All ETZ separately branched from all taurine (African and European), ANK, ASZ, and SUZ cattle. Begait and Bagaria were close to each other, which were located next to the clade of SUZ (Butana and Kenana), while all these four populations were closely related to ASZ. Choke, Semien, Bale, and Boran were also close to each other and they constituted a distinct clade. In this analysis, the taurine and admixed Ankole cattle were clearly separated from the other cattle populations. This phylogenetic topology confirmed the similar genomic differentiation observed previously in the PCA and the admixture analyses.

3.3. Selection Signatures in Ethiopian Cattle

The signatures of selection in ETZ were scanned using methods that can detect selection signatures within a population (Hp and iHS) and between populations (F_ST and XP-CLR). We merged all ETZ samples and ran the pooled heterozygosity analysis in a 100-kb window with an overlapping step size of 50 kb. We detected 154 genomic regions at a ZHp = −3.5 threshold (Figure 7b, Table S3). These genomic regions were found to overlap with 182 protein-coding genes in the Ensemble database (Table S3). The iHS analysis detected 258 candidate genomic regions above the iHS threshold value of 5.0 (p-value of 9.9 × 10⁻³), corresponding to the top 0.5% iHS scores (Figure 7d, Table S4). These candidate genomic regions accounted for a total of 44.45 Mb genomic size and overlapped with 434 protein-coding genes (Table S4). By comparing the pooled ETZ with all non-Ethiopian cattle at ZF_ST = 4 (Figure 7a, Table S5) and at an XP-CLR score = 30 (Figure 7c, Table S6), 152 and 236 candidate genomic regions were identified with 258 and 428 protein-coding genes, respectively.

3.3.1. Comparative Analysis of Selection Sweeps in Ethiopian Cattle

We selected genome regions commonly identified by at least three of the four genome scan methods. This approach detected 16 genomic regions on five autosomes (BTA5, 7, 13, 16, and 19) (

Table 4. Candidate gene regions (Mb) and protein-coding genes in Ethiopian cattle commonly detected by three of the four genome scan methods.

BTA	Region Start	Region Stop	Genes Name	FST	XP-CLR	ZHp	iHS	References
5	48.65	48.75	WIF1	5.78	93	−3.65	-	[48]
	58.55	58.65	OR6C75, OR6C1Q	5.10	31	-	5.86	[49,50]
	112.35	112.45	EP300, L3MBTL2, ZC3H7B, CHADL, RANGAP1	4.29	35	-	5.70	[50,51]
7	50.00	50.10	ENSBTAG00000004415	4.53	102	−3.96	-	[17,52,53]
	50.15	50.25	CTNNA1, LRRTM2	5.21	167	−4.54	-	[17,54]
	50.25	50.35	SIL1	5.77	186	−3.56	-	[17,55]
	50.60	50.70	MATR3, PAIP2	5.69	80	−4.39	-	[17,45,56]
	50.65	50.75	SLC23A1, MZB1, PROB1, SPATA24, DNAJC18, ECSCR, SMIM33, STING1	6.02	287	−4.39	-	[17,53,57]
	50.85	50.95	UBE2D2, CXXC5	6.22	177	−4.70	-	[17,45]
	51.0	51.1	PSD2, NRG2	5.52	336	−4.34	-	[17]
13	57.4	57.5	ENSBTAG00000017475, GNAS	4.11	37	−4.05	5.09	[17,58]
			NELFCD	4.11	37	-	5.09
16	19.5	19.6	USH2A	-	56	−4.09	6.28
19	39.65	39.75	FBXL20	4.48	40	−3.89	-	[17]
	39.75	39.85	CDK12, MED1	3.97	125	−3.89	-	[17,59]
	40.4	40.5	CASC3, RAPGEFL1, WIPF2	5.50	148	−3.63	6.37	[56,60]
			MSL1	5.50	-	−3.63	6.37
	40.3	40.4	THRA, MED24, NR1D1	4.54	43	-	6.37	[61]

). Annotation of these regions with the online Ensembl genome browser revealed 40 protein-coding genes (Table 4). Among them, five candidate genes (GNAS, CASC3, RAPGEFL1, ENSBTAG00000017475, and WIPF2) were detected by all four genome scan methods.

3.3.2. Annotation of Candidate Selection Signatures

We combined the Ensembl id of candidate genes detected by all four genome scan methods and ran them against the online DAVID bioinformatics tool to detect their biological processes, molecular functions, cellular components, and KEGG pathways. Candidate genes with highly similar functions (p ≤ 0.05 and fold enrichment ≥ 1.2) were selected as GO groups relevant to adaptation. We identified gene clusters regulating immunity response, lipid metabolism, and reproduction (

Table 5. GO cluster annotation of candidate genes commonly detected by the Hp, F_ST, iHS, and XP-CLR scan methods.

Term	Count	p Value	Fold Enrichment	Genes
GO:0001501~skeletal system development	46	1.2× 10−7	2.4	FGF18, KIT, RAB33B, HOXC10, HOXC4, HOXD13, HOXD12, HOXD3, IRX5, HOXB7, HOXC9, BMPR2, PDGFRA, HOXC6, ZFAND5, HOXB2, GSC, THRA, FGR, CYP26B1, RARA, ADAMTS12, MED1, MPIG6B, EDNRA, BBS2, HOXD10, HOXD9, EP300, EFEMP1, TIFAB, HOXB3, ACTN3, LHX1, HOXD4, HOXB4, CHADL, HMGA2, HOXB5, HOXD8, HOXB6, GNAS, LY6G6D, HOXB8, WNT5B, HOXB9
GO:0002699~positive regulation of immune effector process	16	2.1 × 10−2	1.9	IL18, IL18R1, KIT, IL23A, SPHK2, STAT5B, FGR, RARA, PRKCZ, TBX21, IL13, IL4, EXOSC3, IL18RAP, MZB1, TRIM6
GO:0002292~T cell differentiation involved in immune response	9	5.1 × 10−3	3.3	IL18, IL18R1, IL23A, RARA, PRKCZ, TBX21, RORA, IL4, LOXL3
GO:0050778~positive regulation of immune response	38	6.8 × 10−3	1.6	IL18, MYD88, IL18R1, STING1, HMCN2, KIT, IL23A, POLR3B, IRAK3, SPHK2, STAT5B, NR1H3, RTN4, FGR, NR1D1, RARA, MED1, PRKCZ, TBX21, CCR7, IL13, IL4, EXOSC3, BAG6, OTUD4, IL18RAP, ALPK1, CHADL, PAWR, CRKL, ELANE, TNIP3, CFD, TRIM6
GO:0042092~type 2 immune response	6	8.4 × 10−3	4.6	IL18, TRAF3IP2, PRKCZ, TBX21, IL4
GO:0019915~lipid storage	9	9.3 × 10−3	3.0	SOAT1, STAT5A, STAT5B, NR1H3, CD36, FAM71F2, PLIN5, EHD1
GO:0071218~cellular response to misfolded protein	5	1.9 × 10−2	4.7	DNAJC18, RNF126, BAG6, AUP1, RNF5
GO:1903706~regulation of hemopoiesis	22	5.0 × 10−2	1.5	IL18, FSTL3, MAFB, IL23A, STAT5A, STAT5B, NUDT21, HSPA9, CYP26B1, RARA, LOX, MED1, PRKCZ, TBX21, AGER, TOB2, CASP8, YTHDF2, IL4, LOXL3, CSF3, HOXB8
GO:0005254~chloride channel activity	9	1.8 × 10−2	2.7	PACC1, ANO6, CLCA1, GABRB3, CLIC1, CLCA4, GLRB, ANO3, CLCA2
bta04915:Estrogen signaling pathway	18	3.5 × 10−4	2.7	TGFA, ITPR2, KRT20, ATF6B, PRKACB, RARA, KRT12, SHC2, KRT40, KRT39, KRT10, SHC4, KRT25, KRT26, KRT27, KRT28, GRM1, GNAS
bta04975:Fat digestion and absorption	9	3.3× 10−3	3.6	PLPP2, AGPAT1, PLA2G2A, CD36, PLA2G5
bta04972:Pancreatic secretion	13	3.7 × 10−3	2.6	ITPR2, PLA2G2A, CLCA1, ATP2B1, CCK, CLCA4, CLCA2, PLA2G5, GNAS
bta04270:Vascular smooth muscle contraction	13	3.9 × 10−2	1.9	ITPR2, PLA2G2A, PRKACB, EDNRA, PLA2G5, GUCY1A2, GNAS

and Table S7). These clusters included T cell differentiation involved in immune response (GO:0002292), positive regulation of immune response (GO:0050778), and type 2 immune response (GO:0042092), which played significant roles in adaptive immunity and homeostasis in response to environmental stresses and diseases [62,63]. The other GO clusters of skeletal system development (GO:0001501) and muscle structure development (GO:0061061) were related to morphological body frame structure and meat production [64]. Candidate genes clustered into the cellular response to misfolded proteins (GO:0071218) and regulation of hemopoiesis (GO:1903706) are thought to be particularly important for environmental stress regulation [65,66].

We also identified KEGG molecular pathways of candidate genes linked to fat digestion and absorption (bta04975), pancreatic secretion (bta04972), and vascular smooth muscle contraction (bta04270) (p < 0.001). They played significant roles in adaptation to heat stress, water imbalance, hypoxia, and nutrient metabolism [67,68,69]. In addition, candidate genes belonging to the estrogen signaling pathway (bta04915) are associated with reproduction, skin development, animal body size, milk yield, and carcass characteristics [58,70].

We further selected candidate genes commonly detected by at least three of the four genome scan methods (Table 4) and clustered them into different functional groups (Table 5). We shortlisted eight candidate genes (CHADL, CTNNA1, MZB1, GNAS, ENSBTAG00000017475, MED1, THRA, and NR1D1). ENSBTAG00000017475 is an uncharacterized candidate gene which overlaps with the GNAS complex locus on BTA13 (57.46–57.47 Mb). Besides these candidate genes, we also selected six genes (IL18, PRKCZ, RARA, KIT, PLA2G2A, and ITPR2) included in many GO clusters and KEGG pathways related to immune response, ion homeostasis, reproduction, and coat color pigmentation. These genes have functions believed to be critical for the adaptation to tropical environments [71,72].

SNP functional annotations were carried out using the variant effect predictor of the Ensembl genome database (VEP), which identified SNPs with missense mutations in ITPR2, CHADL, KIT, STING1, and GNAS (

Table 6. Examples of missense SNPs and their frequencies in Ethiopian cattle populations.

BTA	SNP Id	Position and Type of Allele Substitution	Gene	Amino Acid Change and Position	Exon	Codon Change	Frequency
5	rs482168794	g.83460939T>C	ITPR2	P: M2195T	49/57	aTg/aCg	0.372
5	rs516675337	g.112389173C>G	CHADL	P: V445L	3/5	Gtg/Ctg	0.725
5	rs520244098	g.112389376T>C	CHADL	P: D377G	3/5	gAc/gGc	0.812
5	rs714095032	g.112389392G>C	CHADL	P: P372A	3/5	Ccc/Gcc	0.107
6	rs714949670	g.70205294A>G	KIT	P: D60G	2/21	gAt/gGt	0.217
6	rs109630427	g.70214244T>C	KIT	P: M258T	5/21	aTg/aCg	0.939
7	rs520476700	g.50739789G>T	STING1	P: L201I	6/8	Ctc/Atc	0.928
13	rs211162390	g.57531848A>G	GNAS	P: V144A	1/13	gTg/gCg	0.894

).

Last but not least, we assessed the nucleotide diversity and genome differentiation of genes with missense mutations, e.g., GNAS (Figure 8a) and MZB1 (Figure 8b).

4. Discussion

This study analyzed WGS data from 14 Ethiopian indigenous cattle populations living in and adapted to different agro-ecologies, climatic zones, and topographies in the country. Ethiopia, located in the Horn of Africa, was the entry point of several livestock species into the African continent, in particular, the humped indicine cattle B. t. indicus domesticated in the Indian subcontinent [4,5,8]. We aimed to study the genetic structure and selection signature of Ethiopian cattle populations compared with European taurine cattle, African zebu (e.g., Sudanese zebu), and African sanga cattle populations.

As observed in previous studies, the PCA (PC1) separated the zebu from taurine cattle [4,13]. Interestingly, Ethiopian cattle populations were closely grouped with Sudanese zebu, despite the huge differences in their habitats (e.g., drylands for Sudanese zebu versus dry, hot, and humid lowlands as well as cold and humid highlands for Ethiopian cattle). This result was supported by the population admixture at K = 4. Ethiopian cattle populations shared their major ancestry with other African zebu, followed by Asian zebu and African taurine cattle, while a limited proportion of European taurine ancestry was also observed. Sheko cattle comprised a relatively high African taurine ancestry, followed by Mursi cattle (Figure 3b), and both cattle populations have been reported to be relatively adapted to trypanosome infection [73,74].

The genomic diversity in Ethiopian cattle populations was estimated based on the genome-wide nucleotide diversity, ROH-based genomic inbreeding coefficient, heterozygosity, and population genetic differentiation. Considering the history of introduction and crossbreeding of indicine cattle into the African continent, Ethiopian cattle are expected to have very rich genetic diversity, as illustrated by their higher nucleotide diversity compared to European taurine, African taurine, and Asian zebu cattle (Figure 4 and Table 2). Likewise, previous reports have revealed similarly high nucleotide diversity in several Ethiopian cattle populations based on mitochondrial DNA [21,75], high-density SNPs [76], and WGS data [17,77]. Ethiopian cattle were genetically differentiated from African and European taurine cattle, Asian zebu, and admixed Ankole cattle (Table 3) as a result of geographic isolation, limited introgression from commercial breeds, and restricted movement of cattle within and across the country [13,78]. Low population genetic differentiation among Ethiopian cattle populations suggests that they shared a recent common history of ancestry or admixture. Nevertheless, they have likely evolved rapidly to adapt to contrasting environments such as the hot and dry lowlands and, more importantly, the Ethiopian highlands characterized by a cold and humid environment (Figure 1 and Table 1).

Identification of genomic regions that contain a series of homozygous loci inherited together from a common ancestor can better measure inbreeding level than genetic heterozygosity and pedigree information [35,79,80,81]. We estimated the genomic inbreeding based on ROHs ≥ 0.5 Mb (Figure 5 and Table 2), of which a few ROHs were observed in Ethiopian cattle compared to European and African taurine cattle (Table 2). The accumulation of short ROHs in Ethiopian cattle populations show low genomic inbreeding, while long ROHs in non-Ethiopian cattle populations reflect high genomic inbreeding in the current generation [80,82]. The low genomic inbreeding coefficients observed in Ethiopian cattle populations were also supported by their high heterozygosity, nucleotide diversity (Table 2), and effective population size (Figure 6b, Table S2).

Several candidate genes were commonly detected in Ethiopian cattle using up to four genome scan methods. Among them, THRA, MRD1, STING1, CHADL, CTNNA1, MZB1, NR1D1, and GNAS were identified by at least three of the four methods. Besides these genes, KIT, IL18, PRKCZ, PLA2G2A, and ITPR2 were functionally annotated in many GO clusters (p < 0.05). The DAVID bioinformatics analysis clustered their functions into skeletal system development, immune response, lipid storage, cellular response to misfolded protein, pancreatic secretion, renin secretion, and vascular smooth muscle contraction, which are important genetic controls related to environmental adaptation [72,83].

4.1. Selection Signatures for Immune Response

Candidate genes such as IL18, PRKCZ, STING1, and MZB1 are related to the immune response. MZB1 spans a 2.0 kb region on BTA7. Its function is related to the endoplasmic reticulum calcium ion homeostasis, cellular signal transduction, B cell activation, and T cell proliferation [84]. STING1 is located in the same genomic window as MZB1 based on the Hp, F_ST, and XP-CLR analyses. Its function is related to the innate immune signaling response to protect the host against viral and bacterial infections [85,86]. In the STING1 gene region, we detected rs520476700, a missense SNP at 7:g.50739789G>T with a frequency of 0.93 (Table 6). The biological process of STING1 is related to the immune system development in response to UV light by activating the production of DNA-mediated type I interferon [87], which plays a significant role in cattle adaptation to environmental stresses.

4.2. Selection Signatures for Growth and Production

Other candidate genes (e.g., GNAS, CHADL, and THRA) are related to growth and production traits. They are involved in skeletal system development, carcass characteristics, milk production, and reproduction. The GNAS and ENSBTAG00000017475-associated complex are involved in regulating metabolism responses [88], components of complex cattle performance traits, reproduction, and health [89]. They regulate insulin secretion in pancreatic β-cells relevant to milk yield, milk protein and fat contents, growth, and carcass quality in cattle [89,90,91]. Missense SNPs were found in the first exon of GNAS (rs211162390, 13:57531848A>G) and the third exon of CHADL (rs520244098, 5:g.112389376 T>C) with frequencies in Ethiopian cattle of 0.89 and 0.81, respectively (Table 6). In humans, mutations in GNAS were associated with reduced growth by disrupting the growth hormone-releasing hormone receptor signaling pathways [92]. CHADL is annotated in the GO clusters associated with skeletal system development (p < 1.2 × 10⁻⁷) and immune response (p < 6.8 × 10⁻³) (Table 5). Previous reports showed its involvement in cartilage development and cell modification during skeletal system development [51], which are significantly related to growth and carcass production.

4.3. Selection Signatures for Reproduction

The candidate gene KIT (BTA6: 70.1–70.2 Mb) was found to be associated with skeletal development and the positive regulation of immune responses (Table 5 and Table S7), as well as male reproduction functions in spermatogonia cell differentiation [93] and germ cell maturation [94]. Several reports also indicated that KIT mutations were responsible for the white spotting coat color pigmentation in cattle [95,96] and horses [97] and for coat color in pigs [72]. All these functions are of significance to the adaptation to hot tropical environments [98]. The candidate gene RARA (BTA19: 50.55–50.59 Mb) is involved in the biological processes of skeletal development, immune response, regulation of hemopoiesis, and estrogen signaling pathway (Table 5 and Table S7). The RARA gene is also involved in germ cell organization, meiotic integrity, and spermatogonia cell proliferation, functions critical to spermatogenesis and fertility [99,100].

Another important candidate gene is NR1D1 (BTA19: 40.38–40.39 Mb). GO analysis clustered it into the molecular function of sequence-specific double-stranded DNA binding and the regulation of immune response (Table S7). NR1D1 regulates animal reproduction by controlling circadian clock gene expression [101,102,103] and inducing testosterone production in Leydig cells [104]. Variants identified in its third exon regulate genes involved in lipid metabolism, adipogenesis, gluconeogenesis [105], and macrophage inflammatory response [106]. It is also expressed in response to environmental stress, such as hypoxia, which may lead to pulmonary inflammation [107] and brisket disease in cattle [108].

4.4. Selection Signatures for Environmental Adaptation

GO clusters, such as cellular response to misfolded protein, the stress response to copper ion, and vascular smooth muscle contraction, regulate cell homeostasis in response to environmental stressors such as heat, chemicals, and altitude. ITPR2 is involved in pancreatic secretion, renin secretion, and vascular smooth muscle contraction pathways important to cardiovascular function in response to circulatory homeostasis, calcium signaling pathways, cellular water imbalance, and adaptation to hypoxia [59,109,110]. In addition, DNAJC18 (BTA7: 50.6–50.7 Mb) and BAG6 (BTA23: 27.64–27.65 Mb) genes are found to be responsible for misfolding proteins that play a significant role in heat stress response [111,112].

5. Conclusions

The intensive analysis of whole genome sequences presented here provides new genomic insights into the genomic diversity and selection sweeps of Ethiopian indigenous cattle populations. The population genetic structure and admixture pattern indicate that Ethiopian cattle have a higher zebu ancestry proportion in comparison to other African indicine cattle. Ethiopian cattle are also characterized by their high nucleotide diversity, heterozygosity, and effective population size but limited genetic differentiation and genomic inbreeding. The genome-wide scan for signatures of selection has identified several candidate genes associated with biological functions such as fat metabolism, immune response, water balance, coat color pigmentation, and circulatory homeostasis, which may have facilitated the rapid adaptation of Ethiopian indigenous cattle populations to the diverse and challenging environments of the country.