Phenotypic and genetic characterization of indigenous Tswana goats

Tswana goats that were kept in communal systems in three agro-ecological regions in Botswana were characterized according to phenotypic measurements and genotypic data. Objective measurements for 123 goats included bodyweight (BW), body length (BL), heart girth (HG), height at withers (HW), and tail length (TL), while qualitative traits included coat colour and presence or absence of horns and beards. Age was estimated based on dentition. Hair samples were collected from 48 of the phenotyped animals in the largest region (central region) and genotyped with the Illumina Goat50K SNP chip. Mixed coat colour was predominant and across regions 95% of the goats were horned and bearded. Goats in the northwest region had the lowest BW and significantly higher HG values in all age groups compared with other regions. Goats over four years old in the central and northwest regions were significantly longer in body compared with the ones from the southern region. The average expected heterozygosity and inbreeding coefficient were 0.423 ± 0.03 and 0.009 ± 0.05, respectively. Principal component analysis clustered most animals, with a few outliers. The effective population size has decreased over time and at 13 generations ago was estimated at 266. There were high genetic and phenotypic variations in the indigenous Tswana goats, which should be exploited to increase performance through within-breed selection and structured crossbreeding.


Introduction
Characterization of indigenous livestock species is the key to development of proper strategies for long-term maintenance and use of genetic variation, and for guidance in decisions about future utilization and conservation strategies (Msanga et al., 2012).In southern Africa, studies have been conducted to characterize indigenous goats including Mozambican goats (Garrine, 2007), Tanzanian goats (Madubi et al., 2000;Ngulumaet al., 2018), South African goats (Pieters et al., 2009;Mdladla et al., 2016) and Namibian goats (Els et al., 2004).However, the majority of studies have been limited to microsatellite markers and genomic information on local goat characterization in the region remains scarce.The majority of indigenous populations are not managed through herd book registrations (Van Marle-Köster et al., 2015).Although crossbreeding is important for hybrid vigour and breed complementarity, if it is indiscriminate, it poses a threat to most indigenous goat populations as it results in a decline in the adapted traits, ultimately reducing the efforts at improving food security in rural areas (Rischkowsky, et al., 2007).
In Botswana, there are approximately 1.6 million goats.The main indigenous breed is the Tswana goat, which constitutes 71% of the national goat population (Botswana Agricultural Statistics, 2013).This breed is found in various geographical regions of the country under low-input management systems and contributes significantly to the livelihoods of resource-poor farmers as a source of protein and income (Monau et al., 2017).The Tswana goat is known for its unique adaptive traits such as heat and drought tolerance, and lower disease susceptibility, which make it ideal for production under stressful tropical environments (Katangole et al., 1996;Nsoso et al., 2004).However, it has received little attention in research projects.The last phenotypic characterization studies on Tswana goats were undertaken more than a decade ago (Nsoso et al., 2004), and may not reflect the current situation, because of changes in population and in production systems (Solkner et al., 1998).The only previous effort to characterize indigenous Tswana goats at genetic level was done at a research station, using 12 microsatellite markers on a very small population that was under selection (Maletsanake et al., 2013).This population did not represent the variation that occurs in the country, and failed to show inherent genetic variation of the Tswana goat population.
Although the advantages of microsatellites have been well documented (Singh et al., 2014), genomewide single nucleotide polymorphism (SNP) markers provide new possibilities for genetic characterization and biodiversity studies (Blasco & Toro, 2014).Several SNP assays have been used in the analyses of population diversity and structure of many livestock species (Lin et al., 2010;Kanyile et al., 2015;Makina et al., 2015;Lashmar et al., 2016;Mdladla et al., 2016).The commercial Goat50K SNP panel, which was developed in 2012 (Tosser-Klopp et al., 2014), offers an opportunity for genetic characterization of goats.
The purposes of this study were to phenotypically characterize Tswana goats in three agro-ecological regions of Botswana and to assess the genetic diversity and population structure of this genetic resource in the central region using the Goat50K SNP panel.

