Conservation status and historical relatedness of South African communal indigenous goat populations using a genome-wide single-nucleotide polymorphism marker.

Indigenous goats form the majority of populations in smallholder, low input, low output production systems and are considered an important genetic resource due to their adaptability to different production environments and support of communal farming. Effective population size (N e ), inbreeding levels, and the runs of homozygosity (ROHs) are effective tools for exploring the genetic diversity and understanding the demographic history in efforts to support breeding strategies to use and conserve genetic resources. Across populations, the current N e of Gauteng was the lowest at 371 animals, while the historical N e across populations suggests that the ancestor N e has decreased by 53.86%, 44.58%, 42.16%, and 41.16% in Free State (FS), North West (NW), Limpopo (LP), and Gauteng (GP), respectively, over the last 971 generations. Genomic inbreeding levels related to ancient kinship (F ROH > 5 Mb) were highest in FS (0.08 ± 0.09) and lowest in the Eastern Cape (EC) (0.02 ± 0.02). A total of 871 ROH island regions which include important environmental adaptation and hermo-tolerance genes such as IL10RB, IL23A, FGF9, IGF1, EGR1, MTOR, and MAPK3 were identified (occurring in over 20% of the samples) in FS (n = 37), GP (n = 42), and NW (n = 2) populations only. The mean length of ROH across populations was 7.76 Mb and ranged from 1.61 Mb in KwaZulu-Natal (KZN) to 98.05 Mb (GP and NW). The distribution of ROH according to their size showed that the majority (n = 1949) of the detected ROH were > 5 Mb in length compared to the other categories. Assuming two hypothetical ancestral populations, the populations from KZN and LP are revealed, supporting PC 1. The genomes of KZN and LP share a common origin but have substantial admixture from the EC and NW populations. The findings revealed that the occurrence of high N e and autozygosity varied largely across breeds in communal indigenous goat populations at recent and ancient events when a genome-wide single-nucleotide polymorphism (SNP) marker was used. The use of Illumina goat SNP50K BeadChip shows that there was a migration route of communal indigenous goat populations from the northern part (LP) of South Africa to the eastern areas of the KZN that confirmed their historical relatedness and coincides with the migration periods of the Bantu nation.

Introduction

There are currently 15 South African goat genetic resources listed on the DAD-IS of FAO and 13 on the Domestic Animal Genetic Resources Information System (DAGRIS) of the International Livestock Research Institute, including those listed in DAD-IS. In the country, indigenous goat ecotypes have been used as triple-purpose animals (e.g., skin, milk, and meat); and depending on the region, the animal characteristics, and the geographical isolation, they have begun to diverge into breeds/populations (Visser et al., 2016). These ecotypes have generally been named after their place of origin (e.g., Northern Cape Skilder) and/or their prominent characteristics (e.g., Xhosa lobbed ear) and on the basis of the people who own them (e.g., Nguni) (Morrison, 2007). These ecotypes are widely spread across all major agro-ecological regions of South Africa, displaying adaptability traits to a specific habitat or production environment and represent a significant resource to satisfy present and future demands for sustainable farming in a changing environment.

Improvement of indigenous livestock has been practiced through the introduction of high-performing breeds (exotic and improved breeds), and as a result of indiscriminate mating and breeding, the majority of communal indigenous goat populations are crossbreds (Mdladla et al., 2016a). The majority of the smallholder farmers have small herds or flocks where herd sizes could be less than five animals, with the majority of these herds being non-descript, crossbred or indigenous cattle, sheep, and goats (Mthi et al., 2017; Nyamushamba et al., 2017). The reduction in local indigenous populations suggests a need for the conservation of local genetic resources through the implementation of a national conservation strategy. Various studies conducted in the smallholder communal areas showed average flock sizes of between 1 and 120 goats (Webb and Mamabolo, 2004; Dube, 2015). Furthermore, about 70% of the goats are kept under traditional management systems where the farm structure comprises of about twenty (±20) indigenous goats (Monau et al., 2020). A detailed information on the phenotypic, genetic diversity and population structure of goats ecotype populations become important to guide conservation strategies through utilization of these populations. According to Food and Agriculture Organization of the United Nations, (2015b), conservation and characterization of animal genetic resources is critical because of their contribution to the sustainable livelihoods of rural communities that depend on them for food security. Conservation frameworks should incorporate both genetic diversity and breed merits for prioritizing breeds/populations from community to national level to support breeding programs of current populations.

