Genetic Diversity in Four Populations of Nguni (Zulu) Sheep Assessed by Microsatellite Analysis

Abstract

Zulu sheep are found mainly in the rural KwaZulu-Natal province and the numbers are declining due to indiscriminate inbreeding. There is thus a need for phenotypic and genetic characterisation as a first phase for planning conservation strategies. Zulu sheep populations sampled were from Makhathini research station (MS) (n=33), University of Zululand (UZ) (n=21), a community at KwaMthethwa (KM) (n=32) and from Msinga (EM) (n=33). One European breed Appenninica (AP) was used as out group. Microsatellite analysis using 29 microsatellite loci was used in this study. Among the Zulu sheep, the mean number of alleles per locus was the lowest (3.86) in UZ and the highest (6.24) was realised in EM. The mean values of observed and expected heterozygosity were 0.57 and 0.61, respectively. Neighbour-joining tree showed two main Zulu sheep clusters: the UZ, KM and MS sheep populations clustered together and the second cluster included only representatives from the EM population. The STRUCTURE analysis showed that KM, AP and EM were founded in separate clusters, whereas UZ and MS clustered together. The study demonstrated that there was a common origin of the population from the research stations (MS and UZ populations). It also demonstrated that the EM had a different history for the other three populations. This work suggests that exchange of rams could be useful in reducing inbreeding when considering conservation breeding programmes.

Key Words:

• Local sheep breeds
• STR markers
• Genetic variability
• Breed assignment

Introduction

Nguni sheep are divided into three groups: the Zulu, the Swazi and the Pedi (Kruger, 2011). The arrival of Nguni sheep with the Nguni people to South Africa has been traced back to between 200 and 400 AD (ILRI, 2007; Kruger, 2011; Ramsay et al., 2000). One group of Nguni people (i.e. Zulu) is believed to have arrived on the east coast of KwaZulu-Natal and then dispersed further South depending on the suitable conditions for growing crops (Kruger, 2011). Their sheep are thus referred to as Zulu sheep. This breed is thus found over the entire KwaZulu-Natal province with larger populations in the northern part of the province at places such as Jozini, Msinga (EM), Nongoma compared to southern regions (Mavule et al., 2013). Research flocks were established in Makhathini research station (MS) near the Jozini dam and University of Zululand (UZ) close to Empangeni (Kunene et al., 2009). These research flocks are also used as the reference herds for the Zulu (Nguni) sheep.

Zulu sheep are mainly kept by the rural farmers (Kunene and Fossey, 2006). Their characteristics, such as adaptation to humid and hot conditions, ability to tolerate internal and external parasites, toleration of tick-borne diseases, ability to walk long distances and favourable foraging behaviour are the main reasons why the rural farmers keep these sheep (South African Indigenous Breeds, 2012; Nyamukanza et al., 2010; Kunene and Fossey, 2006; Ramsay et al., 2000). These local breeds do not require expensive concentrate feeds but can utilise low grade crop by-products that are high in roughage (Köhler-Rollefson et al., 2009). Such characteristics are of significance to poorly resourced farmers who cannot afford high input production and thus may be challenged by keeping high yielding sheep breeds.

At the same time, because it is a small breed (Kunene and Fossey, 2006), some of the perception by farmers is that the Zulu sheep breed is not a suitable genetic type for the current market conditions (Mavule, 2013). For this reason, they practice indiscriminate crossbreeding especially with the Dorper and Merino (to increase weight) in areas such as Hlabisa, Nquthu, Babanango, Nkandla and Ulundi (Mavule et al., 2013). This breed is thus under threat and may become extinct if interventions are not put in place. It has been demonstrated that crossbreeding resulted in replacement and dilution of locally adapted breeds with the exotic ones (Köhler-Rollefson et al., 2009). Mavule (2013) reported that the effect of inbreeding may be another factor contributing to the decline of Zulu sheep number in some areas of KwaZulu-Natal. The conservationists in the Department of Agriculture, Forestry and Fisheries (South Africa) have recognised the significance of planning a maintenance programme of the Zulu sheep (DOA, 2006).

