Study area

Blood samples were collected from goats in three counties in Kenya: Nyeri (Mukurweini Sub-County), Meru (Central Imenti Sub-County) and Homa Bay (Homa Bay Town Sub-County) located in the Central, Eastern and Western regions of Kenya, respectively (Figure 1). These areas were selected because they represent the entry points of different imported dairy goat breeds in Kenya. Mukurweini Sub-County lies in the Upper midlands, also known as the main coffee zone, at an altitude of 1460-1710 metres above sea level (masl) and receives 950-1500 mm of mean annual rainfall. Central Imenti is in the upper highlands, at an altitude ranging between 1830-2210 masl and has an average annual precipitation of 800-2600 mm. The Homa Bay Town Sub-County lies in the lower midlands, at 1166 masl and receives an annual rainfall of 1226 mm.

Animal resources and sampling

Goat keeping households were purposively selected based on the following criteria: 1) being a member of dairy goat farmer group and 2) having more than two mature does which matched the breed kept by the farmer group in the said county of study. The herd structures between the breeds and within the county/breed varied among the selected households. Therefore, when a farm had only two mature does, only one doe was sampled. Where more than two does were available, the relationship of the does was confirmed by the farmer to avoid sampling closely related does. Sampling of full and half siblings was avoided. To ensure the representativeness of sampling for each breed, unrelated animals were selected from various farms across the designated counties. The Galla goat breed, however, did not follow the criteria because they were from the government breeding station where breeding records were used to identify the animal to be sampled. Therefore, to minimize sampling from closely related animals within the Galla population, pedigree data were used to select against full and half sibling animals.

A total of 96 animals including three imported breeds (31 Toggenburg, 29 Alpine, 24 Saanen) and one indigenous breed (12 Galla) were incorporated in this study. The Toggenburg and Alpine were found in Eastern and Central Kenya under Meru Goat Breeders Association (MGBA) and Dairy Goat Association of Kenya (DGAK), respectively. Saanen goats were found in Homa Bay under Nyanza Goat Breeders Association (NGBA). All samples were collected from a total of 53 farms in the three counties, Nyeri (18), Meru (19) and Homa Bay (16). Galla goats were sampled from the sheep and goats government station in Naivasha.

Whole blood (10 ml) was collected from the jugular vein into Vacutainer tubes with Ethylenediaminetetraacetic acid (EDTA) as an anticoagulant. The blood was stored at -20°C for two months before genomic DNA extraction. The procedure of blood sampling followed FAO guidelines (FAO, 2012). From each animal, duplicate samples were collected and kept separately during transportation and storage. For each sample, the following information was collected: sex of the animal, basic pedigree information, size of the flock, breed, any relevant phenotypic feature, and a photograph of the goat.