Materials and Method
The use of animals complied with the guidelines approved by the Ethics Committee of Faculty of Natural and Agricultural Sciences at the University of Pretoria (ECO42-15).
The study was carried out in Botswana from December 2015 to March 2016.The climate of the country is mainly semi-arid.Rainfall occurs during the summer months of October to April and is low, unreliable, unevenly distributed, and highly variable from year to year (Makhabu et al., 2002).The highest average temperatures generally occur from October to April and the lowest from May to August.Drought is a recurrent phenomenon, and all rivers in the country are seasonal, except for the Okavango and Linyanti/Chobe rivers in the northwest (Mogotsi et al., 2013).Geographically, soil varies from one region to another, with about two thirds of the country being covered in infertile sandy soils (red and grey desert soils) (Mogotsi et al., 2011).Four agro-ecological regions are distinguished in Botswana, with the southern region being classified as hard veldt, the central region as hard veldt with the dominance of woodland Collospermun mophane trees, the northwest as sand veldt with thick forest, lush green plains and semi-arid shrub savanna trees, and the Ghanzi as sand veldt with shrub savanna trees.
Random sampling was performed to select representative districts, villages and farms, based on knowledge of livestock (particularly goats) from livestock extension officers in each region.In each village, four to five farms were randomly selected and one to five unrelated Tswana goats per farm were measured.To avoid sampling related individuals, farmers stated the origin and familial relationships of individual animals.
Animals were grouped into five age categories based on dentition, namely no pairs of permanent incisors (<14 months), one pair of permanent incisors (15-23 months), two pairs of permanent incisors (24-35 months), three pairs of incisors (36-48 months), and four pairs of incisors (over 48 months), according to Pace & Wakeman (2003).Bodyweight (BW) was measured using a hanging scale.Body length, HG, HW, and TL were measured with a tailor's tape following the standard procedure reported by FAO (2012).Morphological descriptions such as coat colour and presence or absence of horns and beard were also recorded.One hundred and twenty-three goats were measured in total, of which 47, 54 and 22 were located in the southern, central and northwest regions, respectively.No Tswana goats were observed in the Ghanzi region.All goats were kept under extensive systems in communal areas and were not selected for production traits.
Hair samples were collected from forty-eight of the animals in the central region, as this is the largest region (142,302 km 2 ).Hair samples were collected by plucking 50 to 100 hairs from each animal, ensuring intact follicles.The hair samples were kept in labelled envelopes and transported to the University of Pretoria laboratory and then shipped to Animal Genetic Laboratory in France, Labogena DNA platform (Domaine de Vilvert, CS 80009, 78353 Jouy en Josas cedex), for DNA extraction and genotyping.
Procedure frequency was used to analyse qualitative traits and general linear model was used to analyse quantitative traits in Statistical Analysis System (SAS) (SAS Institute, 2009).Owing to low numbers of animals at 0-14 months and 15-23 months, only three age categories were used, that is, 24-35 months, 36-48 months, and above 48 months.Analysis was also confined to females owing to the small number of males (castrates and bucks) that were kept.Fixed effects were considered to significantly affect performance when P <0.05.This statistical model was implemented to analyse BW and body measurements: where Y ijk = observation on body measurement µ = population mean age i = age ( j = 1, 2, 3, age at measurements) region j = region ( k =1, 2, 3) e ijk = normal distributed random error Forty-eight animals from the central region were genotyped using the Illumina Goat50K SNP Bead chip, which contains 53,347 SNPs.Only autosomal SNPs were considered, while SNPs on the sex chromosomes and SNPs with unmapped locations to the latest reference assembly of the goat genome were excluded from the analysis.Plink version 1.07 software (Purcell et al., 2007) was used for analysis.Quality control was performed with these standard thresholds: individual call rate lower than 97%, SNP call rate less than 97%, SNPs with MAF below 0.05, and SNPs that deviated significantly from Hardy-Weinberg equilibrium (P <0.001).After quality control, genetic diversity parameters (observed and expected heterozygosity) and average individual inbreeding coefficient (F IS ) across autosomes were calculated using Plink software (Purcell et al., 2007).A Principal component analysis (PCA) was performed using genomewide complex trait analysis (GCTA version 1.24; Yang et al., 2011).ADMIXTURE version 1.23 (Alexander et al., 2009) was used to investigate population differentiation based on SNP genotype data.A crossvalidation (CV) procedure was performed on 48 animals to choose for optimal K-value with the lowest CV error values.
After quality control, pair-wise linkage disequilibrium (LD) was assessed through the correlation coefficient (r 2 ) using Plink software (Purcell et al., 2007).The command '-r2 -ld-window-kb 2000 -ldwindow-r2 0' was used to calculate the association among SNP pairs up to a distance of 2000 kb.To inspect the LD decay with physical distance, SNP pairs were sorted into bins based on pair-wise marker distance and the average of each bin was calculated.These bins were defined, namely 0-10, 10-20, 20-40, 40-60, 60-80, 80-100, 100-200, 200-500, 500-1000 and 1000-2000 kb.Effective population size (Ne) was determined based on r 2 values at various distances and assuming a model without mutation, as described by Corbin et al. (2010) using SNep version 1.1 (Barbato et al., 2015).Minimum and maximum inter-SNP distances of 0 and 1000 Mb, respectively, were used with 30 distance bins of 50 kb each.Ne estimates were subsequently calculated from r 2 values obtained for the average distance of each bin.

