Details of the cattle populations studied are presented at Table 1. They include 92 non-European introgressed small EASZ from the Western and Nyanza provinces of Kenya (KEASZ) as well as cattle from Uganda, Nigeria, Guinea, Europe, and India. They were genotyped for 777,962 SNPs mapped to the UMD3.1 bovine reference genome using the Illumina BovineHD Genotyping BeadChip.

Quality control (QC) analyses for 735,297 autosomal SNPs were conducted through the check. marker function implemented in the GenABEL package (Aulchenko et al., 2007) for R software version 2.15.1 (R development Core Team, 2012). SNPs with a minor allele frequency (MAF) <0.05 (n = 68,731) or call rate <95% (n = 18,667) were filtered out from the entire dataset. These include 1,712 SNPs that failed both criteria. In total, 649,611 SNPs were therefore retained. The ancestral allelic state for 373,005 SNPs [mean gap size = 6.7 kb and standard deviation (SD) = 12.1 kb] has been reported previously by Utsunomiya et al. (2013) following genotyping of three non-cattle Bovinae species: two B. gaurus (gaur), six Bubalus bubalis (water buffalo), and two B. grunniens (yak), with the fixed allele in the three species considered as ancestral. Only these SNPs were included in the downstream analyses.

Additional QC criteria included a minimum sample call rate of 95% and a maximum pairwise identity-by-state (IBS) of 95%, with the lower call rate animal eliminated from the high IBS pair. One ZS sample was excluded for having a low call rate, whilst 15 samples (two GIR, one NEL, four HOL, four JER, two AG, one AZ, and one WD) were excluded following the IBS criterion.

PCA were conducted using the prcomp function implemented in GenABEL package for R software. These analyses were carried out in three levels; (i) all cattle populations, (ii) all cattle populations excluding the European taurine cattle, and (iii) only the KEASZ, Uganda (UGN) and Nigeria (NGR) cattle.

Admixture analysis using ADMIXTURE 1.23 software (Alexander et al., 2009) with cross-validation and 200 bootstraps for K = 3 was conducted on the whole dataset to determine the European taurine, Asian zebu, and African taurine ancestries at genome-wide level and for each autosome separately. The output files were graphically displayed by the ggplot2 package (Wickham, 2009) for R software.

Rsb analyses (Tang et al., 2007) were conducted between each of the KEASZ, combined UGN cattle populations (AO, KR, NG, ZS) and combined NGR cattle populations (AG, AZ, BJ, OR, SO, WD, YK; Tijjani, 2013; Mbole-Kariuki et al., 2014) with the combined reference cattle populations (NEL, GIR, NDM, MT, HOL, JER) using the rehh package (Gautier and Vitalis, 2012) for R software. The standardized Rsb-values were normally distributed (Supplementary Figure 1), so a Z-test was applied to identify statistically significant SNPs under selection on KEASZ, UGN, and NGR cattle populations. One-sided upper-tail P-values were derived as 1-Φ(Rsb) from the Gaussian cumulative density function Φ. Candidate regions were defined as having five adjacent SNPs, not separated by more than 500 kb, passing the threshold of −log₁₀ P = 4.

iHS analyses (Voight et al., 2006) were conducted on KEASZ, combined UGN cattle and combined NGR cattle populations using the rehh package for the R software. This statistic was calculated for SNPs that passed the QC criteria and exhibited a within-population MAF of at least 0.05, since the algorithm of iHS has a limited power to calculate the statistic for fixed SNPs. As with Rsb, the standardized iHS-values followed a normal distribution (Supplementary Figure 1), so a two-tailed Z-test was applied to identify statistically significant SNPs under selection with either an unusual extended haplotype of ancestral or derived alleles relative to the genome. Two-sided P-values were derived as 1-2|Φ(iHS)-0.5| from the Gaussian cumulative density function Φ. Candidate regions were defined as in Rsb.

As a prerequisite to the Rsb and iHS analyses, fastPHASE 1.4 (Scheet and Stephens, 2006) was used to phase the genotyped SNPs into the corresponding haplotypes using K10 and T10 criteria. Population label information was used to estimate the phased haplotype background.