One requirement for formulating and implementing a conservation programme is information on the genetic structures of the available populations. In the case of Zulu sheep, the lack of historical information as well as pedigree data justifies the use of molecular markers to generate genetic data. Studies done so far on genetic variation of Zulu sheep have been reported by Kunene et al. (2009) and Hlophe (2011). The authors used RAPD markers to determine the genetic variability between three Zulu sheep populations. Furthermore, work by Soma et al. (2012) using microsatellites showed genetic relationships of Nguni sheep to other South African sheep breeds. There is thus limited genetic data available on Zulu sheep population in KwaZulu-Natal.

Microsatellite DNA markers have several advantages over RAPD markers and genotypes of an animal can be determined at many known loci (Dorji et al., 2010). They have been demonstrated to be good tools for measurement of a wide range of population parameters (Farid et al., 2000). These markers have proven to be useful in assessing the genetic variability and the structure of various sheep breeds and to explain the amount of genetic relationship in closely related populations (Arora et al., 2011; Cinkulov et al., 2008; Lasagna et al., 2011; Muigai et al., 2009; Peter et al., 2007). The Food and Agriculture Organization of the United Nations (FAO, 2011) also recommends that microsatellites should be used when livestock populations and breeds are characterised for genetic variation. It is thus important that genetic data for Zulu sheep based on microsatellites is generated and evaluated. The aim of this study was to give a preliminary indication of the genetic structure of four Zulu (Nguni) populations from KwaZulu-Natal using microsatellite analysis.

Materials and methods

Animal sampling and microsatellite analysis

The four studied Nguni populations were reared in ten flocks; the studied samples were taken in eight of them. A total of 119 blood samples of Zulu sheep, of both sexes, were randomly collected from 33 individuals (one farm) at MS, 32 (four farms) from KwaMthethwa (KM), 21 (one farm) from UZ and 33 (three farms) from EM. These areas are located at 28°37’S:31°55’E, 25°27’S:32°10’E, 28°51’S:51°51’E and 28°44’S:30°27’E, respectively. The physical distance between KM and UZ is 40 km, whereas it is 250 km between MS and UZ, and Msinga is 260 km from UZ. In addition, 30 Italian Appenninica breed (AP) were included as an out group.

The geographical locations from which the four Zulu populations were sampled are shown in Figure 1. The GenElute Blood Genomic DNA kit (Sigma-Aldrich, St. Louis, MO, USA) was used to extract the genomic DNA. Twenty-nine microsatellite loci (Table 1) were selected from the list of recommended markers for genotyping analyses in sheep breeds (FAO/ISAG, 2004). The markers were selected based on degree of polymorphism and their position in the sheep genome. The distribution of the markers in the genome involved 19 chromosomes with only four (1, 5, 9 and 17) associated with more than two markers (Table 1). The microsatellite markers were optimised for multiplex PCR amplification using a Biometra TGradient 96 with the conditions reported in previous work by Lasagna et al. (2011). The multiplex PCR products were pooled to allow the analysis of more microsatellites in each electrophoresis. The size of the fragments was determined using an automated DNA sequencer (ABI PRISM 3130 xl; Applied Biosystems, Foster City, CA, USA) and GeneMapper version 4.0 software (Applied Biosystems).

Statistical analysis

