Genetic Diversity of South African Indigenous Goat Population from Four Provinces Using Genome-Wide SNP Data

Genome-wide assessments of the genetic landscape of Farm Animal Genetic Resources (FAnGR) are key to developing sustainable breed improvements. Understanding the FAnGR adaptation to different environments and supporting their conservation programs from community initiative to national policymakers is very important. The objective of the study was to investigate the genetic diversity and population structure of communal indigenous goat populations from four provinces of South Africa. Communal indigenous goat populations from the Free State (FS) (n = 24), Gauteng (GP) (n = 28), Limpopo (LP) (n = 30), and North West (NW) (n = 35) provinces were genotyped using the Illumina Goats SNP50 BeadChip. An Illumina Goats SNP50 BeadChip data from commercial meat-type breeds: Boer (n = 33), Kalahari Red (n = 40), and Savanna (n = 31) was used in this study as reference populations. The Ho revealed that the genetic diversity of a population ranged between 0.39 ± 0.11 Ho in LP to 0.42 ± 0.09 Ho in NW. Analysis of molecular variance revealed variations of 3.39% (p < 0.0001) and 90.64% among and within populations, respectively. The first two Principal Component Analyses (PCAs) revealed a unique Limpopo population separated from GP, FS, and NW communal indigenous goat populations with high levels of admixture with commercial goat populations. There were unique populations of Kalahari and Savanna that were observed and admixed individuals. Marker FST (Limpopo versus commercial goat populations) revealed 442 outlier single nucleotide polymorphisms (SNPs) across all chromosomes, and the SNP with the highest FST value (FST = 0.72; chromosome 8) was located on the UHRF2 gene. Population differentiation tests (PCAdapt) revealed PC2 as optimal and five outlier SNPs were detected on chromosomes 10, 15, 20, and 21. The study revealed that the SNPs identified by the first two principal components show high FST values in LP communal goat populations and allowed us to identify candidate genes which can be used in the development of breed selection programs to improve this unique LP population and other communal goat population of FS, GP, and NW, and find genetic factors contributing to the adaptation to harsh environments. Effective management and utilization of South African communal indigenous goat populations is important, and effort should be made to maintain unique genetic resources for conservation.

