- Open Access
Investigation of the genetic diversity of domestic Capra hircus breeds reared within an early goat domestication area in Iran
Genetics Selection Evolutionvolume 46, Article number: 27 (2014)
Iran is an area of particular interest for investigating goat diversity. Archaeological remains indicate early goat domestication (about 10 000 years ago) in the Iranian Zagros Mountains as well as in the high Euphrates valley and southeastern Anatolia. In addition, mitochondrial DNA data of domestic goats and wild ancestors (C. aegagrus or bezoar) suggest a pre-domestication management of wild populations in southern Zagros and central Iranian Plateau. In this study genetic diversity was assessed in seven Iranian native goat breeds, namely Markhoz, Najdi, Taleshi, Khalkhali, Naini, native Abadeh and Turki-Ghashghaei. A total of 317 animals were characterized using 14 microsatellite loci. Two Pakistani goat populations, Pahari and Teddy, were genotyped for comparison.
Iranian goats possess a remarkable genetic diversity (average expected heterozygosity of 0.671 across loci, 10.7 alleles per locus) mainly accounted for by the within-breed component (GST = 5.9%). Positive and highly significant FIS values in the Naini, Turki-Ghashghaei, Abadeh and Markhoz breeds indicate some level of inbreeding in these populations. Multivariate analyses cluster Iranian goats into northern, central and western groups, with the western breeds relatively distinct from the others. Pakistani breeds show some relationship with Iranian populations, even if their position is not consistent across analyses. Gene flow was higher within regions (west, north, central) compared to between regions but particularly low between the western and the other two regions, probably due to the isolating topography of the Zagros mountain range. The Turki-Ghashghaei, Najdi and Abadeh breeds are reared in geographic areas where mtDNA provided evidence of early domestication. These breeds are highly variable, located on basal short branches in the neighbor-joining tree, close to the origin of the principal component analysis plot and, although highly admixed, they are quite distinct from those reared on the western side of the Zagros mountain range.
These observations call for further investigation of the nuclear DNA diversity of these breeds within a much wider geographic context to confirm or re-discuss the current hypothesis (based on maternal lineage data) of an almost exclusive contribution of the eastern Anatolian bezoar to the domestic goat gene pool.
Goats are multi-purpose animals that produce milk, meat and fiber and also serve other beneficial roles. In particular, they contribute to the economy of farmers living in arid and semi-arid regions, including southern Iran . Although goat products are often cheaper than sheep products in the market place, goats are favored in the most marginal areas of Iran, where they are easier to manage and better adapted to harsh climate and ecological conditions than sheep. According to the latest livestock census, conducted in 2008, the Iranian caprine population is around 25 300 000 animals (http://faostat.fao.org). Iranian goats are mainly reared in traditional systems by small holders. Since nomadic tribes are almost completely economically dependent on animal rearing, these stakeholders play an important role in the conservation of animal genetic resources, especially of small ruminants.
Genetic diversity is an essential component for population survival, evolution, genetic improvement and adaptation to changing environmental conditions . Information on genetic diversity is therefore necessary to optimize both conservation and strategies for the use of animal genetic resources, to meet future market demands and improved production systems. Molecular tools permit the characterization of genetic resources at the DNA level. Because of favorable characteristics, such as abundant number, high polymorphism and co-dominant inheritance, microsatellite DNA markers have been extensively used for a number of applications in livestock genetics, including parentage testing, breed classification, conservation genetics and also to assess genetic variation and structure within and among populations .
This study was undertaken to examine the pattern of microsatellite variation within and among seven Iranian goat breeds. The resulting information may be used in national plans for sustainable improvement and conservation of goat genetic resources. This research was carried out as part of the IAEA-FAO joint program “Characterization of genetic resources in small ruminants in Asia” (D3.10.25), which aimed at developing methodologies, generating information and formulating decision support systems to analyze phenotypic and molecular genetic diversity, develop microsatellite and related technologies, and enable the development and implementation of national and regional strategies for optimum use and conservation of small ruminants in Asia (http://www-naweb.iaea.org/nafa/about-nafa/index.html).
Iranian goats were sampled in six different areas that extended from the north of Iran, in the Alborz Mountains, south of the Caspian Sea, to the far western border of Iran, in the northern Zagros Mountains, to southern Zagros, along the mountain range. Seven indigenous goat breeds were mainly distributed in six provinces: Gilan, Ardabil, Isfahan, Fars, Kurdistan and Khuzestan (Table 1). A maximum number of five samples per flock were collected from an average of 11 flocks per breed (min = 5, max = 18). Two Pakistani goat breeds, collected from the Punjab province, were also included in the dataset for comparison. The geographic distribution of Iranian breeds sampled in this study is depicted in Figure 1 and their typical phenotype in Figure 2.
Microsatellite DNA analysis
The salting-out method  was used to isolate genomic DNA from blood samples of 317 animals from the seven Iranian and two Pakistani goat breeds. Fourteen microsatellite markers were chosen from the list recommended by the FAO . Forward primers were end-labeled with fluorescent dyes (6-FAM, VIC, NED and PET) [see Additional file 1]. Polymerase Chain Reaction (PCR) was carried out on 50-100 ng of genomic DNA in a 15 μL reaction containing 1.5 μL of 10× PCR buffer, 1 μL of 20 mM dNTPs, 0.2 μL of each primer at a concentration of 10 μM and 1 unit DNA polymerase. Each marker was amplified individually. The “Touchdown” PCR protocol used an initial 5 min denaturation step at 95°C, followed by 3 cycles at 95°C during 45 s and 60°C during 1 min, 3 cycles at 95°C during 45 s and 57°C during 1 min, 3 cycles at 95°C during 45 s and 54°C during 1 min, 3 cycles at 95°C during 45 s and 51°C during 1 min and 20 cycles at 92°C during 45 s and 48°C during 1 min, a 45 s extension step at 72°C and a final 10 min extension step at 72°C. Microsatellite genotypes were visualized with the ABI PRISM 3130XL DNA Analyzer (Applied Biosystems, USA) and alleles were scored using GeneMapper™ software Version 3.7 (Applied Biosystems, USA). In total, 11 samples from four populations (KHL, TAL, MKZ and NAI) were genotyped in duplicate, to evaluate data quality and repeatability. To ensure correct genotype scoring, visual inspection was carried out independently by two experienced operators and all conflicting scores sorted out.
Unbiased estimates of genetic diversity (HE), observed heterozygosity (HO) and mean number of alleles (MNA) were calculated using the Microsatellite Toolkit . The probability-test (exact HW test), used to assess deviations from Hardy-Weinberg equilibrium (HWE) for each locus and population and for loci over all populations, was performed with Genepop 4.0 using a Markov chain of 100 000 steps and 1000 dememorization steps, 500 batches and 10 000 iterations per batch . P values from multiple comparisons were corrected using a Bonferroni correction . Null alleles can decrease estimates of genetic diversity and inflate genetic differentiation . To estimate the potential frequency of null alleles (r) for each locus in each breed, we used the EM algorithm of  in the software FreeNA . This method assumes that deviations from HWE do not result from other causes (e.g. from the Wahlund effect). Values of r ≤ 0.2 are expected not to cause significant problems in the analyses . FSTAT program version 2.9.3  was used to estimate Wright’s fixation indices . Standard errors were generated using the jack-knife method over loci and populations. The same software was used to calculate the inbreeding coefficient (FIS) for each population and a pairwise FST distance matrix. The rarefaction technique of El Mousadik and Petit (1996)  was used in FSTAT to calculate allelic richness (number of alleles in a sample of standardized size). Cervus 3.0 was used to calculate the polymorphism information content (PIC) of each locus . The number of effective alleles (Ne) per locus in each population  was calculated with the POPGENE 1.32 software . The average number of effective migrants exchanged per generation (gene flow, Nm) was calculated with the following formula: Nm = (1 - FST)/(4FST) as applied in Genetix 4.05 . This software was also used to estimate unbiased Nei’s coefficient of gene variation (GST) . Beaumont and Nichols’s approach , implemented in LOSITAN  was used to detect loci under selection. This software uses computer simulation to detect loci for which the genetic diversity within (heterozygosity) and between populations (FST) do not conform to the prediction of a simple island model obtained by coalescent simulations . Similarity in FST/HE values for all loci indicates a shared demographic history. Loci showing unusually large amounts of differentiation may mark regions of the genome that have been subject to directional selection, while loci showing unusually small amounts of differentiation may mark regions of the genome that have been subject to balancing selection . All loci outside a 99.5% confidence interval were removed and the mean FST was calculated again. A final run included all loci. The infinite allele model and 95 000 simulations were used in this calculation. Hierarchical analysis of molecular variance (AMOVA) was conducted with geographical location as grouping factor, using Arlequin 3.11. Goat breeds were spatially partitioned into five groups including (A) TAL and KHL, (B) NAI, TUR and ABD, (C)MKZ, (D)NAJ and (E) TED and PAH (Pakistani populations). AMOVA was used to measure the extent of hierarchical genetic differentiation among the locations, among populations within a location, and among individual within a population .
Three approaches were used to analyze the genetic relationships among individuals and populations: (i) genetic distances and dendrograms, (ii) model-based cluster analysis, and (iii) principal components analysis (PCA). A dendrogram was constructed using the Neighbor-Joining (NJ) algorithm  in DISPAN , based on Nei’s genetic distance (DA) . Trees were edited with MEGA4 . The Bayesian model-based clustering method implemented in STRUCTURE software  was used to investigate population structure and define clusters of individuals on the basis of multi-locus genotypes. The number of assumed clusters (K) varied between 2 and 11. For each K, 10 independent runs were performed with a burn-in of 105 and Markov chain Monte Carlo length of 2 × 105 iterations under an admixture model with correlated allele frequencies and no prior information on the population of origin (popinfo = 0). The assignment probabilities were compiled for multiple runs in the program CLUMPP, which addresses multimodality and/or label-switching in run comparisons . We used the Greedy algorithm to increase computational speed, set the pairwise similarity matrix to G’ and ran 1000 random repeats of the data for the determined value of K. The results of STRUCTURE analyses were depicted using the software Distruct . The estimate of the best K was calculated as described by Evanno et al.  using Structure Harvester v.0.6.92 . PCA was performed using XLSTAT software (Addinsoft, Paris) to summarize and visualize the structure of data described by several quantitative variables, while obtaining the uncorrelated factors between them.
A total of 154 genotypes were produced on the 11 animals, which were genotyped twice. All genotypes were identical between replicates indicating a very high repeatability of the genotyping and scoring procedures adopted. A total of 150 alleles were detected at the 14 microsatellite loci in the nine goat breeds. Allele number ranged from five (MAF035) to 18 (ILSTS029) per locus and the average number was equal to 10.7. Most loci displayed a high degree of polymorphism, as revealed by the PIC values that ranged between 0.435 (INRA0132) and 0.851 (MAF70), with a mean of 0.67. Relevant information per locus, such as the range of allele sizes, location on chromosomes, sequence and label of the primers, number of alleles (observed and effective), PIC and deviation from HWE, are presented in Additional file 1 [see Additional file 1]. Seven out of 14 loci deviated from HWE after sequential Bonferroni correction (P < 0.05). Five populations, NAI, TUR, ABD, MKZ and TED, showed deviation from HWE for at least one locus (P < 0.05). Nine of the 126 locus × population combinations revealed significant departures from HWE. Several estimates of the frequency of null alleles (r) were greater than 0.11, i.e. for BMS1494, ILSTS029 and MAF035 [see Additional file 2]. As previously reported, values of r ≤ 0.2 are expected not to cause significant problems in the analysis . The only exception was the ILSTS029 locus in the ABD population, for which r = 0.23. Frequency distributions of two loci (BM1818 and MCM527) were indicative of balancing selection at the 99.5% probability level, whereas directional selection was suggested at the loci MAF035 and ILSTS029 (Figure 3). Since there was a strong suggestion that the MAF035 locus was under selection because of its low heterozygosity and low number of alleles (Table 1; Figure 3), it was excluded from further analyses. Conversely, to avoid too much loss of information, all other markers were retained, including BM1818 that was reported to flank QTL (Quantitative Trait Loci) for fertility and milk traits in cattle . Based on the raw data for 13 markers, the largest mean number of alleles (MNA) was observed in TUR (8.00) and the smallest in NAJ (4.69). NAJ had the fewest samples (N = 20), but this trend remained consistent for corrected allelic richness (Ar), which was also greatest in TUR (5.99) but smallest in NAJ (4.25). The mean of the effective number of alleles per locus and population ranged from 2.65 (NAJ) to 4.00 (KHL) (Table 1). The mean effective number of alleles had a global mean of 3.50 across all loci, which was remarkably lower than the mean observed number of alleles (10.7). The average unbiased expected heterozygosity over all loci ranged from 0.58 (NAJ) to 0.72 (ABD). The overall mean of gene diversity was equal to 0.671 (Table 1). ABD (0.602) and TAL (0.710) had the lowest and the highest observed heterozygosity, respectively (Table 1). The mean FIS across all loci and breeds was equal to 0.052. FIS were significantly greater than zero in NAI, TUR, ABD (P < 0.001) and MKZ (P < 0.01), which indicated inbreeding in these breeds (Table 1).
The global FST obtained by the jack-knife method over loci was equal to 0.062 ± 0.016 and significantly different from zero (P < 0.001). Wright’s F-statistics were calculated for each of the 13 microsatellite loci across the nine breeds [see Additional file 3]. The AMOVA revealed that most of the molecular variance occurred within breeds (92.90%) while it represented 3.73% among geographic groups and 2.07% among breeds within a geographic group (Table 2). Although small, the group/breed components were both statistically significant (P < 0.01 and P < 0.001, respectively). Very similar results were obtained by clustering breeds according to the groups identified in the STRUCTURE analyses.
Genetic relationship and population structure analysis
The FST values between breed pairs ranged from 0.0041 for the TAL-KHL pair to 0.1622 between PAH and MKZ breeds. Pairwise estimates of FST for north (TAL and KHL) and center (NAI, TUR and ABD) of Iran were not significant. For all remaining breed pairs, FST values were highly significant (P ≤ 0.001) (Table 3).The number of migrants per generation (Nm) ranged from 60.23 between TAL and KHL populations to 1.29 between MKZ and PAH populations (Table 3). A gene flow of 3.62 was obtained between the two Pakistani goat breeds, PAH and TED. Nei’s genetic distances (DA) ranged from 0.0521 between TAL and KHL to 0.2760 between NAJ and PAH. The two Pakistani breeds (TED and PAH) showed the lowest DA distance with TUR. Nei’s standard genetic distance (DS) ranged from 0.0108 (TAL-KHL) to 0.4269 (PAH-MKZ). TED and PAH showed the lowest DS distances with the Iranian KHL and TUR goat breeds [see Additional file 4]. A NJ tree was constructed based on DA genetic distances (Figure 4). Most of the bootstrap values were high (> 70%), which indicated that the dendrogram was very robust. According to the NJ tree, Iranian populations showed a clear clustering, in agreement with the traditional breed classification and geographical origin. An exception was the inclusion in the same cluster of MKZ and NAJ breeds in spite of their rather distant distribution areas, adaptation to different climates and different production purposes. However, they are separated by long branches, indicating that, although they share a common ancestry, they still remain well differentiated. In the DA tree, the two Pakistani breeds (TED and PAH) cluster with the northern Iranian breeds (TAL and KHL), even if the branch lengths indicate that they remain rather distant from the Iranian pool and from each other. PCA grouped Iranian breeds in accordance with the NJ trees (Figure 5A).The first component separates populations according to a northwest to southeast gradient, while the second has no clear geographic component. Northern (TAL and KHL) and central (TUR, NAI and ABD) breeds form two groups, while, the western (MKZ) and southwestern (NAJ) breeds clearly separate from each other and also from the other breeds. Inclusion of the two Pakistani breeds (Figure 5B) does not change the configuration of the Iranian breeds. However, in this case, PAH and TED appear to be closer to the central breeds (NAI and TUR) rather than to the northern ones (TAL and KHL) as in the NJ representations. In the STRUCTURE analysis, K = 5 resulted as the most appropriate number of partitions [See Additional file 5 and Additional file 6]. Analysis at K = 5 divided Iranian goats into three clusters formed by populations from the north (TAL and KHL), center (NAI, TUR and ABD) and west (MKZ and NAJ) of Iran. The two Pakistani breeds were assigned to two other distinct clusters (Figure 6). Admixture was particularly evident between central and northern Iranian clusters, with some components also contributed by the Pakistani PAH breed, whereas the western Iranian cluster formed a quite distinct gene pool.
Iran is close to, and in the case of goats, within, the main south west Asian livestock domestication center. In fact, archaeological remains indicate an early goat domestication (about 10 000 years ago) in the Iranian Zagros Mountains , as well as in the high Euphrates valley and southeastern Anatolia . In addition, analysis of mitochondrial DNA of domestic goats and their wild ancestors (C. aegagrus or bezoar) revealed signals of population expansions in wild populations in southern Zagros (Fars Province) and in the central Iranian Plateau (Yazd and Kerman Provinces), likely indicating a pre-domestication management of wild populations . These regions were therefore hypothesized to be the site of origin of one of the mtDNA haplogroups (the C haplogroup) and of “an incipient goat domestication phase”. However, haplogroup C has a modest frequency within the mitochondrial gene pool of modern goats, thus suggesting that these regions may not have contributed much to the molecular variability of domestic goat maternal lines, and definitely much less than eastern Anatolian sites. Therefore, it is highly interesting to investigate the genetic diversity of Iranian farm animals for several reasons: (1) social and economic role that these livestock play within the country, (2) their geographic area of origin and (3) new information on the domestication processes that may arise from the analysis of the nuclear genome of local genetic resources. Particularly intriguing would be to test if nuclear DNA data agree with mtDNA information in estimating a marginal involvement of Iranian and a major contribution of Anatolian gene pools into goat domestication processes.
The analysis of 13 microsatellites resulted in a mean genetic diversity of 0.671 (Table 1). Although comparison with investigations using different marker sets is only indicative, the value we observed is greater than those reported in Swiss goat breeds (0.51 to 0.58) genotyped at 20 microsatellite loci  and in 11 indigenous south east Asian goats analyzed with 25 microsatellites (0.43-0.60) , but it is slightly lower than those reported in Chinese goat breeds (0.77-0.82) analyzed with six microsatellite loci . However, Di et al.  assessed the genetic diversity of nine Chinese cashmere goats, two Iranian goats and one breed from Guinea Bissau using 14 microsatellite markers and found the greatest diversity among the Iranian breeds.
In the present study, expected heterozygosity and allelic richness had the highest values for ABD and TUR (two Iranian goats) with means of 0.72 and 5.99, respectively. We can thus conclude that Iranian goats possess a remarkably high genetic diversity, as expected, for native populations in the vicinity of a domestication center. With the exception of NAJ, the surveyed populations had higher expected than observed heterozygosities. This resulted in positive FIS values that were highly significant for NAI, TUR, ABD and MKZ (Table 1), indicating some level of inbreeding in these populations. FIS reached the remarkable value of 16.6% in ABD. The presence of null alleles probably contributed to the very high FIS value observed in this breed. Conversely the Wahlund effect, did not contribute to increase FIS values, since no population substructure was detected in the four inbred breeds by Bayesian cluster analysis. The investigated goat breeds showed a remarkable difference between the effective and the observed number of alleles (sometimes a decrease of more than 50%) [see Additional file 7], due to a very low frequency of many alleles across loci. This effect may result from the combined effects of the exchange of migrants between populations and of the post-domestication population expansion that can still be detected in traditionally managed populations nearby domestication sites, as opposed to western breeds that have likely lost many rare alleles by genetic drift during the process of breed formation.
An overall mean of GST = 5.9% (based on 13 markers) indicated that within-breed diversity accounts for a large part of the total genetic diversity of the breeds investigated. This observation is confirmed by AMOVA and by the low average pairwise FST value (0.062; Table 2 and Additional file 3). This value is similar to that i.e. 0.069 reported by Canon et al. but lower than values reported for south-east Asian (0.14; ) and Swiss goat breeds (0.17; ). The test for neutrality (Figure 3 and Additional file 8) suggested that some loci are under directional (MAF035 and ILSTS029) and balancing selection (BM1818 and MCM527). Since a panel comprising only a few microsatellite markers is of limited interest to identify selection signatures, here the test was used merely to decide if certain markers were to be eliminated from population genetic analyses to avoid biased results . However, these loci might merit further investigation, since they are potentially associated to traits of interest, e.g. the MAF035 locus has been associated with a QTL for carcass traits (percent lean in carcass and total fat) in sheep  and BM1818 to QTL for milk and fertility traits in cattle .
The distribution of the Iranian breeds described by PCA (Figure 5A) is consistent with the geographical locations of the farms where samples were collected (Figure 1), confirming at the country level that differentiation of diversity in nuclear genomes of goat breeds contains a significant portion of geographic structure, as has already been reported at the continental level . Interestingly, this geographic structure is maintained also in a system of traditional pastoralism, as that present in Iran, in spite of the fact that nomadism and gene flow exist among the populations. The exchange of migrants among populations is, in fact, relevant (Table 3), in particular among the two northern breeds and, to a lesser extent, among the three central breeds. Gene flow is higher within regions compared to between regions and is particularly low between the west and the other regions, due to the topography of the Zagros mountain range. The results of the STRUCTURE analysis are consistent with this interpretation (Figure 6). Western breeds (MKZ and NAJ) form a defined cluster at K ranging from 2 to the most probable value (K = 5). At this value of K, Iranian breeds from the north (TAL and KHL), center (NAI, TUR and ABD) and west and south-west (MKZ and NAJ) form three clusters. The level of admixture is high in breeds from central Iran. It is lower in northern breeds that appear to contain almost identical proportions of ancestral genomes, confirming their high similarity as indicated by genetic parameters in the NJ tree and PCA. The two Pakistani populations constitute two distinct separate gene pools, although they originate from the same area (Punjab province) in Pakistan. Gene flow between these two breeds (Nm = 3.62) confirms STRUCTURE, NJ and PCA analyses and indicates that PAH and TED are distinct, even if some animals from TED seem to have a large portion of their genome in common with PAH. Overall, NAJ and MKZ, although they share a common ancestry at K = 5 (Figure 6), seem to be quite distinct from each other, as indicated by their DA and DS distances [see Additional file 4], low level of gene flow (Nm = 2.71), long branch length in NJ trees (Figure 4) and clear separation in the PCA plot (Figure 5). In fact, common ancestry does not necessarily imply similarity in gene frequencies. Genetics, geographic distance, agro-climatic conditions, phenotype and main use clearly distinguish these two breeds from each other.
MKZ is a breed of the Kurdish areas (Kurdistan province) of Iran (Figure 1). It is a Mohair-producing breed valued for its shiny fine fiber. It is well adapted to withstand the severe winters that occur in western Zagros, with average daily temperatures below 0°C and heavy snowfalls. A recent report indicates that this breed is presently endangered, due to reduction of population size and number of breeding herds. The population size of MKZ was estimated at over 22 000 animals in 1996, but has progressively decreased to around 5000 heads in 2005 . NAJ, a dairy and fleece goat breed from the Arab region (Khuzestan province; Figure 1) is adapted to extremely high temperatures. Morphologically, MKZ and NAJ are clearly different.
The level of population differentiation and genetic structure observed in Iranian goat breeds are clearly different from that observed in Iranian sheep populations (FST = 0.02, unpublished data). This may be due to the massive amount of gene flow occurring in sheep by translocation of superior breeds over a large geographical distance, because of the higher economic importance of sheep compared to goats. Overall, the degree of differentiation at the few microsatellite marker loci used in this study might appear inadequate to represent the degree of differentiation among breeds that is perceived based on physical appearance and other phenotypic traits. However, “neutral” markers such as microsatellites are designed to reconstruct the evolutionary and demographic history of populations and are theoretically “blind” to the effect of natural and anthropogenic selection that is sometimes very effective and rapid in changing morphological and production traits . It has been reported that degree of differentiation in quantitative traits (QST) typically exceeds that observed in neutral marker genes (FST) , suggesting a prominent role for natural selection in accounting for patterns of quantitative trait differentiation among contemporary populations.
Taken together, all the approaches used to analyze genetic relationship among individuals and populations in this study suggested a high molecular diversity in Iranian goats, with varying levels of genetic distinctiveness among breeds and considerable gene flow between breeds from the same geographic region. Between-breed diversity has a strong geographic component. Iranian goat breeds fall into three clusters, northern, central and western, with the western and southwestern breeds relatively distinct from others. Pakistani breeds show some relationships with Iranian populations, even if their position is not consistent across analyses. Pakistan and Iran are neighbors, connected by the Baluchistan region that is shared by the two countries. There is a long history of contact and mutual influence between the two countries. Agribusiness and livestock exchange have been ongoing for ages, so it is not surprising to find some similarity in the genetic background of Iranian and Pakistani goat breeds.
In conclusion, to maintain the present genetic diversity and structure of these breeds, proper strategies of marker-assisted management need to be designed and implemented. Although a decreasing number of MKZ individuals has been noted, none of these breeds seem to be endangered according to the FAO risk classification system . Therefore, breed management plans should emphasize sustainable use and development, rather than conservation per se. One suggested first step is to organize breeders into formal or informal associations, to facilitate development and implementation of genetic resource management strategies. Inbreeding seems to affect some breeds and an organization of breeders may allow for wider exchange of males within breeds, which would address this problem. Conversely, if the breeders express an interest in maintaining genetic purity, gene flow among breeds and regions should be monitored and avoided. Development of more complex strategies would benefit from the analysis of native breeds with high-density marker panels that can distinguish between neutral and selected genomic regions. This additional information would contribute to the decision making process, in particular by identifying patterns of diversity along genomes of neutral (present day) and selected (very near future) genomic regions. Three out of the seven investigated breeds are reared in geographic areas in which mtDNA provided evidence of early domestication. TUR and ABD (southern Zagros, Fars province) and NAI (central Zagros, Isfahan province) fall exactly in the area in which the C haplogroup is observed at high frequency . Interestingly, these breeds are highly variable (Table 1), are placed on basal short branches in the NJ tree (Figure 4), close to the origin of the PCA plot (Figure 5) and, although highly admixed, quite distinct from those reared on the western side of the Zagros mountain range. These observations reveal the necessity for further investigation of goat nuclear DNA diversity within a much wider geographic context, including Turkey, Europe and Asia. Such an investigation would help to clarify the events that occurred in central Zagros and to the west of the Zagros mountain range during domestication, either confirming or re-discussing the current hypothesis based on maternal lineage data of an almost exclusive contribution of the eastern Anatolian bezoar to the domestic goat gene pool.
Tavakolian J: An Introduction to Genetic Resources of Native Farm Animals in Iran. 2000, Tehran: Animal Science Genetic Research Institute Press
Kumar S, Dixit SP, Verma NK, Singh DK, Pande A, Kumar S, Chander R, Singh LB: Genetic diversity analysis of the Gohilwari breed of Indian goat (Capra hircus) using microsatellite markers. Am J Anim Vet Sci. 2009, 4: 49-57.
Fan B, Han JL, Chen SL, Mburu DN, Hanotte O, Chen QK, Zhao SH, Li K: Individual-breed assignments in caprine populations using microsatellite DNA analysis. Small Ruminant Res. 2008, 75: 154-161. 10.1016/j.smallrumres.2007.09.007.
Sambrook J, Fritsch EF, Maniatis T: Molecular Cloning: A Laboratory Manual. 1989, Cold Spring Harbor: Cold Spring Harbor Press
FAO: Molecular genetic characterization of animal genetic resources. FAO Animal Production and Health Guidelines. 2011, 9-http://www.fao.org/docrep/014/i2413e/i2413e00.pdf,
Park SDE: PhD thesis. Trypanotolerance in West African cattle and the population genetic effects of selection. 2001, University of Dublin
Raymond M, Rousset F: GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism. J Hered. 1995, 86: 248-249.
Rice WR: Analysing tables of statistical tests. Evolution. 1989, 43: 223-225. 10.2307/2409177.
Dakin EE, Avise JC: Microsatellite null alleles in parentage analysis. Heredity. 2004, 93: 504-509. 10.1038/sj.hdy.6800545.
Dempster AP, Laird NM, Rubin DB: Maximum likelihood from incomplete data via the EM algorithm. J R Stat Soc B. 1977, 39: 1-38.
Chapuis MP, Estoup A: Microsatellite null alleles and estimation of population differentiation. Mol Biol Evol. 2007, 24: 621-631.
Goudet J: Fstat: a program to estimate and test gene diversities and fixation indices Version 2.9.3. 2001,http://www2.unil.ch/popgen/softwares/fstat.htm,
Weir BS, Cockerham CC: Estimating F-Statistics for the analysis of population structure. Evolution. 1984, 38: 1358-1370. 10.2307/2408641.
Mousadik A, Petit RJ: High level of genetic differentiation for allelic richness among populations of the argan tree [Argania spinosa (L.) Skeels] endemic to Morocco. Theoret Appl Genet. 1996, 92: 832-839. 10.1007/BF00221895.
Kalinowski ST, Taper ML, Marshall TC: Revising how the computer program CERVUS accommodates genotyping error increases success in paternity assignment. Mol Ecol. 2007, 16: 1099-1106. 10.1111/j.1365-294X.2007.03089.x.
Kimura M, Crow JF: The number of alleles that can be maintained in a finite population. Genetics. 1964, 49: 725-738.
Yeh FC, Yang RC, Boyle T, Ye ZH, Mao JX: POPGENE, the User Friendly Shareware for Population Genetic Analysis. 1999, Molecular Biology and Biotechnology Centre, University of Alberta
Belkhir K, Borsa P, Chikhi L, Raufaste N, Bonhomme F: GENETIX 4.05, logiciel sous Windows TM pour la génétique des populations. 2004, Laboratoire Génome, Populations, Interactions, CNRS Université de Montpellier II,http://www.genetix.univ-montp2.fr/genetix/intro.htm,
Nei M, Chesser RK: Estimation of fixation indices and gene diversities. Ann Hum Genet. 1983, 47: 253-259. 10.1111/j.1469-1809.1983.tb00993.x.
Beaumont MA, Nichols RA: Evaluating loci for use in the genetic analysis of population structure. Proc R Soc B. 1996, 263: 1619-1626. 10.1098/rspb.1996.0237.
Antao T, Lopes A, Lopes RJ, Beja Pereira A, Luikart K: LOSITAN: a workbench to detect molecular adaptation based on a FST-outlier method. BMC Bioinformatics. 2008, 9: 323-10.1186/1471-2105-9-323.
Santalla M, De Ron AM, De La Fuente M: Integration of genome and phenotypic scanning gives evidence of genetic structure in Mesoamerican common bean (Phaseolus vulgaris L.) landraces from the southwest of Europe. Theor Appl Genet. 2010, 120: 1635-1651. 10.1007/s00122-010-1282-0.
Holsinger KE, Weir BS: Genetics in geographically structured populations: defining, estimating and interpreting FST. Nat Rev Genet. 2009, 10: 639-650. 10.1038/nrg2611.
Excoffier L, Laval G, Schneider S: Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol Bioinform Online. 2005, 1: 47-50.
Saitou N, Nei M: The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol. 1987, 4: 406-425.
Ota T: DISPAN: Genetic distance and phylogenetic analysis. 1993, University Park, Pennsylvania State University,http://evolution.genetics.washington.edu/phylip/software.dist.html#DISPAN,
Nei M, Tajima F, Tateno Y: Accuracy of estimated phylogenetic trees from molecular data: II gene frequency data. J Mol Evol. 1983, 19: 153-170. 10.1007/BF02300753.
Tamura K, Dudley J, Nei M, Kumar S: MEGA4: molecular evolutionary genetics analysis (MEGA) software version 4.0. Mol Biol Evol. 2007, 24: 1596-1599. 10.1093/molbev/msm092.
Pritchard JK, Stephens M, Donnelly P: Inference of population structure using multilocus genotype data. Genetics. 2000, 155: 945-959.
Jakobsson M, Rosenberg NA: CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007, 23: 1801-1806. 10.1093/bioinformatics/btm233.
Rosenberg NA: Distruct: a program for the graphical display of population structure. Mol Ecol Notes. 2004, 4: 137-138.
Evanno G, Regnaut S, Goudet J: Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005, 14: 2611-2620. 10.1111/j.1365-294X.2005.02553.x.
Earl DA, Von Holdt BM: STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012, 4: 359-361. 10.1007/s12686-011-9548-7.
Martín-Burriel I, Rodellar C, Cañón J, Cortés O, Dunner S, Landi V, Martínez-Martínez A, Gama LT, Ginja C, Penedo MCT, Sanz A, Zaragoza P, Delgado Bermejo JV: Diversity of Iberian cattle genetic diversity, structure and breed relationships in Iberian cattle. J Anim Sci. 2011, 89: 893-906. 10.2527/jas.2010-3338.
Zeder MA, Hesse B: The initial domestication of goats (Capra hircus) in the Zagros mountains 10,000 years ago. Science. 2000, 287: 2254-2257. 10.1126/science.287.5461.2254.
Peters J, Von Den Driesch A, Helmer D: The upper Euphrates-Tigris basin: cradle of agro-pastoralism?. First Steps of Animal Domestication, New Archaeozoological Approaches. Edited by: Vigne JD, Peters J, Helmer D. 2005, Oxford: Oxbow Books, 96-124.
Naderi S, Rezaei HR, Pompanon F, Blum MGB, Negrini R, Naghash HR, Balkız Ö, Mashkour M, Gaggiotti OE, Ajmone-Marsan P, Kence A, Vigne JD, Taberlet P: The goat domestication process inferred from large-scale mitochondrial DNA analysis of wild and domestic individuals. Proc Natl Acad Sci U S A. 2008, 105: 17659-17664. 10.1073/pnas.0804782105.
Saitbekova N, Gaillard C, Obexer Ruff G, Dolf G: Genetic diversity in Swiss goat breeds based on microsatellite analysis. Anim Genet. 1999, 30: 36-41. 10.1046/j.1365-2052.1999.00429.x.
Barker JSF, Tan SG, Moore SS, Mukherjee TK, Matheson LJ, Selvaraj OS: Genetic variation within and relationships among populations of Asian goats (Capra hircus). J Anim Breed Genet. 2001, 118: 213-233. 10.1046/j.1439-0388.2001.00296.x.
Yang L, Zhao SH, Li K, Peng ZZ, Montgomery GW: Determination of genetic relationships among five indigenous Chinese goat breeds with six microsatellite markers. Anim Genet. 1999, 30: 452-455. 10.1046/j.1365-2052.1999.00548.x.
Di R, Vahidi SMF, Ma YH, He XH, Zhao QJ, Han JL, Guan WJ, Chu MX, Sun W, Pu YP: Microsatellite analysis revealed genetic diversity and population structure among Chinese cashmere goats. Anim Genet. 2011, 4: 428-431.
Canon J, Garcia D, Garcıa-Atance MA, Obexer-Ruff G, Lenstra JA, Ajmone-Marsan P, Dunner S, ECONOGENE Consortium: Geographical partitioning of goat diversity in Europe and the Middle East. Anim Genet. 2006, 37: 327-334. 10.1111/j.1365-2052.2006.01461.x.
Negrini R, D’Andrea M, Crepaldi P, Colli L, Nicoloso L, Guastella AM, Sechi T, Bordonaro S, Ajmone-Marsan P, Pilla F, the Econogene Consortium: Effect of microsatellite outliers on the genetic structure of eight Italian goat breeds. Small Ruminant Res. 2012, 103: 99-107. 10.1016/j.smallrumres.2011.08.006.
Cavanagh CR, Jonas E, Hobbs M, Thomson PC, Raadsma HW: Mapping Quantitative Trait Loci (QTL) in sheep. III. QTL for carcass composition traits derived from CT scans and aligned with a meta-assembly for sheep and cattle carcass QTL. Genet Sel Evol. 2010, 42: 36-10.1186/1297-9686-42-36.
Bahmani HR, Tahmoorespur M, Aslaminejad AA, Abassi MA, Ebnabbasi R: Assessment of demographic, geographical and genetic risk in Markhoz goat population. J Anim Vet Adv. 2011, 10: 162-168.
Leinonen T, O’Hara RB, Cano JM, Merila J: Comparative studies of quantitative trait and neutral marker divergence: a meta-analysis. J Evol Biol. 2008, 21: 1-17.
FAO: In vivo conservation of animal genetic resources. FAO Animal Production and Health Guidelines. 2012, 14-http://www.fao.org/docrep/meeting/026/me879e.pdf,
Merila J, Crnokrak P: Comparison of genetic differentiation at marker loci and quantitative traits. J Evol Biol. 2001, 14: 892-903. 10.1046/j.1420-9101.2001.00348.x.
We wish to thank our collaborators at Deputy of Animal affairs and Veterinary section in Jihad-E-Agriculture organizations of Gilan, Ardabil, West Azarbaijan, Kurdistan, Fars, Isfahan and Khuzestan provinces and many local farmers throughout the areas of distribution of the studied breeds. We also thank all the experts of the International Livestock Research Institute (ILRI) for providing excellent technical assistance. This study was supported by grants from the International Atomic Energy Agency (D3.10.25) and Agricultural Biotechnology Research Institute of Iran-North branch.
The authors declare that they have no competing interests.
All authors contributed in designing the experiment, developing methods and drafting the manuscript. All authors read and approved the final manuscript.
Electronic supplementary material
Additional file 1:Characteristics of 14 microsatellite markers used to study nine Iranian and Pakistani goat breeds. The data provided represent the characteristics of the microsatellite markers used to study seven Iranian and two Pakistani goat breeds. Na indicates the number of alleles at each locus, Ne is the number of effective alleles, PIC is the Polymorphic Information Content calculated by Cervus 3.0.3 software, HE the genetic diversity per locus per population, *P < 0.05, *P < 0.01, ***P < 0.001, NS=not significant. (XLSX 12 KB)
Additional file 4:D S and D A genetic distances between nine goat breeds based on 13 microsatellite markers. The data provided represent the D S and D A genetic distances between the breeds analyzed, based on 13 loci. Standard genetic distances (DS) (below the diagonal) and Nei’s genetic distances (DA) (above the diagonal). (XLSX 9 KB)
Additional file 6:Representation of the number of ideal clusters identified by Structure software. The delta K method (Evanno et al. ) was examined to find the most likely K. (PDF 17 KB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.