In South Africa, more extensive research studies on genetic diversity analyses were done using microsatellite markers that were instrumental in providing an insight into the genetic structure and variation among South African goat populations (Kotze et al., 2004; Visser et al., 2004; Kotzé et al., 2014). (Kotzé et al., 2014) observed average heterozygosity of 63% in Kalahari Reds using 18 microsatellite markers, nine of which were used in the study by (Visser et al., 2004). Recently, microsatellites have been used to study genetic variation of the Tankwa feral goat population, which showed it to be highly divergent from the other farmed populations (Kotzé et al., 2014). In spite of their common use in most livestock diversity studies, microsatellites are often criticized for their usual location in the non-coding regions of the genome and for not being directly associated with genes that affect phenotypes. This has led to low-density microsatellites finding little application in studies of the adaptive genetic diversity of local breeds.

The completion of the first draft of the goat genome (Dong et al., 2013) made it possible for the development of high-density markers (Tosser-Klopp et al., 2014). The Illumina goat SNP50K BeadChip includes 53 347 SNPs (Tosser-Klopp et al., 2014) that have found utility in South African population genetic studies in Angora (Lashmar et al., 2016), commercial, indigenous, and village goat populations (Mdladla et al., 2017), as well as investigating genetic adaptation to environmental pressures (Mdladla et al., 2018). The use of the tool has been described in other African countries (Zvinorova, 2017; Onzima et al., 2018) and in specialized breeds (Martin et al., 2016). In South Africa, research work on the use of Illumina goat SNP50K BeadChip to determine the differences of indigenous communal goats at a point of genetic background is limited as compared to the studies where microsatellites were used. Therefore, while much work on South African commercial, indigenous, and village goat populations has focused on genetic studies and investigations on genetic adaptation, less work has focused on the conservation status and historical relatedness of communal indigenous goat populations.

The presence of the extent of an effective population size (Ne) is an important population genetic parameter that has recently received a great deal of research attention (Zhao et al., 2014), determining population demographic development (Deng et al., 2019) and demographic processes such as migration and admixture (Nicoloso et al., 2015), and having profound implications for understanding the architecture of the animal genome (Lashmar et al., 2016; Alvarenga et al., 2018). In addition, Ne is widely regarded as one of the most critical population parameters because it measures the rates of genetic drift and inbreeding as well as affects the efficacy of systematic evolutionary forces such as mutation, selection, and migration (Shin et al., 2013; Silió et al., 2013). It also helps to discover population demographic history and allows for the prediction of the behavior of genetic diversity through time. The Ne is estimated using the r2 coefficient and measures the observed range and the amount of genetic variation within a frame of population genetics (Berihulay et al., 2019). It also provides information on the degree of inbreeding of the population under consideration (Flury et al., 2010). The Ne determines the amount of genetic variation, genetic drift, and linkage disequilibrium (LD) in populations (Liu et al., 2017). Implementation of a national conservation strategy for FAnGR must be based on a better understanding of the degree of inbreeding of the populations, genetic variation, genetic drift, and linkage disequilibrium (LD) in those populations.

An increase in inbreeding (F) over generations leads to a reduction in genetic diversity (Onzima et al., 2018). Higher frequency of homozygous genotypes for deleterious alleles with a reduction in individual performance (inbreeding depression) and lower population viability (Ouborg et al., 2010). Offspring may inherit autozygotic chromosomal segments from both parents that are identical by descent (IBD) when they are inbred, i.e., segments that are derived from a common ancestor (Broman and Weber, 1999). The result is the runs of homozygosity (ROH), also known as the continuous homozygous segments in the genome. The ROHs are contiguous lengths of homozygous segments of the genome where the two haplotypes inherited from the parents are identical (Gibson et al., 2006).

The extent of ROH can be used to estimate the inbreeding coefficient (Bosse et al., 2012; Martin et al., 2016; Peripolli et al., 2018), to disclose the genetic relationships among individuals, usually estimating with high accuracy the autozygosity at an individual and/or population level (Ferenčaković et al., 2011; Ferenčaković et al., 2013). Autozygosity is the homozygous state of identical by descent (IBD) alleles, which can result from several phenomena such as genetic drift, population bottlenecks, mating of close relatives, and natural and artificial selection (Curik et al., 2014; Lashmar et al., 2016). The level of selection pressure on the populations can also be established by the use of ROH (Zhang et al., 2015). Distant from more recent inbreeding may also be distinguished by the length and frequency of ROHs, since the length of IBD segments follows an inverse exponential distribution with a mean of ½ g Morgans, where g is the number of generations from a common ancestor (Howrigan et al., 2011). Effective populations (Ne) and ROH have been studied in humans (Gibson et al., 2006), cattle (Ferenčaković et al., 2011; Marras et al., 2014; Mastrangelo et al., 2016), pigs (Ai et al., 2013; Shin et al., 2013), and sheep (Mastrangelo et al., 2017), but less comprehensively in other livestock species, such as goats, for designing conservation strategies, especially for the communal indigenous goats of South Africa. The objective of this study was to determine the conservation status and historical relatedness of South African communal indigenous goat populations using genome-wide SNP markers.