genotyping platforms [14]. Furthermore, the recent developments in genomic technologies resulting in the availability of the Illumina Goat 50K SNP BeadChip [15] have presented the opportunity to search for genomic regions that may have undergone selection. Such studies in cattle [16,17], sheep [18,19] chickens [20], and pigs [21,22] have each identified genes that have undertaken a positive selection and are likely to contribute directly to phenotypic variation. Furthermore, genetic diversity studies on local and international goat breeds, e.g., Italian [23], Moroccan [24], and South African [20,25], were conducted. Additionally, with regard to goat genetic diversity, more studies on goat 52 K single nucleotide polymorphism chip has been available for some time [15], and has already been used for national and multi-country goat diversity studies [7,[26][27][28]. Using this standardized tool for genotyping, the AdaptMap initiative [11] gathered a dataset that includes genotypes of 148 goat populations which successfully investigated the distribution of goat genetic variation around the world, compared the present-day diversity patterns to those observed in ancient goat samples and in other ruminant livestock species, as well as reconstructed admixture and migration events that have shaped goat post-domestication history.
Studies on population structure and genetic diversity are important for describing the natural selection history and genetic relationships among FAnGR [29]. The genome-wide assessments of the genetic landscape of FAnGRs' is a key aspect of developing sustainable breed improvement strategies and understanding adaptations to extreme environments [30,31]. Additionally, it has value in supporting conservation programs from community initiatives to national policy makers. To facilitate the rapid adaptation to changing environments, maintaining local genetic resources is very much required. The objective of this study was to investigate the population structure and genetic diversity of communal indigenous goat populations from four provinces of South Africa.

Population Description and Sampling
A total of 117 Communal indigenous goat populations were systematically randomly sampled from the Free State (n = 24), Gauteng (n = 28), Limpopo (n = 30), and North West (n = 35) provinces of South Africa (Figure 1). The sampling strategy targeted the major goat producing provinces of South Africa (Figure 1). From each animal, blood samples were collected by jugular venipuncture into 6-mL EDTA vacutainer tubes (Greiner Bio-One, GmbH, Kremsmünster, Austria). The DNeasy Blood and Tissue Kit (Qiagen) was used to extract Genomic DNA, as per the manufacturer's instructions. A 1% agarose gel electrophoresis and Qubit 3.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) were used to determined DNA quality and quantity, respectively. Descriptions of sampling locations and farms have been described in detail in [31].

Genotyping and Quality Control
The Illumina Goats SNP50 BeadChip [15] was used to genotype DNA samples by using the Infinium assay compatible with the Illumina HiScan SQ genotyping platform at the Agricultural Research Council-Biotechnology Platform, South Africa. The BeadChip was developed by the International Goat Genome Consortium (IGGC) with more than 50K SNPs across the whole genome with an inter-SNP spacing of approximately 40 kb [15]. The genotype input file was converted into a PLINK and SNP calling rate using the Illumina Genome Studio v2.0 [32] map/ped input file. Information on the chromosome number, chromosomal position, and Golden Helix SVS v8.3.4 (Golden Helix, Bozeman, MT, USA) were used to update the SNP marker information. PLINK was used to perform individual and SNP quality control procedures [32]. The mind function procedure was used for analyses of individuals with a missing genotype call rate of greater than 5%, which were excluded for further analysis. The three samples that did not meet the quality control threshold (95%) included 2 individuals from Limpopo and 1 from Gauteng, and they were removed from the final dataset. The removal of SNPs with less than 95% call rate (n = 683), and less than 5% minor allele frequency (n = 1228) and HWE exact test (n = 123) resulted in a total of 47,908 autosomal SNPs available for downstream analysis. The average genotyping rate in the remaining 117 samples was 0.99. Additionally, to evaluate population structure and relatedness, commercial populations represented by Boer (n = 33), Kalahari Red (n = 40), and Savanna (n = 31) reference population data obtained from a previous study [33] were included. To prepare dataset 2, the same QC criterion was used, resulting 46,081 SNPs and 195 individuals. In addition, a genome-wide identity-by-descent (IBD) similarity matrix was used to calculate and prune related individuals and the '-indep-pairwise 50 5 0.2-' command was used to remove one of every pair of SNPs with r 2 > 0.2 within 50-SNP sliding windows in PLINK [32] for dataset 3. From the IBD threshold of greater than 0.40, denoting substantial relatedness, twenty-seven animals were removed, therefore data set 3 consisted of 20,873 SNPs and 168 individuals.

Genetic Diversity
Minor allele frequency (MAF) was estimated for each of the four communal indigenous goat populations using PLINK [32]. The fixed alleles (MAF = 0.00), rare alleles (>0.00-<0.05), intermediate alleles (≥0.05-<0.10), and common alleles (≥0.10-≤0.5) were categories for alleles that were used in different bins based on their frequency. PLINK was used to estimate within-population genetic diversity, observed heterozygosity (H o ), and the expected heterozygosity (H E ) [32]. The observed heterozygosity (H o ) estimates for each population were calculated from observed genotype frequencies obtained from PLINK [32] as follows: (N-O)/N (where N is the number of non-missing genotypes and O is the number of observed homozygous genotypes for a given individual). Expected heterozygosity (H E ) estimates for each population were calculated from expected genotype frequencies as follows: (N-E)/N (where N is the number of non-missing genotypes and E is the number of expected homozygous genotypes for a given individual). Allelic richness (R t ) over all loci for each population was computed using ADZE Allelic Diversity Analyzer: Version 1.0.

Population Structure Analysis
AMOVA was executed using dataset 2 and analyzed using ARLEQUIN 3.5 after conversion map/ped file conversion to arp file using PGDSpider. Two groups: The commercial and communal indigenous goat populations were created by doing the significance of variance components for each hierarchical comparison (among populations, among individuals, among individuals within populations). The pairwise fixation indexes (F ST ) [34] were used as a measure of genetic differentiation between population pairs using Golden Helix SVS v8.3.4 (Golden Helix, Bozeman, MT, USA). The Principal Component Analysis (PCA) was carried out using Golden Helix SVS v8.3.4 (Golden Helix, Bozeman, MT, USA), and to infer the most probable number of ancestral populations, ADMIXTURE 1.21 was used [35] using default settings. ADMIXTURE was run from K = 2 to K = 5. To explore the most Sustainability 2020, 12, 10361 5 of 16 probable number of clusters (K), ten-fold cross-validation (CV = 10) was specified, with the error profile obtained thereafter used as described by [35]. The GENESIS software was used for graphical display of the admixture output and was also used to visualize admixture plots [36].

Identification of Signatures of Selection
Identification of signatures of divergent selection between populations using population differentiation (Fst) was calculated using Golden Helix SVS v8.3.4 (Golden Helix, Bozeman, MT, USA) between the Limpopo population versus commercial goat populations: Boer, Savanna, and Kalahari Red. Outliers were identified, and a region was considered to be a high-F ST outlier if it corresponded to the upper 1% of the empirical genome-wide distribution of F ST . The second approach used PCAdapt [37] to find outlier SNPs. We first assessed the optimal Principal Components (PC) from 1 to 10 and then the candidate SNPs significantly correlating to the optimal PCs. All outlier SNPs detected using both methods were mapped to gene-associated regions based on the goat genome annotation (ARS 1) on the BioMart database (www.ensembl.org/biomart). The genes in regions with evidence for selection were searched against the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) pathways and literature search.