Given the polymorphism information content (PIC) for each of the 29 microsatellite loci, the observed and expected heterozygosity in the four Nguni sheep breeds were computed using the MICROSATELLITE TOOLKIT (Park, 2001). The POPGENE 3.2 (Yeh et al., 1999) software package was used to estimate the number of alleles observed at each locus and the mean number of alleles per breed. A test for Hardy-Weinberg equilibrium was done using GENEPOP 4.0 software (Raymond and Rousset, 1995). As according to the algorithm illustrated by Guo and Thompson (1992), a Markov Chain Monte Carlo method (20 batches, 5000 iterations per batch and a dememorisation number of 10,000) was applied to estimate unbiased exact P values. Wright’s F-statistics was estimated to assess the population genetic structure of the four sheep breeds. Fixation indices per locus (F_IS, F_IT and F_ST) were calculated according to Weir and Cockerham (1984) using the FSTAT 2.9.3.2 software package (Goudet, 2002). The F_IS for each population was calculated via bootstrapping using 1000 replicates with GENETIX 4.05 software (Belkhir et al., 1996-2004). Using the locus by locus analysis of molecular variance (AMOVA) procedure, the significance of the fixation indices was tested using ARLEQUIN 3.11 software (Excoffier et al., 2005). Reynolds weighted genetic distance (Reynolds et al., 1983) among populations was calculated using PHYLIP software 3.6 (Felsenstein, 2005). Bootstrap values were obtained with 1000 replicates over the loci. The D_AS genetic distance was estimated using POPULATION software (Chakraborty and Jin, 1993) among individuals. Neighbour-joining methodology was applied, and a tree was graphically built from the inter-individual distances using the MEGA 4 package (Tamura et al., 2007) to highlight the differentiation and potential associations among individuals. STRUCTURE version 2.3.3 (Pritchard et al., 2000) was used to assess the most probable number of partitions in the dataset without the assumption of the breed identities. The assignment of individuals to populations considered an ancestry model with admixture and correlated allele frequencies. One hundred independent runs with 500,000 MCMC (Markov Chain Monte Carlo) iterations and a burn-in of 20,000 steps were performed for 2≤number of clusters (K)≤10 to estimate the most likely number of clusters present in the dataset. The most likely K value best describing the substructure of the populations under study was identified using the K statistic as described by Evanno et al. (2005). The clustering pattern was visualised using the software DISTRUCT 1.1 (Rosenberg, 2004).

Results and discussion

A total of 253 alleles were detected across the 29 investigated loci in the studied breeds (Table 1). The number of alleles at each locus (Table 1) ranged from 6 (SPS115) to 19 (HSC), whereas the mean number of alleles per breed ranged from 3.86 (UZ) to 6.29 (EM) (Table 2). Polymorphism information content (PIC) revealed that the majority of the markers used in this study had high values per locus with the exception of four markers (D5S2, ILSTS11, BM8125 and SPS115) with medium value (<0.5) and one marker (ETH10) with low value (0.37) (Table 1). Because of the significant deviation from the Hardy-Weinberg (H-W) equilibrium observed in three out of four studied sheep populations, the loci TGLA122, ETH225 and ETH10 were excluded from further analysis (Table 1). Mean values of observed and expected heterozygosity were 0.57 and 0.61, respectively (Table 2). These values are similar to those obtained by Soma et al. (2012) on a study involving twenty South African sheep breeds (including Nguni breed) using a panel of 12 microsatellite markers. Besides, the observed and expected heterozygosity values are reasonably high suggesting that the total analysed population is characterised by a noticeable genetic variability. KwaMthethwa had the lowest values (0.55 and 0.56, respectively); whereas EM had the highest values (HO=0.60; HE=0.65). The F_IS value per population was significantly different from zero in MS and EM (Table 2). This index ranged from a minimum value 0.0333 in UZ population to 0.1178 in MS. Relatively higher values in MS indicate a strong presence of inbreeding in this population. Similarly, the difference between the expected and the observed heterozygosity was higher in this population showing excess homozygosity. However, MS had a higher mean number of alleles (5.45) than UZ and KM (3.86 and 4.45), indicating a relatively higher genetic variation. Although the mean number of observed alleles was the highest (6.29) in EM, the inbreeding coefficient was still higher than that for UZ and KM suggesting that although there is a high level of genetic variability, individual animals in this population were inbred. The increased variation for MS and EM populations reflect the sample animals which were widely distributed within the two rural areas, while UZ and MS were research flocks in the two research stations. A small number of alleles per locus in UZ may be the result of small number of founder animals for this population. The level of polymorphism observed in this study followed the recommendation of the FAO/ISAG (1998) on the minimum number of alleles for genetic studies.