Methods

Ethics approval

Permission for the study and ethical approval were obtained from the animal research ethics committees of both the University of Limpopo and the Department of Agriculture, Land Reform and Rural Development. Furthermore, verbal consent was given by the goat owners.

Sample collection and genotyping

A total of 117 communal indigenous goat populations were sampled from the Free State (n = 24), Gauteng (n = 28), Limpopo (n = 30) and North West (n = 35) provinces of South Africa (Figure 1). Additionally, genotypes of communal indigenous goat populations that are under extensive production system (Mdladla, 2016) representing Eastern Cape (n = 20), KwaZulu-Natal (n = 30), Limpopo (n = 30) and North West (n = 20) provinces were included for further analysis.

FIGURE 1

South African map showing the geographical locations of communal indigenous goat populations.

All animals were genotyped using the Illumina GoatSNP50 BeadChip (Illumina Inc., San Diego, CA, United States) using the Infinium assay compatible with the Illumina HiScanSQ genotyping platform at the Agricultural Research Council—Biotechnology Platform in South Africa. A number of quality control measures were applied to all SNPs as follows: SNPs were removed if they had a call rate < 95%, a minor allele frequency (MAF) < 0.05, and if they deviated from the Hardy-Weinberg equilibrium (for a p-value cut-off of 0.00001), had no assigned genomic locations, and on sex chromosomes were also excluded from the analysis. The parameter “–cow” was used to indicate the number (i.e., 29) of autosomes in the goat genome since cows and goats have the same number of autosomes. The final dataset included 47,778 SNPs and 207 individuals.

Genetic diversity indices

Historical and recent effective population sizes (N ) for each breed were estimated with the SNeP (Barbato et al., 2015), which is based on the relationship between linkage disequilibrium (LD), N and recombination rate. The different SNP marker distance bins used for r analysis were used to obtain different estimates of Ne, at different time points by calculating the number of generations (t) in the past as ½c. To verify the accuracy of the coefficient of inbreeding, the genomic coefficient was estimated via two methods. 1) PLINK (Purcell et al., 2007) was used to measure the inbreeding coefficient based on the difference between the observed and expected numbers of homozygous genotypes (F ) using the function –het. The calculation formula was as follows: where L is the number of genotyped autosomal SNPs, E is the number of expected homozygous genotypes, and O is the number of observed homozygous genotypes. The inbreeding coefficient based on the proportion of autosomes covered in runs of homozygosity per individual (F ) was determined using detectRUNs (Biscarini et al., 2018). F was calculated as follows: where L is the total length of ROH on autosomes and L is the total length of the autosomes covered by SNPs, which was 2450.71 Mb. Furthermore, the correlation between F and F was calculated for the four populations.

Distribution of runs of homozygosity

Runs of homozygosity (ROH) were identified in every individual using detectRUNS (Biscarini et al., 2018) using a sliding window of a specified length or number of homozygous SNPs to scan along with each individual’s genotype at each SNP marker position to detect homozygous segments. The parameters and thresholds applied to define a ROH were: 1) a sliding window of 50 SNPs across the genome; 2) the minimum number of consecutive SNPs included in a ROH was 50; 3) the minimum length of a ROH was set to 1 Mb to avoid short and common ROH that occur throughout the genome due to LD (Purfield et al., 2017); and 4) a maximum of five SNPs with missing genotypes were allowed in a ROH to eliminate false positives. ROH was divided into five physical length categories (1–5 Mb, 5–10 Mb, 10–20 Mb, 20–30 Mb and < 40 Mb). The mean number of ROH per individual, the average length of ROH, the total number of ROH per animal, the percentage of chromosomes covered by ROH, and mean ROH were calculated on detectRUNS. The genomic inbreeding coefficient based on ROHs (F ) was also calculated as the sum of the length of the autosome covered by ROHs divided by the overall length of the autosome covered by SNPs as described by (Mastrangelo et al., 2017). The means and standard deviation (sd) of F were calculated as the sum of the lengths of F F F