Quality Control and Genetic Diversity
The mean minor allele frequencies (MAFs) for Gauteng, Free State, Limpopo, and North West were 0.33 ± 0.11, 0.33 ± 0.11, 0.29 ± 0.12, and 0.33 ± 0.11, respectively, with an overall mean of 0.32 ± 0.12 across the Communal indigenous populations. Figure 2 is presenting minor allele frequency distribution for different categories in each population. The number of fixed SNPs (MAF = 0.00) varied from 5 in North West populations to 151 in Limpopo. Polymorphic SNPs were high (<90%) across all populations. The allelic richness for Gauteng, Free State, Limpopo, and North West were 1.99 ± 0.02, 1.99 ± 0.02, 1.99 ± 0.05, and 1.99 ± 0.01, respectively, and with the mean number of alleles per locus being 1.99 ± 0.03 across the communal indigenous populations.

Source of Variation and Population Structure
Analysis of AMOVA revealed variations of 3.39% (p < 0.0001) and 90.64% among and within populations, respectively (Table 1). The PCA of the Communal indigenous goat and commercial goat populations is presented in Figure 3. The analysis ignores population membership, but reveals clear population structures as samples from the same population cluster together. The first two PCs, separated Limpopo population (Cluster 1), and Cluster 2, consisting of Gauteng, Free State, and North West populations, whilst Cluster 3, Cluster 4, and Cluster 5 consisted of commercial goat populations. Additionally, outliers were evident for the Cluster 2 populations in Cluster 3. The Communal indigenous goat populations provided weak sub-structuring with an overlap of populations from Free State, Gauteng, and North West.

