Genetic diversity, population structure and ancestral origin of KwaZulu-Natal native chicken ecotypes using microsatellite and mitochondrial DNA markers

Abstract

South African native chicken breeds are no exception to the declining of local domestic breeds that has long been recognised and found associated with the commercialisation of breeding in domestic animals. The aims of the study were: (i) to provide a comprehensive view of genetic variation in four KwaZulu-Natal indigenous chicken populations (Jozini, Newcastle, Pietermaritzburg, and Port Shepstone), and (ii) to estimate the extent of differentiation of the village populations from three conserved South African indigenous pure breeds (Potchefstroom koekoek, Ovambo, and Venda) by genotyping individuals at 19 autosomal microsatellite loci. Finally, (iii) new information to the history of KwaZulu-Natal indigenous chicken populations was made available by exploring their phylogenetic relationship and their possible maternal origin through the mitochondrial DNA. The results suggested noticeable genetic diversity within and between ecotypes with clear sub-structuring between them. The indigenous populations had high genetic diversity (observed heterozygosity ranging from 0.61 in Pietermaritzburg to 0.70 in Jozini) while conserved populations showed considerable within population inbreeding coefficient (from 0.01 in Potchefstroom koekoek to 0.18 in Ovambo). Median-joining network analyses indicated the dominance of haplogroup E suggesting a likely Southeast Asia and/or Indian subcontinent origin. The presence of haplogroup B and C not only emphasises multiple maternal origin but also highlights genetic introgression of local chickens with commercial genotypes. These results highlighted the importance of local breeds as a genetic reservoir; moreover, the conservation of local breeds may play an important role in the local economy as a source of high-quality products for consumers.

Highlights

• The genetic variation and differentiation of KwaZulu-Natal indigenous chicken populations and pure South African chicken breeds is investigated.

• The village populations showed a noticeable genetic variability with clear sub-structuring between them.

• The results can be used to improve a sustainable breeding and conservation programs to control possible genetic dilution with commercial genotypes.

Keywords:

• Genetic characterisation
• phylogeny
• molecular markers
• indigenous chickens
• South Africa

Introduction

Domestic chickens have long been used for religious activities, decorative arts, and entertainment (Liu et al. 2006). However, their history of domestication and routine of distribution over the whole world is still unstable and open to debate (Zhang et al. 2017). Asia and African continents have a long history of crops and livestock that interchanged between them due to maritime and terrestrial corridors (Boivin and Fuller 2009). Mwacharo et al. (2011) reported that East Africa and Asia experienced important maritime exchanges following the monsoon wind patterns of the Indian Ocean. Cattle, pigs, chickens, broomcorn millet and bananas are some of the livestock and crop domesticates that arrived in Africa from Asia through maritime routes (Boivin and Fuller 2009). Domestic chickens are now recognised as an important component of the rural household livelihood by providing good quality meat and eggs, used as a source of income, and as gifts to strengthen social relationships at a cheaper cost (Weimann et al. 2016).

For centuries, locally adapted chicken breeds have been exposed to unfavourable environmental conditions in Southern African Development Community regions. As a result, they have developed adaptive traits and genetic combinations that allow them to survive under local conditions including water and food availability, climate, and disease outbreaks (Köhler-Rollefson 2004). However, breeding programs in domestic animals have gradually favoured the use of highly productive exotic breeds which has consequently led to decline of population sizes of low performing native breeds (Abebe et al. 2015). There is a concrete risk of extinction of native domestic chicken breeds (FAO 2009; Fathi et al. 2017). Approximately 33% of local chicken breeds globally have been classified as endangered, thus their valuable genotypes could be lost through extinction (Mahammi et al. 2016; Strillacci et al. 2017).

Although local chickens are still widespread in the rural villages of South Africa, their genotype is under threat due to indiscriminate crossbreeding practices (Reardon et al. 1998). Halima et al. (2007) reported that random distribution of exotic chicken breeds by both governmental and non-governmental organisations, resulted in dilution of the indigenous chicken genetic stock, leading to biological diversity loss. Imported chicken breeds over the years have been used in developing countries by genetic programmes to improve productivity (Dana et al. 2011). According to Khobondo et al. (2015) exotic breeds (such as White and Brown Leghorns, Rhode Island Red, New Hampshire, Australorp, and White Sussex) have been introduced to improve local chicken performances.

Despite genetic erosion due to human practice, factors such as diseases and predators are also some of the reasons why these valuable genetic resources are still being lost. However, initiating conservation programmes can be helpful towards protecting their highly diverse gene pool. In animal breeding, decision-making is guided by the knowledge of breed peculiarities (FAO 2007a, 2007b). FAO (2012) reported that genetic characterisation describes and classifies animal breeds and species at molecular level using DNA-based techniques. On the other hand, livestock phenotypic characterisation classifies livestock populations and clarifies their external and production characteristics in a given environment and management, considering social and economic factors that affect them.