ΔAF analysis investigates absolute allele frequency difference between two populations (Carneiro et al., 2014). In this analysis, the mean frequency of the first allele was estimated for KEASZ, combined UGN cattle and combined NGR cattle populations, separately (AF_pop1). Likewise, the mean frequency of the first allele for each SNP was calculated for the combined reference cattle populations (AF_pop2). The standardized values of the ΔAF (AF_pop1–AF_pop2) were normally distributed (Supplementary Figure 1), therefore a Z-test was applied to identify statistically significant SNPs showing higher allele frequency in the first population. Two-sided P-values were derived as 1-2|Φ (standardized ΔAF)-0.5| from the Gaussian cumulative density function Φ. Candidate regions were defined as in iHS and Rsb.

The Stouffer method was used to combine the P-values obtained from Rsb, iHS, and ΔAF, for each analyzed set of populations, i.e., KEASZ, UGN, and NGR, in meta-SS analyses (Whitlock, 2005; Utsunomiya et al., 2013). Each value, for every SNP in each test, was transformed to a Z-score, Φ⁻¹ (1-P-value), where Φ⁻¹ is the inverse cumulative distribution function for a standard normal distribution. Then, the SNP-specific Z-scores were combined together according to the following equation: Z_i = (Z_Rsb + Z_iHS + Z_ΔAF)/k, where i and k are numbers of SNPs and tests, respectively. The resulting Z-scores were referred back to the standard normal distribution to obtain combined P-values [P-value = 1- Φ(Z_i)], where is the cumulative distribution function for a standard normal distribution. Candidate regions were defined as having five adjacent SNPs not separated by more than 500 kb passing the threshold of −log₁₀ P = 4.

A single pool of 10 unrelated KEASZ DNA samples were sequenced using an ABI SOLiD 4 genetic analyser (Supplementary Table 1). SOLiD 2 × 50 bp mate-paired libraries were constructed and sequenced in the Deep Seq facility at the University of Nottingham according to the manufacturer's instructions. Reads were mapped to the UMD3.1 bovine reference genome assembly (Elsik et al., 2009) using the LifeScope Genomic Analysis software 2.5.1 re-sequencing mapping pipeline (http://www.lifetechnologies.com/lifescope). SNPs and indels (insertions/deletions) were called using the diBayes package implemented in LifeScope. A minimum coverage of two uniquely mapped reads and two non-reference allele counts were required to call a variant. Additionally, a minimum read mapping quality of 20 (MAPQ20) and base quality of 20 were also implemented.

Pooled heterozygosity Hp of SNPs detected in the pooled 10 KEASZ full genome SOLiD sequences were calculated on 100 kb sliding windows with 10 kb incremental steps. Window sizes were extended by the number of uncovered bases to improve the accuracy of the calculation and consistency across windows. For each SNP in the window, the number of reads for the most and least frequent allele was counted (n_MAJ and n_MIN, respectively). Hp-values were calculated using the following formula: Hp = 2 ∑ n_MAJ ∑ n_MIN/(∑ n_MAJ + ∑ n_MIN)² (Rubin et al., 2010). The autosomal Hp-values were Z-transformed (ZHp = Hp–mean Hp/SD Hp) and a ZHp ≤ −4 was applied as a threshold to specify windows carrying selective sweep as in Liao et al. (2013). Overlapping candidate windows were merged into a single region.

The Agilent SureSelect^XT target enrichment kit (cat no. G7530-90004) was used for bovine exome sequence enrichment. It covers a total of ~45 Mb of the bovine sequence that composed of the coding regions in the UMD3.1 reference genome (coding regions from Refseq and Ensembl, no UTR “untranslated regions”) and microRNA. The exomes of a further 10 KEASZ samples (Supplementary Table 1) were sequenced in the Deep Seq facility at the University of Nottingham using an ABI SOLiD 5500 genetic analyser. Using the LifeScope Enrichment Sequencing Pipeline, the generated 75 bp reads were mapped to the UMD3.1 bovine reference genome. Supplementary Table 2 summarizes the number of aligned reads, average depth of coverage, and percentage of reference exome for each sample. The same variants calling criteria used in the KEASZ whole genome sequence were implemented for the exome analysis except that reads with MAPQ ≥ 30 and base quality of 30 were used.