Source of Variation and Population Structure
Analysis of AMOVA revealed variations of 3.39% (p < 0.0001) and 90.64% among and within populations, respectively (Table 1). The PCA of the Communal indigenous goat and commercial goat populations is presented in Figure 3. The analysis ignores population membership, but reveals clear population structures as samples from the same population cluster together. The first two PCs, separated Limpopo population (Cluster 1), and Cluster 2, consisting of Gauteng, Free State, and North West populations, whilst Cluster 3, Cluster 4, and Cluster 5 consisted of commercial goat populations. Additionally, outliers were evident for the Cluster 2 populations in Cluster 3. The Communal indigenous goat populations provided weak sub-structuring with an overlap of populations from Free State, Gauteng, and North West. In Figure 4, ADMIXTURE was used to further understand the degree of admixture in the populations and increased K from 2 to 5, where K is the assumed number of ancestral populations. K = 2, Limpopo separated and a subset of Kalahari Red was revealed. K = 3 revealed two sub-groups of the Savanna goats and Boer that shared genetic background with the Savanna. The analysis suggested K = 3 (CV error = 0.65048) as the most likely number of genetically distinct groups, reflecting three In Figure 4, ADMIXTURE was used to further understand the degree of admixture in the populations and increased K from 2 to 5, where K is the assumed number of ancestral populations.   The pairwise mean FST supports the population relationships based on PCA and ADMIXTURE. High FST (FST = 0.084) was shown between the Limpopo population and commercial goat populations (Boer, Kalahari Red, and Savanna). A low FST value (range 0 to 0.038) within the Communal indigenous goat populations was found, indicating a low genetic differentiation between these populations as shown in Figure 6 and additional file 1, Table S2.    The pairwise mean FST supports the population relationships based on PCA and ADMIXTURE. High FST (FST = 0.084) was shown between the Limpopo population and commercial goat populations (Boer, Kalahari Red, and Savanna). A low FST value (range 0 to 0.038) within the Communal indigenous goat populations was found, indicating a low genetic differentiation between these populations as shown in Figure 6 and additional file 1, Table S2.  The pairwise mean F ST supports the population relationships based on PCA and ADMIXTURE. High F ST (F ST = 0.084) was shown between the Limpopo population and commercial goat populations (Boer, Kalahari Red, and Savanna). A low F ST value (range 0 to 0.038) within the Communal indigenous goat populations was found, indicating a low genetic differentiation between these populations as shown in Figure 6 and additional file 1, Table S2.

Candidate Genes under Selection
A total of 442 differentiated SNPs with an FST threshold ≥0.47 were considered to be under selection. The distribution of FST across all chromosomes is shown in Figure 7. Outlier SNPs were found across all chromosomes with varying numbers. In general, bigger chromosomes had a higher number than the smaller chromosomes. Chromosome 8 had the highest number with 37 SNPs, and chromosomes 28 had 4 SNPs.  Table 2 shows the top ten SNP and additional file 1, Table S1 has the full list of outlier SNPs. The SNP with the highest FST value (FST 0.72) was at chromosome 8 and sits on the UHRF2 gene. The 133 candidate genes (Additional file 1, Table S2) were adjacent to these outlier's loci that associated significantly with pathways involved in the immune system, thermoregulation, and longevity (additional file 1, Table S3). The other four important genes from the top ten outlier SNPs that were found include with the FST value (FST 0.69) at chromosome 8 which sits on the GLDC gene, FST value

Candidate Genes under Selection
A total of 442 differentiated SNPs with an F ST threshold ≥0.47 were considered to be under selection. The distribution of F ST across all chromosomes is shown in Figure 7. Outlier SNPs were found across all chromosomes with varying numbers. In general, bigger chromosomes had a higher number than the smaller chromosomes. Chromosome 8 had the highest number with 37 SNPs, and chromosomes 28 had 4 SNPs. Table 2 shows the top ten SNP and additional file 1, Table S1 has the full list of outlier SNPs. The SNP with the highest F ST value (F ST 0.72) was at chromosome 8 and sits on the UHRF2 gene. The 133 candidate genes (Additional file 1, Table S2) were adjacent to these outlier's loci that associated significantly with pathways involved in the immune system, thermoregulation, and longevity (additional file 1, Table S3). The other four important genes from the top ten outlier SNPs that were found include with the