The mean F_IS value, which describes the heterozygote excess compared to the expected value in the subpopulations (the different populations in our study), was 0.0662 (P<0.05). Several markers, showed a small amount of homozygote excess: ETH10, TGLA126, D5S2, ETH225, BM6526, ILSTS11, TGLA122, CSSM66, OarCP20 and OarFCB48 revealed the highest values (Appendix Table). The other markers were characterised by values lower than 0.100.

The mean F_IT index was 0.159 (P<0.05) (Appendix Table). This suggests the presence of a heterozygote deficiency within the total analysed populations. The value of the last mean fixation index, F_ST (0.146) (P<0.05), indicated the existence of segmentation among the subpopulations (Appendix Table). A similar value has been observed in Nigerian sheep breeds as reported by Agaviezor et al. (2013). Pairwise genetic differentiation among populations is shown in Table 3. According to the F_ST genetic distance estimates, the closest populations were UZ and MS (0.09), whereas the most genetically distant populations were KM and EM (0.18). The close genetic distance between UZ and MS was expected because some of the areas where founder sheep for these two populations were purchased were the same. On the other hand, EM and KM are two different communities/rural areas with different micro climates. Forage availability in these areas is also different. Msinga is known for the severe cold spells but also the animals have adapted to fending for themselves in rocky areas with minimum forage in summer and winter. The opposite occurs in the KwaMthethwa area; in these places Nguni sheep graze on dry grass and trees in winter and the forage availability is improved in summer (Nyamukanza et al., 2010). Different adaptation and the selection history over a long time in these two groups may have resulted in the large genetic distance between KM and EM. Genetic relationships based on Reynolds’ genetic distances among the populations were visualised through a neighbour-joining tree (Figure 2), showing two main clusters: the UZ, KM and MS sheep populations clustered together and the second cluster included EM Zulu population. This shows the accuracy of the microsatellite markers in differentiating populations because the EM population proved to be genetically more different from the other three populations. KwaMthethwa is geographically closer to UZ, and UZ and MS have a common history of founder animals as the research stations.

The tree based on inter-individual D_AS distances using the neighbour-joining algorithm (Figure 3) showed four defined clusters for AP, KM, UZ and EM. Makhathini research station exhibited two main clusters, with some individuals misplaced within the tree. Some individuals belonging to the MS population were observed to cluster towards the UZ population implying some genetic similarity. This trend was also demonstrated by the smallest genetic distance between these two populations in Table 3. Although geographically UZ and MS are 250 km apart, these are two research institutions work in collaboration to conserve the Zulu sheep. Some of the founder sheep in UZ were obtained in the same areas where some of the MS founder individuals were reared.

The structure of the populations was analysed following Rosenberg et al. (2002) using a Bayesian approach that inferred the number of clusters (K) present in the population, allowing the detection of differences among populations and hidden substructure within them (Figure 4). The highest K values were obtained for K=4 (Appendix Figure).

When K=2, two clusters were defined. The out group AP formed one cluster, whereas the Zulu populations clustered together. In the case of K=3, two of the Zulu sheep populations (KM and EM) formed a separated cluster, whereas other two Zulu sheep populations (UZ and MS) could not differentiate at this level of STRUCTURE analysis. At K=4, KM, AP and EM were in three separate clusters, whereas UZ and MS clustered together. This suggests that UZ and MS are very similar. However, the existence of clear differences among KM, EM and the other two Zulu sheep populations indicate that there is genetic material in UZ and MS populations that is not found in KM and EM. These results implicate the uniqueness of Zulu sheep populations in the two communal areas (KM as well as EM) which are 300 km apart. Mavule (2013) conducted a survey which included enquiry of the location of the founder flocks. The results of the survey showed that the initial stocks in these two areas had no common origin. To maintain diversity between populations of livestock in the rural farming system, it has been suggested that farmers may exchange rams for breeding among populations that are not closely associated genetically (Kunene et. al., 2009; Mavule, 2013).