To identify the genomic regions that were most commonly associated with ROH, the percentage of the occurrences of a SNP in ROH was calculated by counting the number of times the SNP was detected in those ROH across individuals, and this was plotted against the position of the SNP along the chromosome. This percentage threshold was normalized to 70%, 50%, and 20% of individuals per population to be an indication of a possible hotspot of ROH in the genome. The functions of these genes and pathways in which they are involved were assessed using the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) database and literature search.

Population structure

Principal component analysis (PCA) was calculated and plotted in Golden Helix SNP and Variation Suite (SVS) V8.1 (Golden Helix, Inc. 2012). For the analysis of ancestry proportions (admixture) with K set from 2 to 10, the ADMIXTURE v1.3.0 program (Alexander et al., 2009) was used. The default parameter of PLINK v1.9 (50 SNPs step 5 SNPs, r2 0.5) was used to subject the whole genotype data set to linkage disequilibrium (LD) pruning (Pemberton et al., 2012) prior to use in admixture analysis. Finally, to visualize admixture plots, GENESIS software was used (Buchmann and Hazelhurst, 2014).

Results

Variation of the estimated N at t generations ago (from 12 to 983) is presented in Figure 2. As expected, N decreased progressively across generations, however, N was higher than 150 for all breeds 12 generations ago. The variation in N across generations was smallest for Gauteng (N = 371) and the Free State (N = 386), whilst Limpopo had the highest (N = 723). Ancestral populations exhibited considerably larger N values, with the largest historical N values (N = 5,772).

FIGURE 2

Effective population size (Ne) for the communal indigenous goat populations from Free State, Gauteng, Limpopo and North West.

The average inbreeding coefficient (F ) was lowest in Free State (F = 0.03 ± 0.09), followed by 0.04 ± 0.09 in Gauteng, 0.05 ± 0.01 in the North West and the highest in Limpopo (F = 0.09 ± 0.05). Overall, the inbreeding level was 0.06 ± 0.08. The average F , its range of variation across populations and its distribution are summarized in Table 1. Genomic inbreeding coefficients (F ) based on the distribution of the length of runs of homozygosity by class are described in Table 2 and by chromosome in Figure 3. F differed significantly among populations across the length categories. The genomic inbreeding coefficients of the shortest ROH (0–5 Mb; related to ancient kinship) per animal ranged from 0.02 ± 0.02 in the Eastern Cape population to 0.08 ± 0.10 in the North West population. The F of Eastern Cape, Limpopo populations increased with category size, whilst decreased in Free State. Gauteng F decreased from 0.07 ± 0.09 to 0.05 ± 0.09 for F >20 Mb and increased at >20 Mb. In KwaZulu-Natal, F increased to up to 0.08 ± 0.09 at <20 Mb and decreased for > 20 Mb category.

TABLE 1

Distribution of runs of homozygosity inbreeding coefficients (F ) within each population.

Class (Mb)	Eastern cape (n = 20)	Free state (n = 24)	Gauteng (n = 28)	KwaZulu-natal (n = 25)	Limpopo (n = 55)	North west (n = 55)
0–5 Mb	0.02 ± 0.02	0.08 ± 0.09	0.07 ± 0.09	0.04 ± 0.07	0.03 ± 0.04	0.08 ± 0.10
5–10 Mb	0.02 ± 0.01	0.06 ± 0.09	0.06 ± 0.09	0.05 ± 0.08	0.03 ± 0.04	0.07 ± 0.09
10–20 Mb	0.02 ± 0.01	0.06 ± 0.08	0.05 ± 0.09	0.08 ± 0.09	0.03 ± 0.05	0.07 ± 0.09
>20 Mb	0.03 ± 0.01	0.06 ± 0.04	0.08 ± 0.12	0.06 ± 0.05	0.07 ± 0.03	0.08 ± 0.09

TABLE 2

Number of runs of homozygosity (nROH) and length (in Mb) categorized by ROH length class (ROH1–5 Mb, ROH5–10 Mb, ROH10–20 Mb, ROH20–40 Mb and ROH>40 Mb).

Class (Mb)	Eastern cape (n = 20)	Free state (n = 24)	Gauteng (n = 28)	KwaZulu-natal (n = 25)	Limpopo (n = 55)	North west (n = 55)
0–5 Mb	130	311	390	151	322	645
5–10 Mb	34	146	193	46	100	271
10–20 Mb	8	77	68	34	47	163
20–40 Mb	7	34	19	27	25	72
>40 Mb	1	6	12	6	8	30
Total	180	574	682	264	502	1181