Candidate Genes under Selection
A total of 442 differentiated SNPs with an FST threshold ≥0.47 were considered to be under selection. The distribution of FST across all chromosomes is shown in Figure 7. Outlier SNPs were found across all chromosomes with varying numbers. In general, bigger chromosomes had a higher number than the smaller chromosomes. Chromosome 8 had the highest number with 37 SNPs, and chromosomes 28 had 4 SNPs.  Table 2 shows the top ten SNP and additional file 1, Table S1 has the full list of outlier SNPs. The SNP with the highest FST value (FST 0.72) was at chromosome 8 and sits on the UHRF2 gene. The 133 candidate genes (Additional file 1, Table S2) were adjacent to these outlier's loci that associated significantly with pathways involved in the immune system, thermoregulation, and longevity (additional file 1, Table S3). The other four important genes from the top ten outlier SNPs that were found include with the FST value (FST 0.69) at chromosome 8 which sits on the GLDC gene, FST value  Outlier detection using principal components (PCs) shows the decay of eigenvalues confirmed to the use of K = 2 as optimal (Figure 8a). The Q-Q plot confirmed the expected uniform distribution of the p-values (Figure 8b). Figure 8c shows a Manhattan plot indicating the outlier SNPs that have been detected based on the expected FDR equal to 10%.
The results show that twenty-five outliers were significant to each PC, however, only 5 could be retained (Table 3). The 5 SNPs were found on chromosomes 10, 15, 20, and 21, and were associated with the PC2. In this study, the gene search did not reveal any genes.   Outlier detection using principal components (PCs) shows the decay of eigenvalues confirmed to the use of K = 2 as optimal (Figure 8a). The Q-Q plot confirmed the expected uniform distribution of the p-values (Figure 8b). Figure 8c shows a Manhattan plot indicating the outlier SNPs that have been detected based on the expected FDR equal to 10%. The results show that twenty-five outliers were significant to each PC, however, only 5 could be retained (Table 3). The 5 SNPs were found on chromosomes 10, 15, 20, and 21, and were associated with the PC2. In this study, the gene search did not reveal any genes.