Poultry breeds characterisation methods have gradually been upgraded, from morphological traits to molecular levels (Mtileni et al. 2011a). The use of molecular markers can provide quantified criteria for assessing genetic diversity, either within or between populations (FAO 2009). Microsatellites and mitochondrial DNA (mtDNA) are markers used to investigate genetic makeup for many domestic animals (Groeneveld et al. 2010; Di Lorenzo et al. 2015; Lasagna et al. 2020). These markers complement each other to provide a more relevant reflection of population genetic makeup and its relationship with other populations (Ceccobelli et al. 2015). The study sought to investigate (i) the genetic variation within and between KwaZulu-Natal village indigenous chicken populations sampled in different agro-ecological zones and (ii) the extent of differentiation of the village populations from the conserved South African indigenous pure breeds. White Sussex is a British dual-purpose breed that is adapted to local conditions (Hrnčár et al. 2016) and due to its noticeable high frequencies on studied village farms, it was included in this study in order to estimate its effect on the genetic dilution phenomenon in South African indigenous chicken populations. Furthermore, the control region mtDNA variation of these populations was analysed. Haplotypes were classified in different haplogroups in order to clarify their phylogenetic relationship and to reconstruct their possible maternal origin.

Materials and methods

Ethics statement

All experimental procedures were reviewed and approved by the University of Zululand Research Ethics Committee, Reg No: UZRC 171110-030 PGM 2016/342.

Animal sampling and DNA extraction

A total of 199 blood samples were randomly collected from 4 different regions of KwaZulu-Natal (KZN), including Jozini (JO, n = 30), Newcastle (NE, n = 30), Pietermaritzburg (PM, n = 30), Port Shepstone (PO, n = 30) with 50% representation of males and females per region. Minimum of 3 mL of blood was collected from a wing vein of each animal in a vacutainer^® system tube, containing EDTA as an anticoagulant, then stored at −20 °C until the DNA extraction was made. Genomic DNA was extracted as reported by Selepe et al. (2018). All analysed individuals per region were selected from different farms (minimum of 5 farms per population) that are widely separated by geographical barriers to avoid repeated selection of closely related individuals per region. Selected individuals can then be considered as the representative sample of village populations within each area. Pure breeds of most common South African indigenous chicken were sampled as reference populations for comparison with village populations: Potchefstroom koekoek (KO, n = 20), Ovambo (OV, n = 20), and Venda (VE, n = 19). One exotic pure breed (common in the areas of KwaZulu-Natal), White Sussex (WH, n = 20), was also sampled to trace any crossbreeding.

The details of the studied chicken populations/breeds are reported in Supplementary Table S3.

Microsatellite molecular analysis

A total of 19 microsatellite loci (MCW0248, MCW0111, ADL0268, MCW0020, MCW0034, MCW0222, MCW0037, MCW0295, LEI0094, MCW0098, MCW0078, ADL0278, MCW0067, ADL0112, MCW0216, MCW0123, MCW0330, MCW0165, MCW0069 – Supplementary Table S1) were chosen from a list of 30 microsatellites recommended by the ISAG/FAO (2004) for chicken genetic studies. Markers were selected according to their level of polymorphism and also to their performance in a previous research (Ceccobelli et al. 2015). Multiplex PCR was carried out according to ISAG/FAO (2004) recommendations using an ABI ProFlexTM 3 × 32-well PCR system. The PCR conditions were as follows: initial denaturation step of 5 min at 95 °C, 35 cycles of 30 sec at 90 °C, 45 sec at the annealing temperature of each multiplex PCR, 30 sec at 72 °C and a final extension period of 15 min at 72 °C. Products from Multiplex PCR were pulled under electrophoresis to allow the analysis of several microsatellites. An automated DNA Sequencer (ABI 3500 Genetic Analyser, Applied Biosystems, Forest City, CA, USA) and GeneMapper version 5.0 software (Applied Biosystems, Foster City, CA, USA) were used to determine the size of each fragments.

Microsatellite statistical analysis

Allele frequencies, mean number of alleles, polymorphic information content (PIC) for each microsatellite loci, and observed and expected heterozygosity in the eight populations were estimated using the MICROSATELLITE TOOLKIT (Park 2001). To calculate average allelic richness and richness of private alleles for each population, the rarefaction method (Kalinowski 2004) implemented in HP-RARE version 1.0 software was used, adopting a sample of 13 individuals (Kalinowski 2005). MICRO-CHECKER version 2.2.3 (Van Oosterhout et al. 2004) was used to detect possible genotyping errors due to allelic dropout, stuttering and null alleles. A test for departure from Hardy–Weinberg equilibrium (HWE) was done using the Markov Chain Monte Carlo method (20 batches, 5,000 iterations per batch and a dememorisation number of 10,000) implemented in GENEPOP version 4.0 software (Raymond and Rousset 1995). Levels of significance were adjusted using Bonferroni procedure (Rice 1989). The F_IS for each population was calculated via bootstrapping using 1,000 replicates with GENETIX software version 4.05 (Belkhir et al. 2004). The pairwise F_ST indexes of Weir and Cockerham (Weir and Cockerham 1984) between populations and their associated 95% confidence intervals were estimated using GDA software (Lewis and Zaykin 1999).