FIGURE 3

Distribution of inbreeding coefficients (FROH) based on runs of homozygosity (ROH) for each chromosome.

Chromosomal distribution of inbreeding showed higher inbreeding levels in chromosome 15, in Gauteng chromosome 19, in KwaZulu-Natal in chromosome 19 followed by chromosome 16. In the North West, chromosomes 22 and 19 had the highest inbreeding coefficients compared to other chromosomes.

A total of 3383 ROH were identified across populations, although the frequency and length of ROH differed per population. Among the 207 animals, only 1 animal in the Eastern Cape population was lacking ROH, whilst 206 (99.52%) had at least one ROH longer than 1 Mb. The mean number of ROH per population was 24.36, 23.92, 21.47, 9.13, and 9 in Gauteng, Free State, North West, KwaZulu-Natal, Limpopo, and Eastern Cape, respectively (Figure 4).

FIGURE 4

The genomic length with mean runs of homozygosity (ROH) per individual and the number of ROH for individuals.

The mean ROH length was 7.76 Mb across the population. The longest segment (SNP position 39467151–137516937) was observed in the Gauteng and North West populations and was 98.05 Mb in length (1992 SNPs) found on chromosome 1. In KwaZulu-Natal, Free State and Eastern Cape, the longest segments were found on chromosome 9 (83.54 Mb; 1691 SNPs), chromosome 7 (73.33 Mb; 1449 SNPs) and chromosome 7 (51.81 Mb; 1034 SNPs), respectively. ROH shorter than 5 Mb predominated (n = 1949) across all populations (Table 1), accounting for 57.61% of all detected segments and differed per population. These short segments accounted for 72.22% of the Eastern Cape population, followed by Limpopo (n = 322; 64%), KwaZulu-Natal (n = 151; 57.20%), Gauteng (n = 390; 57.18%), North West (n = 645; 54.61%), and Free State (n = 311; 54.18%).

The relationship between the mean number of ROH and the length of the genome covered by ROH per individual varies considerably among animals and populations. The number of ROH per chromosome displayed a specific pattern with the larger numbers found for chromosome 1, 2, and 3, a number that tended to decrease with chromosome length, and the smallest number on chromosome (Keller et al., 2011) with 44 segments. Chromosome 1 had the highest number of ROH and the Eastern Cape did not have ROH on chromosomes 23 and 28 (Figure 5).

FIGURE 5

Number of runs of homozygosity (ROH) per chromosome identified across all populations.

The proportion of chromosomes covered by ROH is illustrated in Figure 6. Overall, the highest coverage by ROH was observed on chromosomes 1, 2, 6, and 8 at 0.37, 0.31, 0.38, and 0.33, respectively. For the Eastern Cape population, chromosomes 4 and 8 had 0.07%, whilst Free State, Limpopo and North West had the highest proportion on chromosome 1 with 0.06%. Gauteng and KwaZulu-Natal had the highest on chromosome 6 with 0.06 and 0.08%, respectively.

FIGURE 6

Frequency of runs of homozygosity (ROH) per chromosome per population.

To identify the genomic regions that were most commonly associated with ROH, the percentage of SNPs in ROH was assessed by analyzing the frequency of a SNP occurring in those ROH across different individuals (20%), and this was plotted against the position of the SNP along the chromosome (Figure 7). The threshold of 70% and 50% of the individuals did not yield ROH islands across populations, so 20% was used. Several genomic regions were identified that frequently appeared in ROH within individual animals (Table 3). We identified 58 ROH islands at the 20% threshold in the Free State (n = 28) and Gauteng (n = 29) provinces. No ROH islands were observed in the Eastern Cape, KwaZulu-Natal, North West, and Limpopo at the set thresholds. The ROH hotspot with the highest occurrences (SNPs = 149) in Gauteng was located on chr7 (7.69 Mb). Chromosome position, start and end position of ROH, ROH length and number of SNPs within the genomic regions of extended homozygosity are reported in Table 3 and Additional file 1, Supplementary Table S1. A total of 871 genes are inside the ROH islands, which include important environmental adaptation and hermo-tolerance genes such as IL10RB, IL23A, FGF9, IGF1, EGR1, MTOR and MAPK3. An additional file 1 shows the KEGG 292 pathways associated with the genes (Supplementary Table S2) that include Fructose and mannose metabolism, Starch and sucrose metabolism, Vitamin B6 metabolism, Protein export and Phototransduction.

FIGURE 7