The results of Bayesian cluster analysis are summarised in Table 4, where the average q values in each clusters are shown for the different populations. The membership fraction among the breeds ranged between 0.907 in MS and 0.978 in KM. The high average of assignment of the individuals of KM and EM showed the existence of genetic differences compared to the other two Zulu sheep populations.

Results obtained from this genetic analysis could be used to reflect in great detail on the genetic diversity of Zulu sheep populations of KwaZulu-Natal. This is much more elaborate than previously reported (Kunene et al., 2009; Hlophe, 2011).

Conclusions

This work shows that microsatellite markers could successfully be used to investigate the genetic structure of the Nguni (Zulu) sheep populations, providing sufficient detail that could be useful for breeding programmes. Even if few populations were sampled for this research, the microsatellite markers provided sufficient evidence that although there is a substantial level of inbreeding, the populations of Zulu sheep are unique. This may be the basis to implement a ram exchange programme among different areas as a means to promote sustainability and conservation of this breed.

However, further conclusive results can be made after analysing more populations and also in comparison with the Dorper, Damara and Merinos which, in some areas, have been used to crossbreed the Nguni sheep populations. There is an indication, however, that a programme on ram exchange between the communities may improve the sustainability of this breed.

Acknowledgments

This research was financially supported by NRF and MAE, Executive Programme of Scientific and Technological Cooperation between the Italian Republic and the Republic of South Africa for the years 2011-2013 - Projects for the exchange of researchers; and by DAFF, Department of Agriculture Forestry and Fisheries, South Africa. A special thanks to Dr. Vincenzo Landi (University of Cordoba, Spain) for his collaboration in STRUCTURE analysis.

Figure 1. Geographical sampling areas.

Figure 2. Neighbour-joining tree obtained from Reynolds’ genetic distances among the populations of Zulu sheep. UZ, University of Zululand; MS, Makhathini research station; KM, KwaMthethwa; EM, Msinga; AP, out group Italian Appenninica. At each node, bootstrap values above 50% are shown.

Figure 3. Neighbour-joining tree obtained from inter-individual D_AS genetic distances.

Figure 4. STRUCTURE cluster analysis of the sample. UZ, University of Zululand; MS, Makhathini research station; KM, KwaMthethwa; EM, Msinga; AP, out group Italian Appenninica.

Table 1. Microsatellite loci, chromosomal position, number of alleles observed, size range, mean polymorphism information content, and number of populations deviated from the Hardy-Weinberg equilibrium per locus in the four Nguni (Zulu) sheep populations.

Locus	Chr	Na	SR, bp	Mean PIC	HWE pop
BM1824	1	7	169-177	0.52	0
BM6506	1	7	192-208	0.52	1
INRA006	1	10	101-141	0.51	1
OarFCB11	2	13	117-151	0.73	1
OarFCB20	2	15	90-120	0.60	0
OarCP34	3	8	106-130	0.54	1
McM527	5	11	166-192	0.62	0
D5S2^°	5	7	190-206	0.42	0
ETH10^°	5	12	205-211	0.37	3
RM006	5	10	113-141	0.63	1
ETH225	9	13	134-156	0.56	3
CSSM66	9	14	165-225	0.53	1
ILSTS11	9	9	201-288	0.42	1
INRA35	12	11	112-138	0.68	1
BM6526	26	11	151-175	0.63	2
CRSD247^°	14	13	213-259	0.61	0
INRA63	14	10	157-195	0.54	0
SPS115	15	6	234-250	0.47	1
MAF65	15	7	128-144	0.66	0
TGLA126	16	13	119-147	0.69	1
MAF214	16	13	183-267	0.50	0
MAF209	17	12	107-137	0.66	0
BM8125	17	8	108-140	0.42	0
OarFCB48	17	9	141-171	0.51	2
TGLA122	18	13	138-174	0.76	3
OarFCB304	19	15	151-194	0.74	0
HSC	20	19	268-300	0.82	0
BM1818	20	12	256-284	0.63	0
OarCP20	21	10	68-94	0.52	0