The Reynolds weighted genetic distance (Reynolds et al. 1983) among the populations was calculated and a neighbour-joining tree was reconstructed using PHYLIP package version 3.6 (Felsenstein 2005); the dendrogram was depicted using the software package TreeView version 1.6.6 (Page 1996). Population structure across the entire dataset was analysed using a Bayesian approach implemented in STRUCTURE software version 2.3.4 (Pritchard et al. 2000) 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, correlated allele frequencies and defined sampling location for each individual. Ten independent runs with 500,000 MCMC (Markov Chain Monte Carlo) iterations and a burn-in of 200,000 steps were performed for 2 ≤ K ≤ 7 (K= number of clusters) to estimate the most likely number of clusters present in the dataset. The algorithm of Evanno et al. (2005) was adopted in order to evaluate the most probable value of K. STRUCTURE HARVESTER (Earl and vonHoldt 2012), a web-based programme, was used for collating results generated by the programme STRUCTURE. The clustering pattern was implemented in the CLUMPP programme and visualised using the DISTRUCT software version 1.1 (Rosenberg 2004).

To further investigate the genetic structure of each breed when adopting an approach without assumptions about HWE or linkage disequilibrium, Discriminant Analysis of Principal Component (DAPC) was carried out using the method implemented in the ADEGENET package (Jombart 2008) within the statistical package R version 3.3.2 (2020). DAPC was conducted without a posteriori group assignment by inferring the most likely number of genetic clusters (K) using the find clusters function. This procedure utilises K-means clustering to calculate a Bayesian information criterion (BIC) value for each potential value of K (the most likely K has the lowest BIC value) and delineates individual group assignments for DAPC.

Mitochondrial DNA molecular analysis

A subset of DNA samples randomly chosen for each KwaZulu-Natal native chicken ecotypes and pure breeds were sequenced as shown in Table 3. In relation to the complete mitochondrial sequence (accession number NC_040970.1; Liu et al. 2018) control region mtDNA amplification was performed from nucleotide position np 16,741 to np 1279 (1324 base pairs). DNA was amplified by PCR using the following primers: forward 5′-AGGACTACGGCTTGAAAAGC-3′ and reverse 5′-CATCTTGGCATCTTCAGTGCC-3′. DNA amplification was conducted in a 25 μl volume with 3 ul of 25 mM MgCl₂, 2 μl of 10 mM dNTPs, 10 mM of each primer and 0.2 μl of 5 U/μl of Taq^® DNA polymerase. PCR conditions were as follows: denaturation for 5 minutes at 95 °C followed by 35 cycles consisting of 95 °C for 30 s, 60 °C for 30 s, and 72 °C for 75 s. PCR was finalised with a 5 min extension at 72 °C. Purified PCR products were sent to central analytical facilities of Stellenbosch University, for Sanger sequencing with the forward primer 5′-AGGACTACGGCTTGAAAAGC-3′. Electropherograms were aligned, assembled, and compared using Sequencher™ 5.10 software (Gene Codes Inc., Ann Arbour, MI, USA). Fragments of 530 base pairs (from np 21 to np 550) resulted from the standardisation and optimisation of the ecotypes and pure breeds sequencing results. Additional nucleotide bases that exceeded np 550 were excluded from the analysis in order to make uniform our dataset. The mtDNA sequences have been deposited in GenBank with accession numbers MT684350-MT684456. All sequences were compared with the reference sequence (Liu et al. 2018) and classified in haplotypes and haplogroups, following the nomenclature previously published by Miao et al. (2013).

Mitochondrial statistical analysis