Results
Adult goats dominated the flock structure across regions with 29% in the 24-35 month age group, 42% in 36-48 months and 23% over four years old.The young and grower age groups of 0-14 months and 15-23 months accounted for only 2% and 4% of the flocks, respectively.However, the numbers of goats older than 48 months in the southern and northwest regions were significantly (P <0.05) less compared with central region.Across all regions, Tswana female goats were predominant (92%) (Table 1).Almost all goats (95%) were horned and bearded.Seven coat colours were observed across all regions, with mixed colour (28%), brown and white (23%) and black and white (23%) patterns being the most prominent.Evenly coloured goats, that is, brown, black or white, were less common (Table 2).Morphological differences were observed between regions, with HW being significantly (P <0.05) lower at 24-35 months in southern region compared with the northwest.Goats in the northwest region had significantly (P <0.05) lower BW and higher HG values at 36-48 months and at >48 months.Body length measurement was significantly lower in the southern region compared with central region at over 48 months.The tails of goats from the central region were significantly longer than those from the southern region at over 48 months (Table 3).Most morphological differentiation based on region could be observed in the mature goats, older than four years.The average genotyping call rate across 48 animals was 99.6%.No individuals were removed.A total of 5203 SNP markers was removed during quality control, namely 2338 because of low call rate (<0.97), 2728 owing to low MAF (<0.05) and 137 markers that violated HWE (P <0.001).This resulted in 44741 SNP markers that were retained for further analysis.The MAF of SNPs followed a uniform distribution across the chromosomes and averaged at 0.32 ± 0.13.The average observed and expected heterozygosity was 0.419 ± 0.02 and 0.423 ± 0.03, respectively.The average individual inbreeding coefficient (F IS ) was 0.009 ± 0.05.
Breed composition was assessed for the 48 genotypes using a PCA in which most of the animals were clustered, with a few outliers (Figure 1).To investigate whether there was genetic differentiation, an ADMIXTURE version 1.23 (Alexander et al., 2009) was performed and the lowest cross-validation error was at K = 1.The average LD values at given distance intervals for the population are displayed in Figure 2. The highest average r 2 of 0.44 was observed at 0-10 kb.LD declined with increasing distance between SNP pairs, and the most rapid decline was seen over the first 10 kb.The average r 2 for all pair-wise adjacent SNP autosomes was 0.067 ± 0.10, with an average distance between SNP pairs of 255.8 kb.

Distance Intervals (kb)
A graphic representation of the effective population size at each point from 900 to 13 generation ago is given in Figure 3.The results show a progressive decrease in Ne over time with an estimation of 266 animals at 13 generations ago.