To infer about the origin of the variants identified in the KEASZ selected regions, we analyzed 10 N'Dama cattle full genome sequences from Guinea available at the GenBank with the Bioproject accession number PRJNA312138 (Kim et al., 2017), and 10 Muturu cattle genome sequences from South East Nigeria. For the later, genomic DNA was extracted using Macherey-Nagel NucleoSpin® Tissue DNA extraction kit, according to manufacturer's protocol, from 10 ear tissue samples. One hundred and fifty bp paired-end libraries were constructed and sequenced using Illumina Hiseq2500 platform (Illumina, USA). Sequence reads were mapped to the UMD3.1 bovine reference genome assembly using the bwa-mem option of Burrows-Wheeler Alignment tool (BWA) version 0.7.5a (Li and Durbin, 2010; Supplementary Table 3). Genome Analysis Toolkits (GATK) version 3.4.0 (McKenna et al., 2010) was used for base quality recalibration, indel realignment, PCR duplicate removal and variants (SNPs and indels) calling, using Haplotypecaller function, according to GATK best practices recommendations (DePristo et al., 2011; Van der Auwera et al., 2013). A minimum read mapping quality of 50 (MAPQ50) and base quality of 30 were set as two criteria in variant calling.

Protein-coding and RNA genes mapped within the candidate regions, based on the UMD3.1 bovine reference genome annotation, were processed using the functional annotation tool implemented in DAVID Bioinformatics resources 6.7 to determine the over-represented (enriched) functional terms (Huang da et al., 2009a,b). An enrichment score of 1.3, which is equivalent to the Fisher exact test P = 0.05, was used as a threshold to define the significantly enriched functional terms in comparison to the whole bovine reference genome background. The list of genes mapped on the UMD3.1 reference bovine genome was obtained from the Ensembl Genes 86 database (Flicek et al., 2013) using the BioMart tool (Kinsella et al., 2011). The bovine Quantitative Trait Loci (QTL) and their UMD3.1 genome coordinates were downloaded from the cattle QTL database (http://www.animalgenome.org/cgi-bin/QTLdb/BT/index). The intersectBed function from the BedTools software was used to overlap these QTL with the identified KEASZ candidate regions (Quinlan and Hall, 2010). Candidate regions lacking annotated genes (protein-coding, RNA, and pseudogenes) based on UMD3.1 bovine reference genome were also identified. The intersectBed function from the BedTools software was also used to overlap these candidate regions with transcription factors binding sites identified previously on the bovine reference genome by Bickhart and Liu (2013).

Variants (SNPs and indels) in the genes within the genome-wide SNP and Hp overlapping candidate regions were annotated using the variant effect predictor tool (McLaren et al., 2016) based on Ensembl variants 86 databases. Comparisons with the previously discovered bovine variants listed in the dbSNP database (Sherry et al., 2001) classified these variants into KEASZ-specific (novel) and general bovine variants. The identified missense, splice region SNPs, frameshift and splice region indels were cross-checked with variants identified in the KEASZ exome data and with N'Dama and Muturu full genome sequences.

Candidate genes within the overlapping genome-wide SNP and Hp candidate regions were selected based on their biological function, e.g., immunity, reproduction, fertility, heat tolerance and anatomical development. Variants (SNP and indels) within these genes were defined and the biological effect of the missense SNPs were predicted by the online tool PolyPhen-2 (Adzhubei et al., 2010).

LAMP software version 2.4 (Sankararaman et al., 2008) was used to estimate the Asian zebu and African taurine ancestry proportions of the high-density genotyped SNPs in the KEASZ samples. The genome-wide autosomal zebu ancestry proportion of 70% was obtained from the admixture proportions α of the ADMIXTURE analyses. An estimated number of 500 generations was set for the beginning of the zebu-taurine admixture in light of our current knowledge of zebu arrival on the continent, assuming a generation time of 6 years (Keightley and Eyre-Walker, 2000). A uniform recombination rate of 1 cM = 1 Mb was assumed. The average excess/deficiency in Asian zebu ancestry at each SNP (ΔAZ) was calculated by subtracting the average estimated Asian zebu ancestry of the SNP from the average estimated Asian zebu ancestry of all SNPs. The median ΔAZ for the significant SNPs of KEASZ meta-SS analysis within candidate regions was considered.

The PCA plot on all cattle populations using the filtered autosomal SNPs dataset (Figure 1A) reveals the already described triangle-like 2-dimentional global organization of cattle genetic diversity (Gautier et al., 2010). The first component, which accounts for 20.6% of the total variation, separates the European (HOL and JER) and African (NDM and MT) taurine populations from the Asian zebu cattle (NEL and GIR). The second component, which accounts for 4.9% of the total variation, separates the African taurine cattle from the other cattle populations. All the East (KEASZ and UGN) and West (NGR) African cattle populations analyzed are positioned between the African taurine and Asian zebu cattle supporting zebu × taurine genome admixture.

The PCA conducted on Asian zebu (NEL and GIR), African taurine (NDM and MT), East African (KEASZ and UGN) and West African (NGR) cattle populations (Figure 1B) shows genetic differentiation within the East African cattle (KEASZ and UGN) cattle and between the East and West African (NGR) cattle. The first component, which explains 13.9% of the total variation, separates Asian zebu from African taurine. Whilst, the second component, which explains 2.6% of the total variation, separate West African zebu cattle (NGR) and the East African (KEASZ and UGN) cattle. A PCA on KEASZ, UGN, and NGR cattle only, separates first the East (KEASZ and UGN) and West (NGR) cattle populations (PC1 = 2.99%), while PC2 (2%) reveal a much higher genetic heterogeneity within and across East African cattle populations compared to the West African ones (Supplementary Figure 2).

The admixture analysis indicates a relatively even admixed genome of Asian zebu and African taurine ancestries across KEASZ animals with an estimated genetic proportion of 0.7 ± 0.01 SD and 0.30 ± 0.01 SD, respectively. The same even admixed proportion is present across the Karamajong zebu cattle from Uganda, while the Serere zebu has an estimated 0.74 ± 0.04 SD and 0.24 ± 0.04 SD zebu × taurine genetic proportions, with minor European taurine ancestry (0.02 ± 0.01 SD). The genome of the Ankole cattle has a larger African taurine genetic proportion (0.48 ± 0.01 SD). Nganda cattle has an admixed genome of 0.57 ± 0.04 SD Asian zebu and 0.36 ± 0.03 SD African taurine ancestries, it also carries low level of European taurine ancestry (0.07 ± 0.05 SD) (Figure 2). This pattern of admixture is also observed in the cattle of Nigeria but with higher average African taurine ancestry (0.35 ± 0.01 SD) and lower zebu ancestry (0.65 ± 0.01 SD) than the East African cattle, except Ankole. No European introgression was detected in the West African zebu cattle analyzed.

Variation in Asian zebu ancestry is observed among the autosomes of the East (KEASZ and UGN) and West African (NGR) cattle populations (Supplementary Figure 3, Supplementary Table 4) with increase or decrease (mean ± one SD) in Asian zebu ancestry proportions (Table 2). A strong correlation in the zebu ancestry autosomal variation is observed when comparing KEASZ with East (UGN) and West (NGR) cattle populations (P < 0.00001; Supplementary Figure 4).

LAMP software 2.4 (Sankararaman et al., 2008) estimated the mean Asian zebu ancestry proportion for all autosomal SNPs in KEASZ to be 0.76 (SD = 0.14). Based on this estimation, the mean ΔAZ for autosomal SNPs is 0 (SD = 0.14). The majority of the candidate sweep regions show high zebu ancestry proportion, but similar to the mean Asian zebu ancestry proportion (within one SD from the mean). For the 32 East African candidate regions (Table 3), eight regions reveal substantial ΔAZ, a single region shows deficiency and seven show excesses (more than ± one SD from the mean ΔAZ). Moreover, six East and West African candidate regions (Table 3) demonstrate substantial ΔAZ. Two of these regions show deficiencies and four show excesses of Asian zebu ancestry. When the genome-wide SNP and Hp analyses overlapping KEASZ candidate regions were considered, 13 regions demonstrate substantial ΔAZ. One region, which is shared between East and West African cattle, show deficiency in Asian zebu ancestry. Whilst 12 regions show excess in Asian zebu ancestry, in which six are specific to East African cattle and one is shared between East and West African cattle (Table 7).

Archaeological and genetic evidences so far indicate that the history of African zebu cattle started with the introgression of Asian zebu to the native African taurine populations ~4,000 and 1,300 years ago (Epstein, 1971; Hanotte et al., 2002). The coordinates of the African cattle samples (KEASZ, UGN, and NGR) in the PCA plots (Figures 1A,B) indicate zebu × taurine admixture level in their genome. This zebu introgression appears even across animals within population on the autosomes of KEASZ and also in West African cattle from Nigeria as well as in Ankole and Karamojong zebu cattle from Uganda (Figure 2). At the contrary, uneven European taurine introgression, likely of recent origin following ongoing crossbreeding of indigenous African cattle populations with exotic cattle breeds to improve their productivity (Mwai et al., 2015), is observed in the genome of Nganda and Serere zebu cattle from Uganda (Figure 2).

Interestingly, we also observe positive significant correlation in Asian zebu ancestry autosomal proportion between the KEASZ and the other East and West African cattle populations examined. This is in agreement with the known history of zebu cattle on the continent (Hanotte et al., 2002) and it supports a common ancestry for the African cattle examined here. Our results also support an East African origin for the West African zebu cattle as previously showed in Hanotte et al. (2002) with migration of admixed zebu × taurine cattle populations.

This study is the first to our knowledge that exhaustively investigated the genome of an indigenous African cattle population for signatures of positive selection using both high density genome-wide SNP and full genome sequence information. The outputs of these analyses represent a follow up of our previous identifications of signatures of selection on KEASZ genome using the lower density Illumina BovineSNP50 BeadChip v.1 (Bahbahani et al., 2015). Here, we have validated 14 candidate regions out of 24 regions previously identified to be targeted by positive selection in KEASZ (Supplementary Table 17; Bahbahani et al., 2015). The remaining 10 regions might have been false positives. Indeed, in our previous study, we used a genomic tool characterized by high European taurine ascertainment bias and low genome coverage. This issue has been addressed in this study by using high-density HD SNP array and full genome sequence data. Moreover, we previously selected candidate regions based on only two SNPs above the significant threshold. Here we have been using much more stringent criteria with a minimum of five SNPs above the threshold to define a candidate region under positive selection (see Materials and Methods Section).

We also used a different strategy to increase the power of detecting genomic signatures of positive selection. First, instead of analyzing regions defined by each genome-wide SNP analyses (Rsb, iHS, and ΔAF) separately, a composite statistical approach “meta-SS analysis” was conducted at the autosomal level as in Utsunomiya et al. (2013) to combine the SNP-specific P-values for each test into a single index.

Simulation data have shown that combining the signals from different tests into a single statistic increase the power of defining genomic regions under selection (Grossman et al., 2010). Most importantly, coalescent simulations and accurate calibrated demographic models, which are lacking in African cattle, are not required by meta-SS in comparison to the original Composite of Multiple Signals (CMS) method proposed by Grossman et al. (2010). Secondly, the reference cattle populations were pooled into a single population as in Bahbahani et al. (2015). As suggested by Gautier and Naves (2011), this pooling approach increases the haplotype diversity in the reference populations and it breakdown population-specific linkage disequilibrium (LD) which may result from genetic drift. Thirdly, information from both genome-wide SNP genotypes and full genome sequence were used to adequately cover the KEASZ genome, as well as, to address any breed ascertainment bias associated with the commercially available SNP arrays (Matukumalli et al., 2009). Overlapping candidate regions from these two approaches further support the identification of the candidate regions.

Both selection and genetic drift may have shaped the genome of the KEASZ and the other populations examined here. Distinction between the two is difficult, and fixation or near fixation of allele through genetic drift may lead to false positive candidate regions for signatures of positive selection (Qanbari and Simianer, 2014). Comparison of the results between different cattle populations, the unraveling of the zebu or taurine origin of the selected regions and the function of genes within selected regions may here provide further information.

In this study, we observe a small number of candidate regions for positive selection shared between KEASZ HD SNP and/or genome sequencing information and the East and West African cattle population examined here (Tables 3, 7). While it may be argued that these overlapping candidate regions for positive selection are a consequence of common genetic backgrounds, the low number of shared candidate regions as well as the identification of the same regions in tropical-adapted cattle populations with different population histories, e.g., Creole cattle from Guadeloupe (Gautier and Naves, 2011), zebu × taurine admixed cattle from West Africa (Gautier et al., 2009; Flori et al., 2014; Xu et al., 2015), Brahman (Ramey et al., 2013; Xu et al., 2015) and Gir (Liao et al., 2013; Perez O'Brien et al., 2014) cattle (Tables 3, 7, and Supplementary Table 11) strongly supports the role of selection rather genetic drift for these shared regions. It underlines the importance that these genome regions may play a role in the adaptive traits of African tropical-adapted admixed cattle.

Examining the probable ancestral origin of the candidate regions reveals a subset of candidate regions with substantial ΔAZ (Tables 3, 7). Most of these regions show excesses of Asian zebu ancestry, e.g., BTA 7: 51.4–53.4 Mb and BTA 13: 47.5–48.1 Mb, indicating that the indicine haplotypes are more likely to be under selection in the African admixed cattle populations than the taurine. This is perhaps not surprising considering the predominant zebu genomic background in EASZ (Mbole-Kariuki et al., 2014). However, we also observed candidate regions showing substantial excess of African taurine ancestry, e.g., BTA 3: 76.1–76.4 Mb and BTA 19: 3.3–3.8 Mb (Table 3). These are present in chromosomes with overall low level of African taurine ancestry adding further support for selective pressure rather than genetic drift for their presence (Table 2).

Examining the potential biological functions under positive selection reveal several different significantly enriched biological pathways likely under selection in KEASZ given their importance for a cattle population living in a challenging tropical environment (e.g., bovine adaptive and innate immunity, response to hormone stimuli, intermediate filaments, and keratins pathways). Additionally, genes and QTL related to regulation of bovine immunity, fertility and reproduction, anatomical development, and heat stress have also been found within the identified KEASZ candidate regions (Supplementary Tables 9, 12, 18).

Indeed, innate and adaptive immune genes may be expected to be primary targets of selection in African cattle that are exposed to a diversity of pathogens and associated physiological stresses in their surrounding environment, e.g., endoparasites, haemoparasites, and bacteria (de Clare Bronsvoort et al., 2013; Murray et al., 2013; Thumbi et al., 2014). Examples of candidate genes related to this category are: C-C chemokine receptor type 7 precursor (CCR7) and granulocyte macrophage-colony stimulating factor (CSF3). Upon binding to two chemoattractants: CCL19 and CCL21, CCR7 is involved in maturating dendritic cells and hence activate T lymphocytes (Marsland et al., 2005; Forster et al., 2008). This receptor has also demonstrated a role in regulating innate immunity by attracting macrophages to sites of infection (van Zwam et al., 2010). The multifunctional cytokine (CSF3) acts as a positive regulator for macrophages to induce their antimicrobial effects (Grabstein et al., 1986; Tarr, 1996). Interestingly, several trypanotolerance QTL, identified by Hanotte et al. (2003), overlap with the identified KEASZ candidate regions (Supplementary Table 19). These QTL might indicate an undocumented level of trypanotolerance in KEASZ, as it has already been shown in other East African cattle populations (e.g., Orma Boran, Sheko and Mursi cattle; Dolan, 1987; Mwangi et al., 1993; Bahbahani and Hanotte, 2015).

Our results also show that besides tolerance to disease challenges, fertility and reproduction traits have been also selected in the African zebu × taurine admixed cattle. Examples of these candidate genes are: the retinoic acid receptor α subunit (RARA), relaxin/insulin-like family peptide receptor 2 (RXFP2), spermatogenesis associated 24 (SPATA24), and sperm associated antigen 7 (SPAG7). The retinoic acid receptor, which is expressed in sertoli cells in the seminiferous tubules, plays a role in maintaining retinoic acid signal to induce spermatogonia differentiation (Wolgemuth and Chung, 2007). RXFP2, which has been found to be under selection in admixed Creole cattle (Gautier and Naves, 2011) and Gir cattle (Liao et al., 2013), plays a role in testicular descent development (Gorlov et al., 2002; Agoulnik, 2007; Park et al., 2008; Feng et al., 2009). Candidate genes such as spermatogenesis-associated 24 (SPATA24) and sperm-associated antigen 7 (SPAG7) can also be classified into the cattle fertility and sexual reproduction category due to their role in spermatogenesis. Such signals may be the legacy of selection for fertility in hybrid populations between two cattle lineages (zebu and taurine), with a common ancestry perhaps as old as half a million years ago (MacHugh et al., 1997). This requires further investigation. QTL related to bovine reproduction and fertility, e.g., sperm motility and calving ease, are also target of positive selection in KEASZ.

An important gene category identified within the candidate regions is the heat stress category. Genes and QTL within this category might be targeted by natural selection to adapt to the tropical environmental condition (Hansen, 2004). Two heat shock protein genes (HSPA9 and HSPB9) and two members of the DnaJ family (DNAJC8 and DNAJC18) are mapped within three of the identified candidate regions (Table 8). Heat shock proteins have critical roles in maintaining protein folding under heat stress (Parsell and Lindquist, 1994; Coleman et al., 1995). The members of DnaJ family act as cofactors for other heat shock proteins (Hsp70) to maintain protein folding (Kampinga and Craig, 2010).

Genes and QTL related to anatomical development were also identified within the KEASZ candidate regions (e.g., LEMD3, LOX, and RXFP2). These genes are important to maintain optimum growth and development. LEMD3 and LOX are associated with the development of different organs, such as heart (LEMD3), lung and blood vessels (LOX) (Maki et al., 2005; Ishimura et al., 2008). The candidate region harboring LEMD3 has also been found to be under selection in Brahman cattle (Ramey et al., 2013). The role of RXFP2 in testicular descent development (Gorlov et al., 2002; Agoulnik, 2007; Park et al., 2008; Feng et al., 2009) can classify this gene into the anatomical development, fertility and reproduction as well as the heat tolerance categories. This gene has also been associated with the horn phenotype in sheep (Johnston et al., 2011; Kijas et al., 2012) and a study by Johnston et al. (2013) has demonstrated an association between variants of this gene and reproductive success and survival rate in Soay sheep from St. Kilda.

In addition, various production traits QTL, e.g., marbling score, milk fat percentage, and milk fat yield have also been found within the KEASZ candidate regions. Given that several of the identified KEASZ candidate regions overlap with regions under positive selection in commercial dairy and beef cattle breeds (Tables 3, 7, Supplementary Table S11; Larkin et al., 2012; Kemper et al., 2014; Qanbari et al., 2014), this may illustrate possible human selection for some production traits, e.g., milk yield, has taken place in KEASZ at least in the past.

In parallel to candidate regions harboring genes, the 39 gene desert candidate regions on KEASZ are targets of further research to define their biological roles. Although these regions do not contain any of the bovine transcription factors binding sites identified by Bickhart and Liu (2013), they may still harbor unannotated regulatory elements and/or genes targeted by positive selection. These regions may also be transcribed to generate long non-coding RNA transcripts (≥200 nucleotides), which could be further validated by RNA sequencing. This type of RNA molecules has critical roles in regulating the expression of neighboring genes at transcriptional and post-transcriptional levels (Mercer et al., 2009; Wang and Chang, 2011).

Several genomic variants, in which some are unique so far to KEASZ, have been identified within the KEASZ overlapping genome-wide SNP and Hp candidate regions. A subset of these variants, such as missense SNPs, framshift indels, and splice region SNPs/indels, can be considered as primary target of positive selection due to their functional roles on the corresponding genes. Interestingly, within the selected candidate genes four missense SNPs, a splice region SNP and a frameshift indel have been identified. These variants need to be confirmed in larger sample size of tropical-adapted (e.g., EASZ) and non-tropical-adapted (e.g., European taurine) cattle to further support their role as causative variants under selection. Two of the missense variants are predicted to show non-benign effects on their genes calling for their effect to be validated, e.g., through gene editing approaches (Carlson et al., 2012; Wang et al., 2013). Moreover, variants on the other genes within the candidate regions or in non-coding regions following their possible roles in regulating gene expression may also be targets of selection.

Based on available sequences of KEASZ and African taurine cattle (N'Dama and Muturu), we were able to infer the possible ancestral origin of some of these functional variants (Asian zebu or African taurine origin) indicating the role of the zebu × taurine admixture as a selective force shaping the genome of KEASZ. However, this origin assignment was only partially successful due to the low number of sequenced samples and the unavailability of a full indicine de novo reference genome.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer MB and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.