Manhattan plot of occurrences (%) of a SNP in ROH across populations.

TABLE 3

Regions of the ROH islands at 20% across Gauteng and Free State populations by length.

Population	Position (chr: SNP 1: SNP 2)	nSNP	Length (Mb)
Gauteng	3:54056151:55070625	22	1.01
Gauteng	15:7897116:9088595	28	1.19
Gauteng	9:76995912:78321101	30	1.33
Gauteng	8:53874425:55257845	34	1.38
Gauteng	19:27892936:29379671	33	1.49
Gauteng	24:7056817:8823189	37	1.77
Gauteng	1:122842771:124683449	34	1.84
Gauteng	15:43908135:45815862	33	1.91
Gauteng	21:30374560:32314473	45	1.94
Gauteng	14:51115518:53088572	38	1.97
Gauteng	25:29189441:31200199	43	2.01
Gauteng	2:53997579:56229103	44	2.23
Gauteng	11:85024580:87384856	47	2.36
Gauteng	8:60260123:62640344	51	2.38
Gauteng	7:29457232:31898006	46	2.44
Gauteng	11:45829634:48502615	52	2.67
Gauteng	1:17079512:19890667	60	2.81
Gauteng	11:95518465:98340959	43	2.82
Gauteng	19:31681378:34691035	57	3.01
Gauteng	1:28400001:31438929	61	3.04
Gauteng	6:15586763:18769652	60	3.18
Gauteng	3:62095554:65301414	62	3.21
Gauteng	8:55933868:59160945	67	3.23
Gauteng	11:50602742:54104265	72	3.50
Gauteng	15:9294238:12900069	70	3.61
Gauteng	22:35705179:39950926	81	4.25
Gauteng	7:45438722:51833004	126	6.39
Gauteng	7:60187788:67883046	149	7.70
Freestate	25:28641340:31924694	69	3.28
Freestate	7:48193142:49278113	19	1.08
Freestate	20:21671108:22811848	22	1.14
Freestate	15:30999298:32153016	30	1.15
Freestate	20:4743026:6211601	32	1.47
Freestate	6:11346351:12878262	32	1.53
Freestate	20:19597454:21199311	37	1.60
Freestate	14:73908424:75554996	30	1.65
Freestate	22:36144242:37987490	36	1.84
Freestate	6:107507640:109477575	40	1.97
Freestate	25:25914911:28121785	48	2.21
Freestate	10:73412812:75622539	47	2.21
Freestate	4:94620783:96868758	55	2.25
Freestate	10:70170521:72471193	48	2.30
Freestate	1:81289909:83642225	50	2.35
Freestate	18:6154744:8563747	54	2.41
Freestate	1:145037345:147663352	55	2.63
Freestate	8:62328607:64991613	62	2.66
Freestate	11:36630830:39372620	53	2.74
Freestate	17:32672670:35522389	54	2.85
Freestate	8:94585706:97590188	52	3.00
Freestate	12:57673212:60742379	61	3.07
Freestate	16:5685888:9117064	69	3.43
Freestate	10:48160488:51636811	64	3.48
Freestate	18:20922080:24775492	76	3.85
Freestate	12:63521630:67564200	87	4.04
Freestate	10:36878056:43333293	124	6.46
Freestate	3:33683337:40785725	140	7.10

For population structure analysis, further quality control parameters were effected in PLINK v1.90 (Peripolli et al., 2018): to remove linked SNPs using the --indep-pairwise 50 5 0.5 command and related individuals (IBS > 0.65). The database for population structure included 32886 SNPs and 189 individuals across populations. The plot of the first (PC1) and second eigenvectors (PC2) (Figure 8) shows weak differentiation among the clusters of admixed populations. PC1 shows clusters 1 and 2 consisting of populations from KwaZulu-Natal and Limpopo, respectively. The Eastern Cape population (Cluster 3) revealed a separate cluster (with three outliers). Clusters 4 consists of the majority of the North West and Gauteng, whilst cluster 5, was Free State.

FIGURE 8

A plot of principal components (PCA) analysis showing differentiations among the clusters of admixtures across six communal indigenous goat populations.