The study was conducted in strict accordance with the recommendations of the Institute of Primate Research (IPR) Ethical Guidelines on Animal Care and Use of Laboratory Animals (https://grants.nih.gov/grants/olaw/guide-for-the-care-and-use-of-laboratory-animals.pdf).

The protocol was approved by the committee on the ethics of Animal Experiment of Egerton University of Egerton in Kenya (ISERC/03/2020). A qualified veterinary officer collected the whole blood following FAO guidelines (FAO, 2012) to reduce pain and discomfort to a minimum.

DNA extraction and genotyping

Genomic DNA was extracted from the whole blood using the DNeasy Blood and Tissue Kit (Qiagen®, Hilden, Germany) according to the manufacturer’s instructions. Ten µl of each ten randomly selected samples were subjected to a preliminary estimate of the DNA quality and quantity on a 1 % agarose gel electrophoresis. Secondary quantification and purity analysis of the DNA were confirmed using one µl for each sample on both Nanodrop Spectrophotometer (Nanodrop® ND-2000) and Qubit ® dsDNA BR (Broad-Range) Assay Kit on the Qubit 3.0 fluorimeter (Invitrogen). The extraction and quality control check of genomic DNA was done at Kenya Agricultural and Livestock Research Organization biotechnology laboratory in Kabete, Kenya.

The DNA samples were genotyped using the GoatSNP50 Bead Chip (Illumina, Inc. San Diego, CA 92122 USA), developed by the International Goat Genome Consortium (IGGC), which features 53347 SNPs across the whole genome with inter-SNP spacing of approximately 40 kb (Tosser-Klopp et al., 2014). The genotyping was outsourced to Neogen Europe Limited in Scotland (https://genomics.neogen.com/en/).

SNP quality control and data analysis

The SNP genotype quality control process was applied to raw reads for both merged (all breeds) and then per breed using PLINK v1.07 (Purcell et al., 2007). First, individuals with a missing genotype call rate of more than or equal to 10 % were removed from further analysis using the mind function in PLINK with default settings. The remaining individuals were then exposed to further filtering. SNPs with less than 95 % call rate, Minor Allele Frequency (MAF) of less than 0.05 and P<0.001 Hardy Weinberg Equilibrium (HWE) were excluded from downstream analyses. The SNP data set used for downstream analysis is accessible from the Mendeley Digital Repository (https://doi.org/10.17632/hhb9rhdzzt.1).

Basic genetic diversity indices, which include the proportion of polymorphic markers, inbreeding coefficient, observed (H_o) and expected (H_e) heterozygosity were calculated within breeds using PLINK (Purcell et al., 2007). The proportion of polymorphic SNPs (P_N) offers the fraction of the total SNPs that showed both alleles within each population. The P_N was calculated as the proportion of SNPs with more than 1 % MAF within each breed. The MAF is the approximate frequency of the second most common allele per breed. The output from PLINK for observed and expected heterozygosity per animal within breed was subjected to further calculation to get an average estimate of H_o and H_e per breed. The heterozygosity values were calculated by getting the average of all SNPs (that is the sum of all heterozygosity values averaged over the total number of SNPs passed the quality control).

The population structure and relatedness were estimated by principal components analysis (PCA) using the R package SNPRELATE (Zheng et al., 2012) and admixture proportion inference using model-based clustering ADMIXTURE 1.3.0 software (Alexander, Novembre, & Lange, 2009). The PCA analysis allowed for visual investigation and solid quantitative summaries. The admixture analysis inferred the proportions of ancestry within the populations by use of prior defined K-values matching the assumed number of ancestral populations. The admixture procedure employs a maximum-likelihood based method by converging the ancestry proportions and allele frequencies that maximize the likelihood function. The most optimal population structure was determined by cross-validation error procedure (McVean, 2009) with assumed ADMIXTURE runs from K = 2 to K = 4. The K-value with the lowest CV error was selected as the optimal value. A phylogenetic tree based on Reynolds genetic distances representing relationships among goat breeds was visualized using iTOL software (Letunic & Bork, 2019).

Genetic diversity

The four indices of genetic diversity (polymorphic SNPs, mean allele frequency, observed and expected heterozygosity and inbreeding coefficient) were calculated within each breed (Table 2). The assessment of the proportion of SNPs that exhibited both alleles within each breed indicated high levels of diversity. The percentage of within-breed polymorphic SNPs ranged from 94.6% to 80.7%. The highest values of polymorphic loci were found in Alpine (94.6%) and Saanen (92.2%) while the lowest proportion was found in the Galla breed (80.7%). Across all loci, the lowest MAF was found in Galla (0.291) and the highest in Alpine (0.323).

Results revealed differences in genetic diversity between breeds. The expected heterozygosity was, in all cases, higher than the observed heterozygosity (H_e > H_o). The Alpine had the lowest observed heterozygosity (H_o = 0.558 ± 0.026) while Toggenburg had the highest (H_o = 0.580 ± 0.032). Inbreeding coefficients for all the breeds were negative and ranged between -0.013 (Toggenburg) and -0.042 (Galla).

Population structure analysis

Principal components analysis was used to cluster goats and explore the association among individuals and breed groups. In Figure 2, the principal component 1, which accounts for 15.2% of the total variance, separated Galla breed from the other three breeds. The principal component 2 accounts for 14.1% of the total variance, split the goat breeds into four clusters (Alpine, Saanen, Galla and Toggenburg clusters). One outlier, was, however, observed for the Saanen population mixing with Alpine population.

To examine admixture between the breeds, model-based clustering was performed and the most likely number of genetic population (cluster or K) between the four goat breeds were deduced using ADMIXTURE cross-validation procedure (McVean, 2009). The K-value with the lowest CV error was K = 4 and was selected as an optimal number of ancestral populations (Figure 3).

A population structure plot (Figure 4) showed proportions of ancestral populations for all breeds (Alpine, Galla, Saanen and Toggenburg) for K = 2 to K = 4. At K = 2, Galla goats were separated from the other three goat breeds (Toggenburg, Saanen and Alpine). Moreover, Galla goats largely do not carry ancestral components present in Saanen, Alpine and Toggenburg goats (shown in light blue, Figure 4). At K = 3, Alpine and Saanen goats carry ancestral components largely absent from either the Galla or Toggenburg goats. At K = 4, Galla had the lowest level of admixture, whereas Toggenburg, Alpine and Saanen demonstrated some signs of admixture with Galla.

The proportions of individuals in each breed in the four most likely clusters estimated by ADMIXTURE is shown in Table 3. 71 % of Alpine breed were allocated to cluster one, 97 % of Galla were assigned to cluster two with one percent (1 %) of its genome assigned to cluster one, three and four, 84 % percent of Saanen were in cluster three with seven percent (7 %) of its genome assigned to cluster one. On the other hand, 78 % of Toggenburg were assigned to cluster four with 17 %, and three percent (3 %) of its genome allocated to cluster two and one respectively.

Breed relationships were evaluated by computing the genetic distance between all pairwise combinations of individuals (D) from the average proportion of shared alleles. Based on the calculated Reynolds genetic distances, a phylogenetic tree was constructed to represent breed clustering (Figure 5). The results revealed five clusters for the four populations (Alpine, Saanen, Toggenburg and Galla). Some Toggenburg goats were found to be grouped together with Alpine, forming the fifth cluster.

Genetic diversity

Livestock has been exposed to various forces that contributed to the genetic diversity underlying phenotypic dissimilarities ever since domestication. These forces include natural selection, artificial selection for specific traits, migration, genetic drift and inbreeding (Andersson & Georges, 2004; Groeneveld et al., 2010). Genetic drift plays an important role during short-term evolution in situations where populations are reproductively isolated (Laval, SanCristobal, & Chevalet, 2002).

Genotyping with the GoatSNP50 Bead Chip revealed some levels of diversity within the goat breeds in this study. In each breed, fewer than 80% of SNPs exhibited polymorphisms, and heterozygosity ranged from 0.558 to 0.580 (Table 2). A large number of polymorphic SNPs were detected for Alpine (94.6) and Saanen (92.2) breeds; this was expected because sequenced data from Alpine and Saanen were included in the 50K SNPs panel discovery (Tosser-Klopp et al., 2014). Other results from various authors using different numbers of samples and goat breeds showed >93 % of polymorphism (Onzima et al., 2018; Visser, Lashmar, Marle-Köster, Poli, & Allain, 2016). However, it is difficult to compare and conclude the estimates of SNPs stated as polymorphic by other authors, because the number of samples genotyped per breed and proportion of genotyped samples used for SNP discovery varied.

The diversity amongst the four breeds showed Galla had the lowest diversity among sampled individuals in comparison with the other breeds. The polymorphic variance is dependent on the history of each breed. As opposed to Galla goats, the three introduced breeds showed greater diversity, likely as a result of crossbreeding with local goats (Galla and small East African). Therefore, each breed may contain genetic contributions from various breeds, thus revealing higher polymorphisms than Galla goats. In contrast, Galla goats sampled from the government breeding station with detailed pedigree records still maintained levels of diversity (P_N = 81%). The Galla population had been introduced in the Naivasha Sheep and Goat Station in the early 1970s during a sheep and goat project funded by FAO (Palian & Racokzi, 1976). The population is registered with the Galla Goat Breeders Society of Kenya (GGBSK) and the Kenya Stud Book. The population is inspected every year by inspectors from GGBSK using the Galla goat breed standards. Therefore, the station has maintained pure breed Galla goats which produce meat, milk and reproduce under harsh conditions while maintaining or conforming to the set standards of excellence defined by the GGBSK, where the objective is genetic improvement of target traits while controlling the level of inbreeding. The negative inbreeding coefficient in this study can effectively be taken as zero values, which means that there is no inbreeding observed in the reference populations. It could also mean that many heterozygotes were observed although the sample size for the four breeds was small. The increased heterozygosity could be due to random mating within the herd rather than random differences between herds.

The observed heterozygosity was lower than the expected heterozygosity (H_o < H_e) in the three breeds apart from Toggenburg, which recorded the same value for both observed and expected heterozygosity. The difference between the observed and expected heterozygosity was small, which may not be due to inbreeding but a Wahlund effect (Garnier-Géré & Chikhi, 2013). The observed heterozygosity in the current study was from a sample of individuals from a structured population even though all sub-divisions were in Hardy-Weinberg equilibrium. Over and above the semi and intensive production systems practised by smallholders, there is the presence of artificial selection, gene flow and non-random mating, hence not holding the law of HWE in these populations. In this study the observed and expected heterozygosity for Alpine (H_o = 0.558; H_e = 0.564), Saanen (H_o = 0.559; H_e = 0.573) and Toggenburg (H_o = 0.580; H_e = 0.580) were higher than those stated in Canada for Alpine (H_o = 0.385; H_e = 0.388), Saanen (H_o = 0.379; He = 0.382) and Toggenburg (H_o = 353; H_e = 336) (Brito et al., 2017). Moreover, Saanen in Italy recorded the same trend as in Canada (H_o = 0.41; H_e = 0.41) (Nicoloso et al., 2015). Differences in effective population sizes, length of isolation, selection and breeding management practices in the various production system may be the cause of these variances.

Toggenburg and Galla breeds had the highest expected heterozygosity. This could be explained by the types of crossbreeding programmes practised by farmers keeping Toggenburg breeds resulting in an admixed population. Organized breeding strategies using artificial selection are practised for Galla goats under the government breeding station resulting in genetic variability and lack of inbreeding for the populations.

Population structure and relationship

Principal component and population structure analyses confirmed distinctiveness among the goat breeds (Saanen, Galla, Toggenburg and Alpine) according to their geographic regions in Kenya. This can be explained by the demographic history of these breeds that have been reared for a long time in separate geographic locations (Ahuya et al., 2009; Peacock et al., 2011). Although goats from each breed clustered separately, model-based clustering revealed some signs of admixture and genetic links between Alpine, Toggenburg, Saanen and Galla.

The results (Figure 4 and Table 3) of this study indicate that Kenyan Alpine goats were the most admixed breed with about 14 % of its genome derived from Galla, while ten and six percent of its genome is resulting from Saanen and Toggenburg respectively. It is worth noting that Saanen were introduced in the sampling region (Nyeri County) already before the Alpine were imported in the late 1970s. Therefore, the 10 % of Saanen genes in the Alpine genome may be a result of Saanen being one of the Kenyan Alpine ancestors. According to Waineina, Okeno, Ngeno, and Ilatsia (2019), lack of breeding stock was one of the challenges Alpine farmers were encountering, thus driving them to source breeding animals from local markets, friends, neighbours and commercial farms notwithstanding their undefined genetic composition. Furthermore, the increase in demand for dairy goats in the country has resulted in several farms setting up nucleus flocks with a significant proportion of the crossbred flocks as a source of breeding material for distribution to lower cadre farmers (Ahuya et al., 2006; Bett, 2009; Ogola et al., 2010). Through such arrangements, most of the breed-types have migrated to other areas apart from their original entry in the country (Mburu et al., 2014; Peacock, 2007). Toggenburg and Alpine goats shared some linkage with Galla goats, 17 % and 14 %, respectively. This was expected because Galla goat was used as the founder population for crossbreeding with Alpine and Toggenburg breeds (Ahuya et al., 2009; Bett et al., 2011; Mburu et al., 2014; Peacock et al., 2011; Shivairo, Matofari, Muleke, Migwi, & Lugairi, 2013).

As expected, Galla was the least admixed breed, in agreement with the history of this breed as the first indigenous goat for which a breed society was formed in Kenya. Moreover, the particular population in this study has been managed in seclusion within the government farm, and only animals registered within the society are allowed into the population. The Galla breed displayed isolation by distance and seemed to be at equilibrium under dispersal and genetic drift. In comparison with the other breeds in this study, Galla arrived in their current locations long before these breeds were introduced in Kenya and that is why there has been sufficient time for isolation by distance to take effect and, that long distance gene dispersal is sufficiently common to prevent genetic divergence.

The phylogenetic analysis categorized the breeds into five clusters (Figure 5). The outcomes show a clear differentiation of Galla, Saanen, some Toggenburg and Alpine. A group of some Alpine and Toggenburg, however, remained clustered together, which may be attributed to the adjacent regions of the breeds (Figure 1). Lack of differentiation in some of Alpine and Toggenburg breeds signified a high level of genetic resemblance and low divergence, which may be a result of gene flow among Alpine and Toggenburg breeds. Common ancestry, short domestication history, lack of selection pressure and movement of the goats may play a role in lack of differentiation in varied geographically separated populations. Furthermore, in Kenya, as well as other parts of Africa, goats are also used for religious and other cultural ceremonies such as payment of dowry and gifts (Herrero et al., 2013; Mbuku et al., 2015; Ogola et al., 2010). Therefore, some of the Alpine and Toggenburg breeds clustering together may be a result of movement of breed animals between the communities in those two regions due to the forementioned cultural ceremonies. As mentioned earlier, one of the criteria for selecting the goat keeping households in this study was them being members of a dairy goat farmer group association (DGAK, MGBA, NGBA). The associations are responsible for buck rotation among the group members, maintaining the purtity of the breed and providing technical backstoping. However, the results indicate a need to technically strengthen the Dairy Goat Association of Kenya for Alpine and the Meru Goat Breeders Association for Toggenburg, because urgent management efforts are essential to improve on breeding aspects, utilization and conservation of the various goat genetic resources.

All Saanen goats formed one cluster in the phylogenetic analysis. Indeed, the long distance (over 450 km) between the regions where Saanen and the rest of the breeds are kept may be the barrier to gene flow from other breeds. Through adaptive hitchhiking, natural selection can play an essential role in shaping this variability (Andolfatto, 2001). Therefore, the observed genetic divergence of Saanen from Alpine, Toggenburg and Galla breeds could have been contributed by random genetic drift and natural selection for adaptation to their environment/region.

Genetic uniqueness can be determined from the magnitude of genetic distances and phylogenetic relationships between populations if supporting indications such as genetic history, records of production, reproduction and on adaptation are lacking (Eding & Laval, 1999; Tosser-Klopp et al., 2014; Zheng et al., 2012). Embracing this principle with respect to the results of this study, Saanen and Galla breeds seem to be the most genetically distinct among the populations sampled, and can be categorised as important genetic resources. It will be interesting to enlarge this breed level investigation in later studies through addition of all Kenyan goat breeds to better appreciate the genetic relationship among them.