Discussion
Indigenous goats are valuable genetic resources, particularly for resource-poor smallholder farmers, owing to their ability to thrive in diverse geographical environments with limited inputs (Huson et al., 2014).The flock structure of the indigenous population indicated that the flocks generally consist of mostly female animals.Farmers prefer to sell males for income, slaughter them for home consumption, or give them as gifts (Monau et al., 2017).Mdladla et al. (2017) reported a similar observation in various provinces of South Africa.Additionally, farmers prefer to keep exotic bucks such as Boer goats for crossbreeding, which contributed to lower numbers of young (i.e.0-14 and 15-23 months) indigenous goats being kept.This is in contrast with a previous study, which reported a significant number of indigenous Tswana males of 0-12 months in the flocks (Nsoso et al., 2004).A decrease in the number of indigenous breeding bucks is one of the factors that threaten the existence of an indigenous breed (Mandal et al., 2014).Tswana goats could be endangered or vulnerable, and it is important to strategize and conserve indigenous Tswana goats for sustainable utilization and future breeding programmes.
Peculiar morphological characteristics of Tswana goats include a mixture of coat colour patterns and the presence of horns and beards (Katangole et al., 1996).Farmers could unintentionally be selecting towards these traits as an environmental adaptation mechanism for thermoregulation and predation.Predation is one of the major constraints affecting goat production in Botswana (Monau et al., 2017), and animals use their horns for protection, while mixed coat colour patterns act as camouflage (Hagan et al., 2012).Coat colour is also useful in protecting deep tissue against excess exposure to solar shortwave radiation in tropical zones (Castanheira et al., 2010).The predominant mixed coat colour could be an advantage during the common seasonal temperature fluctuations in the country.The presence of beards is perceived to be associated with superior reproductive traits such as high conception, high prolificacy and high fertility rates (Gatew et al., 2015).Large-scale studies are needed to establish the true effect of these qualitative traits on adaptation, performance and the overall productivity of Tswana goats.
Bodyweight and body sizes of Tswana goat have decreased compared with the previous report (Nsoso et al., 2004).Nsoso and colleagues (2004) reported an average BW of 41.7 ± 0.5 kg and heart girth of 80.5 ± 0.3 in mature (>36 months) female Tswana goats.The difference could be due to a lack of availability of feed resources inequality and quantity.The current survey was undertaken during drought and there were substantial shortages of feed and water that could have major impacts on physical growth of animals.However, the indigenous mature Tswana goats had superior performance in terms of BW and body measurements compared with other indigenous goat breeds in South Africa (BW: 28.96 ± 0.37 kg, BL: 60.85 ± 0.29 cm) (Selolo et al., 2015) and Tanzania (BW = 28.97 ± 0.52 kg, BL = 51.60 ± 0.40 cm) (Nguluma et al., 2016).More studies on the effect of vegetation on performance should be performed for comparison in and across agro-ecological regions.
The availability of SNP arrays provides an opportunity to investigate current genetic structure and diversity of livestock for effective selection and conservation strategies (Groeneveld et al., 2010).This study is the first attempt to genetically characterize the indigenous Tswana goat, using genome-wide SNP markers.The overall informative SNPs (92.5%) observed in the study were higher than those reported by Lashmar et al. (2015) on Angora goats (82%) and Mdlala et al. (2016) on various South African indigenous goat breeds (87.1%), but comparable with Lashmar et al. (2016) on dairy goats (92.3%).The variation could be attributed to the threshold applied during QC.For instance, a SNP call rate of 98% was used for Angora goats (Lashmar et al., 2015) whilst in this study a SNP call rate of 97% was used.The results indicate, however, that the SNP chip could be used successfully for genomic studies on indigenous goats.
The results revealed high levels of genetic diversity (0.423) which is comparable with estimates by Kim et al. (2016) on the indigenous Barki goats (0.40) of Egypt and by Mdladla et al. (2016) on South African indigenous goats (0.41).Indigenous goats kept under communal systems are exposed to natural selection, where animals become genetically adapted for survival in their natural environments, while maintaining high within and between population genetic variability (Kim et al., 2016).Previous studies noted that the goat ancestor (bezoar) had broad genetic diversity, which has been conveyed to the modern goat, which has not undergone intensive selection, as has been experienced in species such as cattle (Gerbault et al., 2012).
Individual inbreeding coefficient estimates in this study were lower than the estimates reported (0.12 ± 0.16) by Maletsanake et al. (2013) using microsatellite markers on a Tswana goat population kept at an experimental farm under a controlled breeding programme.The difference could be because the samples in this study were drawn from a large population of communal areas, and care was taken to select unrelated individuals.
Principal component analysis revealed that most animals in the central region were clustered, with a few outliers.The reason for few outliers could be lack of selection and uncontrolled breeding practised in communal management systems.The low value of cross validation (K = 1) in this study is similar to the report of Kominakis et al. (2017) in Frizarta sheep.This indicates that there is no genetic differentiation or inferred clusters in the population.However, a more extensive sampling of Tswana goats covering various geographical regions of Botswana is necessary to assess the demographic history and relationship of individuals.
The average LD (0.067 ± 0.10) reported in this study is slightly lower than the report of Mdladla et al. (2015), which reported values of 0.09 ± 0.12 on Zulu, Venda and Xhosa and 0.11 ± 0.14 on South African Tswana goats.LD values, however, are influenced by factors such as history and structure of the population, sample size and strictness of SNP filtering (i.e.threshold of MAF and HWE) (Bohmanova et al., 2011).Khatkar et al. (2008) pointed out that studies with relatively small sample sizes are subject to bias and loss of accuracy.This bias may vary with inter-marker distance.Therefore, it would be interesting to confirm the LD results in this investigation with a larger number of genotyped animals.
The pattern of LD decay with distance in this population is consistent with the report of Brito et al. (2015) on various breeds of goats.Compared with other indigenous species, linkage disequilibrium declined more slowly in the indigenous Tswana goat population than in the indigenous cattle population studied by Makina et al. (2014).The low level of long range LD indicates lack of selection on the population or large effective population size in the recent past (Brito et al., 2015).The persistence of LD over short distances, however, indicates that more markers will be needed to obtain the same power to detect association (Meadows et al., 2008).LD decay determines the power of QTL detection in association mapping studies and helps to determine the number of markers required for successful association mapping and genomic selection (Mastrangelo et al., 2014).Meuwissen et al. (2001) and Qanbari et al. (2010) proposed an r 2 value of ≥0.2 as being useful for association studies.Further advancement on dense SNP panel for goats would be essential for high power association mapping and genomic selection efficiency in future breeding programmes of Tswana goats.
The observed values of Ne in this study were consistent with the report of Mdladla et al. (2015) for South African ecotypes.Large effective population size is mostly correlated with natural selection and multiple breeding objectives (Meuwissen & Woolliams, 1994), which are the basic characteristics of smallholder communal goat production (Monau et al., 2017).It is important to manage this effective population size to preserve genetic variability in the indigenous Tswana goat population for optimal utilization and sustainable development programmes.

Figure 1
Figure 1 Genetic relationship among 48 Tswana goats sampled from the central region of Botswana (PCA1 vs. PCA2) PCA: Principal component analysis

Figure 2
Figure 2 Average r 2 values at given distances (kb) for Tswana goat population

Figure 3
Figure 3 Average effective population size for central region Tswana goats, from 900 to 13 generations ago.Ne: Effective population size

Table 1
Flock structure of indigenous Tswana goats in surveyed regions of Botswana abc Different superscripts within a row differ significantly (P <0.05)

Table 2
Qualitative traits of indigenous Tswana goats in various agro-ecological regions of Botswana

Table 3
Mean (± SE) for bodyweight and body measurements of female Tswana goats ab Means with different superscripts within a column, age group and trait differ significantly (P <0.05)