For the further understanding of the degree of admixture within the populations, the ADMIXTURE 1.3 (Alexander et al., 2009) software was used for K = 2 to 10 hypothetical ancestral populations, only K = 2 to K = 6 is plotted since no further clusters were observed after K = 6 (Figure 9). Assuming two hypothetical ancestral populations, the populations from KwaZulu-Natal and Limpopo are revealed, supporting PC 1. The genomes of KwaZulu-Natal and Limpopo share an origin but have substantial admixture from the Eastern Cape and North West populations. K = 3 as the most likely number of genetically distinct groups within our populations, presenting the lowest cross-validation error (0.4617), reflecting a clear cluster of the Limpopo populations. The Free State, Gauteng, and North West showed similar genetic heterogeneity patterns with a considerable level of admixture. The North West revealed a high level of within-population genetic differentiation as there are individuals closer to the Eastern Cape and another subpopulation closer or clustering with Free State. This is also in agreement with the second PC coordinate analysis in showing genetic heterogeneity within the population. Moreover, with the increment of the value of K (K = 4 to K = 6), Free State and Gauteng show a higher level of genetic heterogeneity than the other populations.

FIGURE 9

Clusters inferred from ADMIXTURE at K = 2–6.

Discussion

In this study, the genetic diversity and population structure of goat populations in South Africa were revealed. The main links to the development of commercial goat populations are represented by indigenous goat populations and may potentially be relevant as a future source of untapped adaptable genetic material (Webb and Mamabolo, 2004). Therefore, improving our understanding of within-ecotype relationships among the major goat producing provinces in South Africa offers a rare opportunity to enhance efficient use of the breeds and implement conservation programs. This study investigated the indications for population status on inbreeding and runs of homozygosity in the indigenous goat population. Data from a previous study (Mdladla et al., 2018) enabled a broad geographical coverage of South Africa and represented populations from the major goat producing provinces within the country.

Effective population size (N ) is a crucial population genetic parameter because of its relationship to the loss of genetic variation, increases in inbreeding, the accumulation of mutations, and the determination of the accuracy of genomic selection (Goddard, 2009; Berihulay et al., 2019). Gauteng had the smallest estimated N among the population, and Limpopo had the highest. It was also observed in other studies that effective population size (N ) showed a reduction to 132 in the Kingdom of Eswatini and the highest in South Africa 12 generations ago (Webb and Mamabolo, 2004). It is recommended that to prevent a reduction of the adaptive value in populations, N values between 50 and 100 animals (Meuwissen and Goddard, 2007) and 50 to avoid inbreeding depression. A study by (Colli et al., 2018), reported a large N in local goat breeds, such as those from Africa, Spain, and Central-Southern Italy, local goats breeds and a small Ne in the Angora, Boer, Nubian, Cashmere, Saanen, and Alpine populations.

The rapid increase pattern in Ne may also include bottlenecks associated with domestication, selection, and breed formation, as well as the endangerment of the breed (Shin et al., 2013). A study by (Mdladla, 2016) based on SNP data and using the same method, reported large N for all investigated breeds (ranging from 140 to 348). Furthermore, a study by (Mdladla et al., 2016a) revealed that the ecotype goat had a slightly higher effective population size than the Tankwa and commercial goat populations across generations. From a conservation standpoint, the indigenous goat population should be top priority in the population studied due to their diminishing effective population size and increased inbreeding coefficients.

Runs of homozygosity (ROH) can disclose the genetic relationships among individuals, estimating with high accuracy the autozygosity at the individual and population levels and can elucidate about selection pressure events (Purfield et al., 2012). If long ROH accumulates in the genomes of some individuals, they could seriously impact the overall biological fitness (Manunza et al., 2016). Therefore, it was an important objective to investigate and understand the level of homozygosity among the populations. In this study, only 1 animal was lacking ROH, whilst 206 (99.52%) had at least one ROH longer than 1 Mb. The genomic inbreeding coefficients (F ) values found in the study for the Ethiopian goats were F values (Nandolo et al., 2019). Similar results were found by (Nandolo et al., 2019) with more African goats (Cameroon, Ethiopia, Kenya, Madagascar, Malawi, Mali, Mozambique, Nigeria, Tanzania, Uganda, and Zimbabwe) using a clustering algorithm. In this study, differences in terms of total number and length of ROH were short (>5 Mb) were more abundant (57.61%) (Onzima et al., 2018) reported results that showed lower than the F for Kenya, Uganda, and Mozambique goat breeds when the Goat 50K BeadChip was used. On the other hand, the Eastern Cape population showed very low amounts of ROH. This has been suggested to be consistent with recent admixture in the individuals of Chinese cattle (Xu et al., 2019). Long segments were abundant in the North West population (Islam et al., 2019). A recent study revealed a high mean ROH in the long length category (>30 Mb), and their study suggested that inbreeding is more recent and is indicative of demographic decline. Considering the extensive management systems of goats in the region, these results might be likely even if some researchers have argued that such extensive systems may lead to inbreeding, especially where goats are on an extensive system or are shepherded with other flocks for some part of the year (Tefera et al., 2004; Rumosa Gwaze et al., 2009b). The lower inbreeding levels in African goats could be due to the openness of the breeding systems in most of Africa that leads to the loose definition of livestock breeds in the region (Manunza et al., 2016).