Genetic diversity indices such as haplotype diversity (Hd), nucleotide diversity (π), and average number of nucleotide differences (k) were estimated by DnaSP 5.10 software (Librado and Rozas 2009). Analysis of molecular variance (AMOVA) and pairwise F_ST calculations were performed using the Arlequin software package version 3.5 (Excoffier and Heidi 2006). Significance testing of the values was performed using 1000 coalescent simulations. Intrapopulation and interpopulation comparisons were performed based on the number of pairwise differences between sequences and calculated using an Arlequin integrated R script (http://www.rproject.org/).

In order to assess the potential phylogeographic origin of the samples, reference sequences from Miao et al. (2013) haplogroups were included (Supplementary Table S2). Network 5.0 software (http://www.fluxus-engineering.com/sharenet.htm) was used to draw maximum parsimony median-joining (MP) network plots.

Results

Microsatellite polymorphisms

A total of 161 alleles were observed with an average of 8.47 alleles per locus across the 19 microsatellite loci for the eight studied South African populations (Supplementary Table S1). All studied markers were found to be polymorphic. Level of polymorphism fluctuated between markers with number of alleles per locus ranging from 4 (MCW0222 and MCW0098) to 13 (MCW0034 and LEI0094). The expected frequencies of heterozygotes per locus ranged from 0.51 (MCW0216) to 0.75 (MCW0165) while that observed stretched from 0.43 (MCW216) to 0.77 (MCW0248).

Mean number of observed alleles (MNA), allelic richness (R), private allelic richness (PR), mean observed (H_O) and expected heterozygosity (H_E), and inbreeding coefficient (F_IS) per population are presented in Table 1.

Table 1. Genetic diversity of eight studied chicken populations/breeds.

Sample type	Population/breed	Sample Size	Breed acronym	MNA ± SD	R(PR)^a	HO ± SD	HE ± SD	dHWE^b	FIS [IC95%]
	Jozini	30	JO	5.11 ± 2.28	4.69(0.44)	0.70 ± 0.02	0.68 ± 0.03	0	−0.030 [−0.106 to −0.001]
Ecotypes (Village populations)	Newcastle	30	NE	4.84 ± 1.95	4.32(0.17)	0.67 ± 0.02	0.65 ± 0.02	0	−0.038 [−0.121 to −0.002]
	Pietermaritzburg	30	PM	4.63 ± 1.50	4.13(0.12)	0.61 ± 0.02	0.66 ± 0.02	2	0.074 [0.006 − 0.103]
	Port Shepstone	30	PO	5.32 ± 1.80	4.84(0.27)	0.68 ± 0.02	0.72 ± 0.02	4	0.067 [−0.002 − 0.091]
Mean (SD)				4.98 ± 1.48	4.50(0.25)	0.67 ± 0.02	0.68 ± 0.02	1.50	0.018 ± 0.060
	Potchefstroom koekoek	20	KO	3.05 ± 0.78	2.97(0.05)	0.50 ± 0.03	0.51 ± 0.03	1	0.013 [−0.124 − 0.079]
South African conserved pure breeds	Ovambo	20	OV	4.42 ± 1.61	4.27(0.33)	0.55 ± 0.03	0.67 ± 0.03	5	0.181 [0.068 − 0.231]
	Venda	19	VE	3.53 ± 0.96	3.47(0.54)	0.53 ± 0.03	0.59 ± 0.02	3	0.117 [−0.029 − 0.196]
Mean (SD)				3.67 ± 1.12	3.57(0.31)	0.53 ± 0.03	0.59 ± 0.03	3.00	0.104 ± 0.084
Conserved exotic pure breed	White Sussex	20	WH	3.42 ± 1.50	3.30(0.09)	0.47 ± 0.03	0.53 ± 0.05	0	0.113 [−0.001 − 0.164]

Mean number of observed alleles (MNA), allelic richness (R) and private allelic richness (PR) per population, mean observed (H_O) and expected heterozygosity (H_E), and inbreeding coefficient (F_IS) per population with confidence interval (IC).

^a 13 individuals

^b Number of loci that deviated from Hardy–Weinberg equilibrium (after Bonferroni correction, p < .0026, Rice 1989).

JO: Jozini; NE: Newcastle; PM: Pietermaritzburg; PO: Port Shepstone; KO: Potchefstroom koekoek; OV: Ovambo; VE: Venda; WH: White Sussex.

The village populations had higher mean number of alleles per locus compared to the conserved flocks. The MNA per population extended from 4.63 (PM) to 5.32 (PO) and has an across all the village population average value of 4.98 ± 1.48 alleles per locus. Mean allelic richness values (mean value of 4.50) were similar in all village populations, varying within a short range, from 4.13 in PM to 4.84 in PO, assuming a minimum sample size of 13 individuals.

The expected and observed average heterozygote frequencies within populations across loci were 0.68 (ranging from 0.65 to 0.72) and 0.67 (ranging from 0.61 to 0.70), respectively. The conserved flocks showed lower values of expected and observed heterozygote frequencies (0.59 and 0.53, respectively) than those recorded for the village populations. The MICRO-CHECKER analysis suggested that there were no indications of scoring error due to stuttering, large allele dropout or evidence for null alleles. Following Bonferroni correction, 1.50 of the 19 loci (a mean value in the total ecotype samples) were noted to have deviated from HWE. JO and NE populations were at HWE for all examined loci.

The inbreeding coefficient (F_IS) estimated in the village populations ranged from 0.074 (PM) to −0.038 (NE).

Genetic differentiation, distance and phylogeny

Table 2 demonstrates pairwise genetic differentiation (F_ST) between the populations that were studied. Values presented showed that the longest distance relationship was between VE and WH (0.371), while PO and NE (0.094) were the closest. The village populations had a relatively closer genetic distance to KO than to other breeds, JO was closely related (0.179), followed by NE (0.192) while PO and PM had the longest distances (0.231 and 0.236, respectively) to KO.

Table 2. Pairwise F_ST among studied breeds/populations (with confidence interval of 95%).

	JO	NE	PM	PO	KO	OV	VE
NE	0.120 [0.183–0.061]
PM	0.127 [0.191–0.065]	0.136 [0.202–0.063]
PO	0.104 [0.169–0.044]	0.094 [0.137–0.048]	0.150 [0.218–0.086]
KO	0.179 [0.287–0.073]	0.192 [0.284–0.094]	0.236 [0.312–0.165]	0.231 [0.304–0.159]
OV	0.202 [0.287–0.115]	0.228 [0.290–0.162]	0.249 [0.328–0.170]	0.219 [0.283–0.159]	0.274 [0.391–0.167]
VE	0.231 [0.307–0.144]	0.251 [0.323–0.169]	0.284 [0.346–0.217]	0.254 [0.318–0.185]	0.329 [0.454–0.200]	0.334 [0.423–0.239]
WH	0.244 [0.329–0.153]	0.230 [0.329–0.137]	0.293 [0.369–0.208]	0.263 [0.318–0.203]	0.274 [0.401–0.168]	0.349 [0.446–0.257]	0.371 [0.484–0.255]

JO: Jozini; NE: Newcastle; PM: Pietermaritzburg; PO: Port Shepstone; KO: Potchefstroom koekoek (KO); OV: Ovambo; VE: Venda; WH: White Sussex.

The consensus neighbour-joining phylogenetic tree constructed using Reynolds’ genetics distance (Figure 1), revealed three main clusters among the eight chicken populations/breeds. The first cluster comprised PO, NE, PM, and JO. The second cluster had OV and VE, while the last cluster comprised KO and WH. The conserved breeds (KO, OV, VE, WH) showed longer branches than the four KwaZulu-Natal populations studied.

Figure 1. Genetic relationship among the four studied populations and the four reference breeds, using Reynolds’ genetic distance according to the neighbour-joining algorithm. JO: Jozini; NE: Newcastle; PM: Pietermaritzburg; PO: Port Shepstone; KO: Potchefstroom Koekoek; OV: Ovambo; VE: Venda; WH: White Sussex.

Genetic structure and admixture analysis

Using Bayesian approach, individuals were clustered to 2 ≤ K ≤ 7 for admixture analysis and the results of STRUCTURE are presented in Figure 2. This analysis was made to assess identification of variation among populations and also to reveal hidden substructures within them. At the lowest K-value (K = 2), all ecotypes clustered together and the conserved population (KO, VE, OV) formed one cluster with the exotic breed (WH). KO and WH showed a homogenous population structure, forming one common structure for every K-value. At K = 3 and 4 JO, PM and OV, and VE separated and formed independent clusters respectively. Similarity between populations appeared to drop dramatically at 5 ≤ K ≤ 7 for every K value each population clustered separately. NE and PM were found to have few individuals from the PO population at K = 7. The highest similarity coefficient was at K = 7, and was found as the most probable clustering (Supplementary Figure S1(a)).

Figure 2. Clustering of the four studied populations and the four reference breeds with STRUCTURE analysis. JO: Jozini; NE: Newcastle; PM: Pietermaritzburg; PO: Port Shepstone; KO: Potchefstroom Koekoek; OV: Ovambo; VE: Venda; WH: White Sussex.

Although ecotype populations and pure breed flocks appeared to cluster independently at K = 2, while partial admixture is noticeable between OV and all ecotypes (Figure 2). The map presented in Figure 3 was used to highlight the distribution of membership coefficients according to population and geographic location. This was done using the Q-matrix averaged over similar run for K = 7.

Figure 3. Representation of sample size (indicated by diameter of the pies) and of relative admixture distribution (indicated by different colours) inferred by STRUCTURE at K = 7. JO: Jozini; NE: Newcastle; PM: Pietermaritzburg; PO: Port Shepstone; KO: Potchefstroom koekoek; OV: Ovambo; VE: Venda; WH: White Sussex. Source: figure taken from http://www.d-maps.com and adapted for illustrative purpose only.

The Bayesian Information Criterion (BIC) statistic generated by Discriminant Analysis of Principal Components (DAPC) showed that the optimal number of clusters in the data set is K = 8 (Figure 4, Supplementary Figure S1(b)), performing one more cluster by DAPC than those generated by STRUCTURE. In the DAPC analysis, 60 PCs of the PCA were retained as input to discriminant analysis, accounting for approximately 93% of the total genetic variability. The scatterplot of the first three components of the DA (Figure 4) showed extensive sharing of genetic variation among KwaZulu-Natal ecotypes. In particular, the plot showed genetic proximity and a particular affinity among NE, PO, PM, JO. While VE appeared clearly distinct from the other populations. Using the grouping function obtained in the discriminant analysis, a high proportion of individuals were found to be correctly assigned to their original group: JO (97%), NE (97%), PO (97%). The lowest score was observed for PM (90%).

Figure 4. Scatterplot of the first two principal components of Discriminant Analysis of Principal Component (DAPC) using populations as a posteriori clusters. The individuals are assigned to populations a posteriori, that is, after determining the number of clusters by the programme, instead of forcing them into known populations. Populations are labelled inside their 95% inertia ellipses and dots represent individuals. JO: Jozini; NE: Newcastle; PM: Pietermaritzburg; PO: Port Shepstone; KO: Potchefstroom koekoek; OV: Ovambo; VE: Venda; WH: White Sussex.

Mitochondrial DNA diversity

Nucleotide diversity, polymorphic site, number of haplotypes, diversity of haplotypes and number of sequences are presented in Table 3. A total of 107 mtDNA control region sequences spanning from nucleotide position (np) 21 to np 550, including 530 bp of the mtDNA control region, and belonging to eight South African chicken populations/breeds, were analysed. All the sequences revealed the presence of 36 polymorphic sites (S) ranging from six (JO, PO, NE) to 19 (OV) (Table 3). Nucleotide diversity (π) across all ecotype individuals was estimated at 0.00902, whereas between pure breeds π was equal to 0.01004. The haplotype diversity (Hd) was determined at 0.734 for ecotypes and 0.723 for pure breeds. The total number of haplotypes found was five with PM having the highest (four) and JO having the lowest (two) values. The ecotype diversity analysis showed the highest haplotype diversity (0.654) in PM and the lowest (0.485) in NE. Among pure breeds, the haplotype diversity was 0.723, with OV the highest (0.895) and WH the lowest value (0.473).

Table 3. Genetic diversity within the studied chicken ecotypes/pure breeds based on mitochondrial DNA variation.

	ID	N	π	nh	Hd	S
Ecotypes (Village populations)
Jozini	JO	14	0.00597	2	0.527	6
Pietermaritzburg	PM	13	0.01108	4	0.654	13
Port Shepstone	PO	10	0.00294	3	0.644	6
Newcastle	NE	12	0.00549	2	0.485	6
All ecotypes		49	0.00902	5	0.734	14
South African conserved pure breeds
Potchefstroom koekoek	KO	14	0.01190	6	0.868	16
Ovambo	OV	15	0.01154	8	0.895	19
Venda	VE	15	0.01013	4	0.657	12
Exotic conserved pure breed
White Sussex	WH	14	0.00657	3	0.473	10
All pure breeds		58	0.01004	10	0.723	22

N: number of sequences used; π: diversity of nucleotide; nh: number of haplotypes; Hd: diversity of haplotype; S: number of polymorphic sites.

Mitochondrial DNA phylogeny analysis

To investigate the relationships of haplotypes, we constructed a network plot for South African village chickens with the pure chicken breeds reared in the same geographic area (Figure 5). The median-joining network analysis was based on 11 haplotypes of this study along with 55 reference sequences corresponding to the 13 haplogroups (A–Z) defined by Miao et al. (2013). The median-joining network revealed a clear pattern among village chickens and pure breeds. The median-joining tree exhibited a star-like distribution profile centred on the predominant haplogroup (E), indicating a recent exponential population growth. In clade B, D and E, the core haplotypes occurred in ecotype chickens, whereas, A and C haplogroups are a peculiarity of pure breeds.

Figure 5. Median-Joining network. Black circles refer to reference haplotypes. Circles are proportional to the numbers of samples. JO: Jozini; NE: Newcastle; PM: Pietermaritzburg; PO: Port Shepstone; KO: Potchefstroom koekoek; OV: Ovambo; VE: Venda; WH: White Sussex.

According to the analysis of pairwise genetic distances between populations, PM was the most genetically distant and it was also the most heterogeneous breed even when comparing to village populations and to the conserved South African indigenous pure breeds. Furthermore, JO and NE resulted being the most closely related (Supplementary Figure S2). The maternal genetic sub-structure was also assessed within and between studied populations by analysis of sequence variations (AMOVA) following imposition of two hierarchical groupings (Table 4). A huge within-population variation contributed 76% to the total maternal variance, while the 5.24% of the observed variation was accounted to differences between the South African ecotypes belonging to this study and pure breeds reared in the same geographic areas, with an overall Φ_ST value of 0.240 (p ≤ .001).

Table 4. Hierarchical analysis of molecular variance (ecotypes vs pure breeds).

Source of variation	df	Sum of square	Variance component	Variance (%)	Fixation index^a	p Value
Among groups (ecotypes vs pure breeds)	1	18.185	0.15462	5.24	ΦCT = 0.052	.17791
Among populations within groups	6	58.602	0.55335	18.76	ΦSC = 0.198	.000***
Within populations	101	226.387	2.24146	76.00	ΦST = 0.240	.000***

^aΦ_CT = variation among groups divided by total variation, Φ_SC = variation among sub-groups divided by the sum of variation among sub-groups within groups and variation within sub-groups, Φ_ST = the sum of variation groups divided by total variation.

***p ≤ .001.

Discussion

Comprehensive knowledge of population structure and the distribution of genetic diversity in livestock are essential aspects for improving selection designs and breeds, supporting the safeguard of biodiversity, enhancing the efficient use of the breeds, and implementing conservation programs adapted to local conditions. In this study, the level of diversity of KwaZulu-Natal village indigenous chicken populations from different South African agro-ecological zones, their genetic distinctiveness from conserved South African indigenous pure breeds and the phylogenetic relationship and maternal origin were reported.

Microsatellites

All 19 microsatellites used in this study displayed a minimum of four alleles per locus with mean number of alleles (8.47) above 4, making them effective for genetic diversity judgement (Barker 1994). The mean values of H_O and H_E were 0.59 and 0.63 (above 50), respectively (Supplementary Table S1), highlighting noticeable genetic diversity across eight studied chicken populations. Ecotypes had a higher mean number of alleles and allelic richness than conserved flocks. Allelic richness is a measure of genetic diversity and an indicative of a population’s long-term adaptability and persistence potential (Özdemir and Cassandro 2018). Results showed that the village populations are more genetically diverse than the conserved pure breeds and current values are in range with those reported by Van Marle-Köster et al. (2008) and Mtileni et al. (2011a) on South African populations. Relatively higher genetic diversity has been reported previously on local chickens (Wimmers et al. 2000; De Marchi et al. 2006; Muchadeyi et al. 2007; Mtileni et al. 2011a; Özdemir and Cassandro 2018). Muchadeyi et al. (2007) and Mtileni et al. (2011a) suggested that continuous gene flow due to immigration and emigration between subpopulations could be the reason for high genetic diversity possessed by scavenging populations sharing the same geographic zone. The observed heterozygosity values in all the ecotypes were above 60%, indicating that all analysed ecotypes are characterised by noticeable genetic diversity. Additionally, this high observed heterozygosity can further be attributed to the large number of alleles detected. Conserved South African indigenous chicken breeds observed heterozygosity range (50–55%) suggest that these populations are well managed to protect their genetic variability in conservation flocks (Table 1).

Conserved South African breeds had greater genetic distance among their populations, as evidenced by high F_ST values between them (Table 2). This was expected because they are managed under an intensive production system and are subjected to closed populations with minimum immigration and emigration. In addition, these populations had restricted gene flow accompanied by non-random mating with a limited gene pool which Muchadeyi et al. (2007) reported as militating against genetic diversity. Sumreddee et al. (2019) reported that inbreeding is very sensitive to high allele frequency and actually reduces it. Low allelic frequency within conserved populations (OV, VE) may be the result of excessive inbreeding within these populations, as confirmed by high inbreeding coefficients (Table 1). Similar to findings described by Van Marle-köster et al. (2008), conserved flocks showed a high inbreeding coefficient compared to village flocks.

Animals respond to the environment and adapt phenotypically and genetically for their survival (Benatchez 2016). South African ecotypes appear to have moderately high F_ST values between them, suggesting clear sub-structuring. These results are in contrast to Muchadeyi et al. (2007) findings with Zimbabwean single structured ecotypes. The authors reported that Zimbabwean local chicken populations are not sub-structured across their agro-ecological zones and not locally adapted and restricted to their local zones. However, the study by Lyimo et al. (2014) on Ethiopian ecotypes’ genetic diversity reported that Ethiopian ecotypes appeared to cluster in three distinct groups related to their geographic location, thus complementing findings of the current study. Results from this study indicate that the analysed South African village population’s genotypes are highly diverse with well sub-structured ecotypes.

STRUCTURE-based analysis further supports the high differentiation among ecotypes of the South African chickens. The most probable clustering was found at K = 7 where each flock separated independently, suggesting clear sub-division between all the studied South African populations. These results suggest that all populations from different regions of KwaZulu-Natal showed different allelic frequencies, which is a reflection of the larger gene pool of KZN indigenous chicken populations.

Mitochondrial DNA

Vilakazi (2018) reported high phenotypic variation among South African indigenous chicken populations on traits such as plumage colour, skin colour, shank colour, comb type, ear lobe colour, and eye colour. Considerable ecotypes genetic diversity presented by the current study substantiate findings by Mtileni et al. (2011a) and Mtileni et al. (2011b) using microsatellites and mtDNA genetic markers respectively. However, our mtDNA sequence analysis revealed higher genetic diversity from conserved than village population thus contradicts with our previous findings using microsatellite markers. According to Fratini et al. (2011), microsatellites and mitochondrial DNA genetic markers have different molecular and evolutional properties thus may disclose different distribution patterns of recorded genetic diversity. The haplotype diversity in the conserved pure South African breed flocks and village flock ecotypes occur at equal frequencies, and this reveals the positive effect of controlled genetic flow in a population and further highlights the interest and the need for conservation. The lowest haplotype diversity in the conserved exotic pure breed (WH) is in agreement with Muchadeyi et al. (2008) findings using mtDNA, in which Zimbabwean ecotypes possessed greater genetic diversity compared with exotic commercial white egg layer. Fluctuation of genetic variability presented by the studied indigenous chicken ecotypes of KwaZulu-Natal agrees with the suggestions of Delany (2003) and Muchadeyi et al. (2007) that population isolation through geographic location can lead to population sub-structuring and gene pool status variability.

Four of the five haplogroups observed in the current study (A, B, D and E), were previously reported by Mtileni et al. (2011b). Although this work reveals for the first time haplogroup C in South African populations, it was previously described in other African countries such as Kenya, Madagascar (Herrera et al. 2017), Nigeria (Lasagna et al. 2020) and Ethiopia (Al-Jumaili et al. 2020). It is known that haplogroup E is one of the haplogroups showing ubiquity on chicken mtDNA profile globally. According to Mwacharo et al. (2011) observations on domestic chickens in East Africa, haplogroup E was only observed in Ethiopia and Sudan. However, the current study discovered equal representation of all studied populations in haplogroup E suggesting it as a dominant maternal lineage. This is in agreement with Mtileni et al. (2011b) previous findings which also revealed the dominance of haplogroup E in South African chicken populations. Although there is still uncertainty with the origin of this haplogroup, according to the suggestions by Oka et al. (2007) haplogroup E originated in Southeast Asia. Furthermore, the authors stated that this haplogroup was first introduced into the Indian subcontinent before reaching to other regions. Recent study stated that haplogroup E remains by far the commonest in the African continent, likely originated from the Indian subcontinent following a terrestrial route along the Nile River Basin and/or through the coastal areas of the Horn of Africa (Al-Jumaili et al. 2020). Liu et al. (2006) proposed that haplogroup B and C originated from Yunan province and/or surrounding areas and their biological history seems to be closely related (Mwacharo et al. 2011). Muchadeyi et al. (2008) highlighted the relationship of both haplogroup B and C with commercial egg layers (brown and white) and commercial broilers including chicken from the Northwest Europe respectively. The presence of these haplotypes accompanied with its lack of diversity could be revealing the recent genetic introgression of local populations with commercial genotypes. The predominance of haplogroup D has been reported in Africa and in South African regions in particular (Muchadeyi et al. 2008; Razafindraibe et al. 2008; Mtileni et al. 2011b). Liu et al. (2006) reported that this haplogroup was common in jungle fowl and gamecocks from Indonesia, India and Japan.

Conclusions

Our results highlighted noticeable genetic variability within and between KwaZulu-Natal local chickens. Microsatellites revealed clear substructuring between studied local populations, substantiating their local adaptation to respective agro-ecological zones. Ecotypes genetic introgression by the studied exotic breed is still minimum and unnoticeable thus should be further eluded. Mitochondrial DNA sequences highlighted that South African chickens can be assigned to five haplogroups with multiple potential maternal linages. Moreover, a possible recent genetic introgression of local chickens with commercial genotypes is highlighted by mtDNA; thus, there is a need for sustainable breeding and conservation programs to control possible genetic dilution. These results can improve knowledge on structure of indigenous chicken populations. This information will contribute to breeders’ organisations constitution and can be included in the FAO DAD-IS database. Moreover, these findings can enhance the local breeds productivity in order to obtain high quality products for the consumers.

Author contributions

Conceptualisation. S.C., H.L., E.L., and N.K.; Methodology. N.N., S.C., and I.C.; Software S.C. and I.C.; Validation S.C., H.L., E.L., and N.K.; Formal analysis S.C. and I.C.; Investigation. N.N., and S.C.; Resources: H.L., E.L., and N.K.; Data curation N.N.; Writing—original draft preparation N.N.; Writing—review and editing. S.C., H.L., E.L., and N.K.; Supervision E.L., and N.K; Project administration E.L., and N.K; Funding acquisition E.L., and N.K.

Acknowledgments

The authors thank the three anonymous referees for their valuable comments to the manuscript and their constructive suggestions.

Disclosure statement

The authors have declared that no competing interests exist.

Additional information

Funding

This research was financially supported by Italian Ministry of Foreign Affairs (to EL) and National Research Foundation (to NWK) – ‘Executive programme of scientific and technological co-operation between the Italian republic and the republic of South Africa for the years 2014–2017 – Projects for the exchange of researchers’. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.