Chr, chromosomal position; Na, number of alleles observed at each locus; SR, size range; mean PIC, mean polymorphism information content per locus; HWE pop, Hardy-Weinberg equilibrium per locus.

°Relative to cattle linkage map (not mapped in Ovis aries).

Table 2. Sheep populations studied, sample size of each population, mean number of observed alleles, mean observed and expected heterozygosity, and inbreeding coefficient per population.

Population	Sample size	MNA±SD	HO±SD	HE±SD	FIS
UZ	21	3.86±1.41	0.56±0.021	0.58±0.032	0.0333
MS	33	5.48±2.47	0.56±0.016	0.63±0.027	0.1178^a
KM	32	4.45±1.88	0.55±0.017	0.56±0.033	0.0348
EM	33	6.29±2.16	0.60±0.017	0.65±0.023	0.0787^a
Mean value		5.02±1.98	0.57±0.02	0.61±0.03	0.0662

MNA, mean number of observed alleles; SD, standard deviation; HO, mean observed heterozygosity; HE, mean expected heterozygosity; F_IS, inbreeding coefficient; UZ, University of Zululand; MS, Makhathini research station; KM, KwaMthethwa; EM, Msinga.

^aSignificantly different from zero (P<0.05).

Table 3. Estimates of pairwise F_ST distances between the analysed breeds.

Population/breed	UZ	MS	KM	EM	AP
UZ	-	0.09	0.13	0.12	0.19
MS		-	0.12	0.10	0.16
KM			-	0.18	0.22
EM				-	0.17
AP					-

UZ, University of Zululand; MS, Makhathini research station; KM, KwaMthethwa; EM, Msinga; AP, out group Italian Appenninica.

Table 4. Proportion of membership of the five sheep populations for K=4.

Population/breed	Clusters
	1	2	3	4
UZ	0.011	0.030	0.954^°	0.005
MS	0.015	0.032	0.907^°	0.045
KM	0.005	0.978^°	0.012	0.005
EM	0.011	0.004	0.008	0.977^°
AP	0.984^°	0.003	0.007	0.006

UZ, University of Zululand; MS, Makhathini research station; KM, KwaMthethwa; EM, Msinga; AP, out group Italian Appenninica.

°Contributions higher than 0.900.

APPENDIX

F-statistics indices per locus in the whole population.

Locus	FIT	FST	FIS
BM8125	-0.014	0.033^*	-0.048
BM1818	0.066	0.023	0.044
BM1824	0.256^**	0.277^**	-0.028
CSRD247	0.220^**	0.146^**	0.086
CSSM66	0.281^**	0.177^**	0.126^**
ILSTS11	0.335^**	0.232^**	0.135^*
INRA6	0.231^**	0.206^**	0.031
MAF65	0.114^**	0.029^*	0.087^*
SPS115	0.154^*	0.194^**	-0.050
TGLA122	0.177^**	0.049^**	0.134^**
BM6506	0.149^**	0.187^**	-0.046
ETH225	0.289^**	0.091^**	0.218^**
ETH10	0.634^**	0.406^**	0.384^**
INRA63	0.140^*	0.142^**	-0.002
TGLA126	0.405^**	0.123^**	0.322^**
BM6526	0.343^**	0.210^**	0.169^**
OarCP20	0.221^**	0.119^**	0.116^*
OarCP34	0.135^*	0.162^**	-0.032
OarFCB304	0.100^**	0.112^**	-0.013
RM006	0.052	0.076^**	-0.026
D5S2	0.369^**	0.189^**	0.222^**
HSC	0.114^**	0.083^**	0.034
MAF209	0.126^**	0.076^**	0.055
McM527	0.240^**	0.165^**	0.089^*
OarFCB11	0.141^**	0.090^**	0.056
OarFCB48	0.155^**	0.061^**	0.101^*
FCB20	0.170^**	0.248^**	-0.104
MAF214	0.223^**	0.202^**	0.027
Total	0.159^*	0.146^*	0.015^*

*P<0.05;

**P<0.01.

ΔK calculated as in Evanno et al. (2005) from K=3 to K=9.