One of the main advantages of genomic coefficients is the availability of chromosomal inbreeding coefficients (McQuillan et al., 2008; Mastrangelo et al., 2016). ROHs, representing the level of genomic autozygosity, are continuous homozygous segments at the individual and population levels that can be used as a measurement of inbreeding; more in-depth ROHs are the result of demography, natural and artificial selection, and inbreeding (Shi and Manley, 2007). In this study, we do not discuss in detail all the genomic regions associated with ROH but focus on some selected regions that show associations with several specific traits related to livestock breeding. We identified five genes reported to be associated with the important traits of goats (Figure 6) identified by the selection signature. Overall, the highest coverage by ROH was observed on chromosomes 1, 2,and 6, respectively. Gene INHA, located on chromosome 2, was reported as a candidate gene for litter size in goats (Hou et al., 2012). Significant QTL for milk production traits such as milk yield and milk protein have been reported on chromosome 2 in sheep (Garcia-Gámez et al., 2013). Gene INHA, located on chromosome 2, was reported as a candidate gene for litter size in goats. The PPP1R36 and Heat Shock Protein A2 (HSPA2) (CHI10, 26.402–26.719 Mb) identified in these communal indigenous goats are involved in heat stress response and, in other studies, HSPA2, DNAJC24, and DNAJC13 are associated with the heat shock family of genes (Shi and Manley, 2007). The presence of multiple genes associated with heat stress would seem to suggest that the trait is under intense selection pressure in tropically adapted breeds (Onzima et al., 2018).

In accordance with (Nothnagel et al., 2010), these regions in humans, when they are present in more than 50% of the individuals of a population, can indicate a strong selection occurrence. The occurrence of ROH hotspots in genomic regions that harbor candidate genes may be involved in selection pressure in response to production and environmental conditions. This study identified 58 ROH hotspots in Gauteng and Free State populations and revealed 871 genes and 292 KEGG pathways. This threshold did not yield any results, and only about 20% of ROH islands were detected in the communal indigenous goat population. The 20% threshold has been used in indigenous Chinese pigs, Jinhua (Xu et al., 2019). ROH islands can be defined as genomic regions with reduced genetic diversity and, consequently, high homozygosity around the selected locus that might harbor targets of positive selection and are under strong selective pressure (Pemberton et al., 2012; Peripolli et al., 2018). The minimum expected length of homozygous DNA segments is based on the time frame of approximately 25 generations, over which goats are believed to have been characterized into separate breeds (Onzima et al., 2018).

This study explored the population genetic structure of the indigenous goat population in the context of all South African goat populations. In accordance with our earlier studies (Lashmar et al., 2016; Mdladla et al., 2017), the principal component analysis (PCA), and the ADMIXTURE analyses based on the SNP array and sequence data sets, capitulated the major genetic division among the South African goat populations from two large geographic regions: Eastern Cape and Limpopo. The KwaZulu-Natal and Limpopo populations’ genomes shared an origin yet with significant admixture from the Eastern Cape and North West populations. Some signals of admixture and underlying genetic relationships among the populations were generated by analysis of population admixture (Hou et al., 2012). A migration route of ancient goat from the northern part of South Africa to the eastern areas of the KZN, during their migration periods of the Bantu nation, is supported by this study as observed.

Conclusion

The results of this study indicated a greater negative impact of inbreeding in recent times, which is important for planning conservation strategies. It was revealed that the occurrence of high N and autozygosity varied largely across ecotypes in communal indigenous goat populations at recent and ancient events when a genome-wide SNP marker was used. The use of Illumina goat SNP50K BeadChip shows that there was a migration route of communal indigenous goat populations from the northern part (Limpopo province) of South Africa to the eastern areas of KwaZulu-Natal, which confirmed their historical relatedness and coincides with the migration periods of the Bantu nation. The communal-traditional indigenous goat farming system and adaptation to different climatic conditions had an influence on the results in this study. The study deepened the understanding of the conservation status and selection mechanisms of goats in a communal indigenous goat production setting. For effective conservation programs and utilization of South African communal indigenous goat populations, effort should be made to establish a conservation program for the unique genetic resources of indigenous goat populations.