Discussion
The ability to adapt to a changing environment is the genetic diversity of indigenous FAnGR, as they provide information to support improvement and conservation efforts. To design genetic improvement and conservation strategies at national and global levels, information on within and between genetic diversity and population structure becomes paramount [38]. In this study, population structure and genetic diversity of communal indigenous goat populations from the goat producing provinces of Free State (FS), Gauteng (GP), Limpopo (LP), and North West (NW) provinces were investigated.
Communal indigenous goat populations are genetically diverse, exclusive, and are adapted to particular production systems, and their value of diversity has been strongly emphasized [33,39]. The overall mean of the minor allele frequency (MAF) was 0.32 ± 0.12 across the communal indigenous goat populations. It has been indicated that the communal indigenous goat populations were more diverse than the commercial goat populations with regards to the observed within-population diversity measures. These results support the findings by [33] on the Communal indigenous goat population, where goat ecotypes of Venda and Zulu showed the maximum number of alleles per locus.
The number of fixed SNPs (MAF = 0.00) varied from five in NW populations to 151 in LP populations. Polymorphic SNPs were high (<90%) for all populations. The mean MAFs were estimated on 49,942 SNPs and 117 (<0.95 call rate) individuals for FS, GP, LP, and NW, respectively, with an overall mean of 0.32 ± 0.12 across the communal indigenous goat populations. The results are comparable to the results from [25], who recorded 43,759 (82.03%) observed for the South African Angora breed and [31] observing the overall proportion of informative SNPs (87.09%). In this study, a level of polymorphic SNPs was at a 95% call rate threshold, which was similar to the study of [33], whilst in [25]'s study, it was a 98% call rate threshold, which was higher for the Angora goat breed. The importance of the differences in the SNP quality control thresholds has to be emphasized, as it has an influence on the resultant proportion of useful SNPs.
Analysis of the distributions of alleles across populations is important for evaluating genetic diversity and population relationships [40]. Allelic richness was identical to allele frequency, implying that there was no bias based on sample size. The low levels of allelic diversity indicate that the communal indigenous goat populations are experiencing low levels of nonrandom mating when bred and not maintaining allelic diversity over the entire populations. The level of allelic diversity in this study was similar to that reported for nine Chinese cashmere goats, but lower than in the Chinese cashmere [40] and Swiss [27,41] goat populations. The high mean number of alleles per locus and expected heterozygosities indicated that the overall gene diversity was high in the population.
A key point to design/update breeding programs and conservation strategies is represented by genomic characterization of the genetic diversity of breeds. Communal indigenous goat populations from the four provinces of FS, GP, LP, and NW have not been comprehensively studied for population structure and diversity, thus genetic characterization of these populations is of paramount importance. Heterozygosity was calculated to determine genetic variation ranging from 0.42 ± 0.09 in NW and 0.39 ± 0.11 in LP. Ref. [6] reported heterozygosity values ranging from 0.387 in Swazi goats and 0.41 in Tswana goats, whilst [33] reported similar values for communal indigenous goat populations. Ref. [23] showed levels of H o ranging from 0.35 to 0.41 for 14 Italian goat breeds using the same SNP Beadchip, whilst [42] revealed low values in Nubian goat populations (0.338) and high values in (0.413) rangeland goat populations. Management practices that caused breed segregation were reported to be responsible for the low H o that characterizes the Malawi populations, whereas the high H o observed for populations in East Africa most likely mirror admixture by cross-breeding [26]. An increase in heterozygosity due to pastoralism and nomadism was frequently reported in sheep [43,44] and is consistent with specific geographical patterns in Iranian goats [45]. These similarities between the studies in heterozygosity values may be attributed to the overlap of geographical locations and the use of population types under similar production environments used. More similarities in H E and H o in some of the populations may be due to the use of populations/breeds which were previously not used during the development of the 50K SNP panel, as it was noted in other Southern African indigenous goat populations such as the Binga, Chipinge, Matopo, and Shurugwi goat populations of Zimbabwe [46]. This can suggest that these values represent the genetic diversity levels of the communal indigenous goat populations.
In the study, it was shown that populations can be exploited through appropriate breeding strategies to improve productivity, as large within-population variation (90.64%) was observed. It was revealed that the level of the among-population genetic variation was lower than that reported in commercially developed breed populations including Angora with 16.12% [47,48]. This high within-population variation is supported by the population structure analysis. The PC1 and PC2 show the highest level of genetic heterogeneity and weak sub-structuring in the FS, GP, and NW provinces, which might be triggered by continuing crossbreeding with improved breeds to improve their productivity practiced in these regions (farmer communication). Limpopo Communal goat populations and commercial breeds had the highest F ST (0.084), while GP populations had a low F ST value (range 0 to 0.038). Indigenous breeds in the Southern African regions have shown high genetic diversity, but weak population sub-structuring [10]. The highest genetic distance (F ST ) was observed to be higher than 0.25, moderate to between 0.05 and 0.25, and the lowest estimate below 0.05 when microsatellite markers were used [42]. Many authors in sheep studies for genetic distance among most of the populations recorded F ST values (0.101), (0.042) and almost negligible (<0.05) and/or moderate (0.05 < F ST < 0.25) [49,50]. Significant genetic distance estimates among populations was revealed by some of the authors. Among the goat populations, there is relatively low to moderate genetic sub-differentiation. An indication of significant differentiation among populations is considered to be a fixation index (F ST ) of about 0.15 [42].
The commercial goats from different geographical areas and breeding backgrounds provided stronger sub-structuring. According to [33], samples were taken from a farm with closed populations of Savanna and Kalahari Red in Griquastad, Northern Cape, and these populations were used as reference populations. The Kalahari Red and Savanna populations that were found to be kept in commercial farms, corporations, and research stations showed levels of admixture. In shaping the differentiation of breeds, isolation by geographical distance can play an important role, however, in the populations under study, management and breeding strategies seemed to have more effect. The PCA, ADMIXTURE, and pairwise F ST results are in agreement, thus emphasizing the exclusive genetic diversity and clustering populations according to the breed type and production system of communal indigenous goat population of Limpopo.
Searching for genomic regions that show high levels of differentiation between populations or by looking for regions of low genetic diversity within a population is possible by detection of signatures of selection from genotype data [51]. In this study, 442 and 5 SNPs based on the F ST and PC analysis, respectively, were identified. The 5 SNPs were related to PC2, which showed communal indigenous goat population of Limpopo as having a unique genetic background. The SNP with the highest F ST value (F ST 0.72) was identified to be at chromosome 8 and sits on the UHRF2 gene. Ref. [52] reported that selected genes in other studies on sheep were mostly enriched in biological processes, including metabolic processes, regulation of biological processes, developmental processes, and reproductive processes. An additional four important genes from the top ten outlier SNPs were discovered with the F ST value (F ST 0.69) at chromosome 8 sits on the GLDC gene, F ST value (F ST 0.68) at chromosome 6 sits on the NDST3 gene, and F ST value (F ST 0.68; 0.67) at chromosome 13 sits on the CFAP61 and CUBN gene, respectively. The candidate genes found in this study were involved in pathways involved in metabolism, disease response, thermoregulation, and longevity. In other studies that investigated coat color, Ref. [51] detected regions on chromosomes 5, 13, and 18 in the black vs. white comparison, regions on chromosomes 5, 9, and 13 in the white vs. black and red comparison, and regions on chromosomes 8, 22, and 29 in the red vs. white comparison. Analysis of the data by [51] for the three groups of solid coat colors detected signatures of selection near at least five genes that are known to be involved in the color and pattern definition of the coat: ADAMTS20, MC1R, ASIP, SOX18, and TIMP3. Two of these genes (MC1R and SOX18) were specifically detected in the comparison between solid black and solid white individuals, whereas the remaining three were identified in both solid black vs. solid white and solid black and red vs. solid white groups.
The 133 candidate genes were adjacent to the outlier's loci that associated significantly with pathways involved in the immune system, thermoregulation, and longevity. In the study, chromosome 8 had the highest number with 37 SNPs, and chromosomes 28 had 4 SNPs. The results show that twenty-five outliers were significant to each PC, however, only 5 could be retained, which were found on chromosomes 10, 15, 20, and 21, and were associated with the PC2, while the gene search did not reveal any genes. Ref. [53] has emphasized that studying differences in different genes could elucidate the underlying genetics of the differences between Sudanese goat breeds in regard to bone formation and body measurement characteristics, which is applicable to the South African communal indigenous goat populations. Significant roles in the processes of oocyte maturation, regulation of follicular growth, and ovulation in livestock is played by other enriched pathways, such as bile secretion (hsa04976), retrograde endocannabinoid signaling (hsa04723), the hippo signaling pathway (hsa04390), the oxytocin signaling pathway (hsa04921), and ovarian steroidogenesis (hsa04913) [50].
Genetic diversity observed between these communal indigenous goat populations with more attention on Limpopo populations provide a unique opportunity to conserve the population. The study confirms the observation that communal indigenous goat populations kept closer to the metropolitan cities have different production values, and therefore market seems to be more of a breeding factor for GP communal indigenous goat population. Crossbreeding with commercial goat populations, especially the use of the dominant Boer goat bucks, was in evidence. The impact of crossbreeding if not given attention by policy makers will have a negative impact on the genetic diversity of the communal indigenous goat populations and to their genetic merit, and thus conservation programs will be affected.

Conclusions
The communal indigenous goat populations showed a weak sub-structuring with admixed individuals with evidence of gene flow and shared genome ancestry with commercial goat populations, while the Limpopo population separated to form its own cluster. Further research is required to identify the genomic regions which are associated with different important environmental adaptation and economical traits in communal indigenous goat populations using the recent updated Illumina Goat 50K SNP BeadChip in other goat producing provinces of South Africa. The patterns of genetic diversity and population structure showed that effective management and utilization of South African communal indigenous goat populations is important, and effort should be made to maintain these unique genetic resources for conservation.  Funding: Research grants from the Department of Agriculture, Land Reform and Rural Development's Directorate Genetic Resources and Directorate Policy Research and Support funded the fieldwork data used for this paper from the project "Characterisation of indigenous goats" as well as full genome sequences of indigenous goats. The first author would like to acknowledge support by the University of Limpopo for allowing him to be their Ph.D. student.