According to existing evidence, Iran is one of the main sites of sheep domestication in the world (Zeder, 1999). There are 27 native breeds of sheep in Iran, making this country one of the most diverse areas in terms of genetic resources of sheep in the world. The Baluchi animals account for approximately 30% of all Iranian sheep population raised for meat and wool production. The Baluchi breed is medium size with fat-tailed and the white coat color with black muzzle. The keeping place of the Baluchi breed is the dry and hot climate in the eastern part of Iran (Yazdi et al. 1997). One of the most important quantities in population genetics is the effective population size (Ne) that it is especially important in monitoring endangered populations and species. Ne parameter is the number of individuals of an idealized population that determines the amount of genetic drift, or the rate of inbreeding and genetic diversity in the population under consideration (Crow and Kimura, 1970; Falconer and Mackay, 1996). The effective population size can be defined as both effective population size per generation (Ne) and effective number of breeders (Neb). The artificial selection for more wool, milk and meat production, giving more attention to high-yielding breeds, use of genetically modified animals, extensive use of unplanned crossings lead to decrease sheep native breeds in world (Moradi et al. 2017). Within a population with small effective size, the genetic variation and consequently the genetic gain is limited. Also, the risk of loss of genetic diversity and extinction are more in breeds with small population size. The probability of the two alleles at any locus being identical by descent (IBD) was indicated with a parameter called the inbreeding coefficient (F). Selective methods using the breeding values estimated based on animal model, especially within small populations, can lead to an unfavorable event called inbreeding depression (Falconer and Mackay, 1996). This phenomenon affects growth, production, health and reproduction of inbred animals as well as decreases genetic diversity of population (Vozzi et al. 2007). Therefore, it is necessary to regularly monitor the effective population size and the inbreeding coefficient of a breed to ensure conservation of genetic diversity. Traditionally, the effective population size and inbreeding coefficient have been estimated from the pedigree information on the different breeds of sheep (Tahmoorespur and Sheikhloo, 2011; Ghafouri-Kesbi, 2012; Mokhtari et al. 2015; Gholizadeh and Ghafouri-Kesbi, 2016; Rochus and Johansson, 2017). However, completeness of the pedigree is necessary for obtaining of the Ne and F estimates using pedigree records (Uimari and Tapio, 2011). Recently, genomic data obtained from the high-density SNP chip provide an opportunity for genome-wide estimation of effective population size and inbreeding. The whole genome data can reflect the actual homozygosity percentage of the genome even on incomplete pedigrees. Also, genomic data can consider the relationships of very distant common ancestor. Therefore, genome-based estimates are expected to be more accurate than pedigree-based estimates (Keller et al. 2011). There are four indirect methods for estimating contemporary Ne from genetic data: the temporal (Waples, 1989), linkage disequilibrium (Waples and Do, 2008), allele rarefaction (Hedgecock et al. 2007) and heterozygote-excess (Pudovkin et al. 1996) methods. Because all these methods have a low level of precision and their advantages and disadvantages are different; it is recommended to combine these methods altogether and then integrate the obtained estimates (Zhdanova and Pudovkin, 2008). In recent years, three methods have been suggested for computing inbreeding coefficients using genotype data: 1) FGRMbased on the variance of the additive genotype following VanRaden (2008); 2) FHOMbased on the excess of homozygosity following the method outlined by Wright (1948); and 3) FUNI based on the correlation between uniting gametes following advocated by Wright (1922). Estimations of effective population size have been obtained mostly in dairy cattle. Qanbari et al. (2010) estimated the Ne in German Holstein cattle approximately 103 for four generations ago. Rodríguez-Ramilo et al. (2015) reported that the Ne was 74 in Spanish Holstein. However, few studies have used SNP chip data for estimating the Ne and inbreeding coefficients in sheep. Zhao et al. (2014) estimated the effective population sizes of 7 prior generations in Sunite, German Mutton Merino and Dorper sheep, reporting Ne of approximately 207, 74 and 67, respectively. Moradi et al. (2017) calculated the Ne value in six Iranian sheep breeds from 4 generations, which ranged from 9 to 89. Al‑Mamun et al. (2015) estimated the genomic inbreeding coefficients in five sheep breeds in Australia and did not observe any significant inbreeding in the pure breeds. There are no reports of Ne and inbreeding coefficients based on SNPs chip data in Baluchi sheep. Therefore, the objective of this study was to calculate the genomic estimates of population effective size and inbreeding coefficients in Baluchi sheep using density SNP markers.
MATERIALS AND METHODS
This study has been performed with the approval of Abbasabad Sheep Breeding Station of Mashhad, Iran and Iran National Science Foundation (INSF) (number 89,002,382). All institutional and national guidelines were followed for the care and use of laboratory animals.
Animal resources and genotyping
In this study, we used 96 Baluchi sheep (93 females and 3 males) collected from Abbasabad Sheep Breeding Station, Mashhad, Iran. The conventional industry practices were observed for keeping of animals. The mating period was between August and September and included at most three estrous cycles (51 days). The Illumina OvineSNP50 BeadChip (at the Sci-Life Lab in Uppsala, Sweden), included 54,241 SNPs, was used for genotyping the animals.
Data quality control
Because the inheritance patterns of loci on the X and Y chromosomes are different in male and female animals, first the markers on sex chromosomes were excluded from the analysis. Using PLINK software version 1.07 (Purcell et al. 2007), the quality control criteria including genotyping frequency > 95%, minor allele frequency > 0.05 and Hardy–Weinberg equilibrium P > 0.001 were performed on SNP data.
The expected heterozygosity (He) and the observed heterozygosity (Ho) were calculated for measuring of genetic diversity. Using PLINK software, He and Ho were calculated for each SNP separately and then averaged to find the average genetic diversity for Baluchi breed.
Estimation of effective population size
In the present study, we calculated the Ne in two different time scales: Ne in short-term and in long-term. The short-term Ne measures the contemporary Ne in the current generation and the long-term Ne infers the historical Ne in several generations past on an evolutionary time scale (Waples, 1991; Wang, 2005).
Estimation of the contemporary Ne
The contemporary Ne was calculated using NEESTIMATOR software version 2 (Do et al. 2014) with two different methods: The heterozygote-excess method and the LD method. The effective population size was separately estimated for each chromosome and then averaged to find the average Ne for the Baluchi breed.
The heterozygote-excess method
The observed frequency (Hjobs(i)) and the expected frequency (Hjexp(i)) of heterozygotes having allele i were calculated. The Hjobs(i) was obtained by counting heterozygotes in the population. The Hjexp(i) was calculated following the method of Zhdanova and Pudovkin (2008):
pi: frequency of allele i at each polymorphic locus j.
Nj: number of samples having data at locus j.
The D index for excess or deficiency of heterozygote excess was given by:
The D index was taken as the weighted mean of all Dj(i). Then the effective number of breeders was calculated following Zhdanova and Pudovkin (2008):
The LD method
The default parameters and the random mating model were used for running NEESTIMATOR software using LD method. The estimate of Ne was calculated following the method of Waples and Do (2008):
R2′ called R2′-drift: corrected R2.
S: sample size.
Estimation of the historical Ne
LD values (r2) calculated between each pair of loci following Hill and Robertson (1968). Corbin et al. (2012) found the relationship between LD and Ne that enables the users to include corrections for sample size and uncertainty of the gametic phase. This equation is:
Nt: effective population size t generations ago calculated as t = (2f (ct))-1 (Hayes et al. 2003).
ct: recombination rate defined for a specific physical distance between markers.
α: correction for the occurrence of mutations (Ohta and Kimura, 1971).
r2adj: LD value adjusted for sample size and the gametic phase. r2adj was calculated following the methods of Weir and Hill (1980):
n: number of individual sampled, β = 2 when the gametic phase is known and β = 1 if the phase is not known. In this study, the relationship between the recombination rate (c) and the physical distance (d) was distinguished in according with Sved (1971):
A binning system was implemented in order to obtain averaged r2 values that reflects LD for 0.01 to 10 Mb distances and subsequently obtains Ne in several generations ago. SNeP software version 1.1 (Barbato et al. 2015) was used for estimation of Netrends across generations ago.
Estimation of genomic inbreeding coefficients
In this study, three estimates of inbreeding coefficients (FGRM, FHOM, FUNI) were calculated using SNP chip data. The FGRM estimate was calculated based on the variance of the additive genotypes (VanRaden, 2008). FGRM was computed as follow:
pi: observed fraction of the first allele at locus i, hi = 2pi(1 - pi).
xi: number of copies of the reference allele (Yang et al. 2011).
The FHOM was calculated based on the excess of homozygosity using the method from Wright (1948):
O (hom) and E (# hom): observed and expected numbers of homozygous genotypes in the sample, respectively (Yang et al. 2011).
The FUNI was calculated based on the correlation between uniting gametes in according with Wright (1922):
RESULTS AND DISCUSSION
Data quality control
First, non autosomal SNPs were removed. 213 SNPs for call rates less than 95%, 6865 SNPs for minor allele frequency (MAF) less than 0.05 and 24 SNPs for Hardy–Weinberg (HWE) less than 10-6 were removed. After the quality-control, 41453 SNPs were included in the final analysis.
Analysis of genetic diversity using the average expected and observed heterozygosity was shown in Table 1. The He and Ho estimates were 0.374 and 0.383, respectively; this indicates that there is high genetic diversity in Baluchi sheep. More studies on genetic diversity in sheep have been carried out using microsatellites markers. For instance, Esmaeilkhanian and Banabazi (2006) reported the He within five Iranian sheep populations ranged from 0.744 to 0.847. Vajed Ebrahimi et al. (2017) investigated the genetic diversity in 14 sheep breeds in Iran and reported that the He value was the lowest in Kermani (0.75) and the highest in the Arabi breed (0.88). Microsatellites markers have a high level of polymorphism and heterozygosity, however, these markers are not able to cover the whole genome (Moradi et al. 2017). Currently, with developing the high density markers, it is possible to estimate the genetic diversity in the whole genome. Al-Mamun et al. (2015) used SNPs chip data for investigation of genetic diversity in five sheep breeds in Australia. They found the range of the observed heterozygosity in the pure breeds between 0.3 in Border Leicester (BL) to 0.38 in Merino (MER). Moradi et al. (2017) used genome wide SNP data in some Iranian sheep breeds found the average He and Ho ranged between 0.36-0.37 and 0.37-0.43, respectively. Mohammadi et al. (2017) using SNP chip data in Zandi sheep estimated the average He and Ho as 0.393 and 0.407, respectively. The He and Ho estimates for Baluchi sheep in the present study were consistent with the results obtained by Al-Mamun et al. (2015), Moradi et al. (2017) and Mohammadi et al. (2017). The average of MAF (<0.05) and percent of SNPs with deviation from HWE (<0.05) for Baluchi sheep were obtained less than the studied breeds in research of Moradi et al. (2017).
Effective population size
SNPs information used in estimation of effective population size and genomic inbreeding coefficients was shown in Table 2. The highest and lowest number of examined loci were on chromosomes 1 and 24, respectively. The average percent of polymorphic SNPs (%95.36) was lower than the research done by Al-Mamun et al. (2015). The highest and lowest percent of polymorphic SNPs were obtained on chromosomes 18 and 10, respectively. We used SNPs chip data for estimating Ne in Baluchi sheep. The effective population size estimate for each chromosome is shown in Table 3. We used two methods for estimation of Ne in the current generation. The average Ne estimated based on the LD method was 25, higher than average Ne of 15 estimated based on heterozygote-excess. An explanation for the variation in Ne estimates can be the different nature of the calculation of these methods as mentioned in the material and methods. Heterozygote-excess method estimates the effective number of breeders (Neb), while LD method obtains effective the population size per generation (Ne) (Zhdanova and Pudovkin, 2008). It is reasonable that number of breeders is lower in Heterozygote-excess method.
Table 1 The mean of expected (He) and observed (Ho) heterozygosity, minor allele frequency (MAF) and percent of alleles with deviation from Hardy-Weinberg equilibrium (HWE) in genotyping data after cleaning in Baluchi sheep
Table 2 Summary of single nucleotide polymorphisms (SNPs) information used in estimation of effective population size and genomic inbreeding coefficient in Baluchi sheep
Using pedigree information, Mokhtari et al. (2014) showed that the effective numbers of founders in Iran-Black sheep was 13, while the numbers of founders that contributed to the reference population was 150. Ne estimate based on LD in our study was close to that obtained by Mokhtari et al. (2014) at 28. Most studies about the effective population size in sheep have been conducted based on pedigree information. For instance, Ghafouri-Kesbi (2012) estimated the value of Ne in Afshari sheep as 50. Tahmoorespur and Sheikhloo (2011) estimated Ne in Baluchi sheep equal to 134. Leroy et al. (2013) calculated the average Ne in 40 different sheep breeds approximately 100. However, the pedigree is not always complete and errors may occur in registration that affect Ne estimates. There are a few reports on using of SNPs chip data in estimating the effective population size in sheep. For instance, Mohammadi et al. (2017) estimated the average Neb in Zandi sheep equal to 69. Moradi et al. (2017) calculated the Ne for 4 generations ago ranged from 44 to 89 in Iranian native sheep breeds. The highest and the lowest historically Ne were found in the Zel (89 heads) and Afshari (44 heads) breeds, respectively. Al-Mamun et al. (2015) reported that BL breed had the smallest estimated Ne (140) whereas MER was much more diverse (348). Zhao et al. (2014) calculated the Ne in Sunite, German Mutton Merino and Dorper sheep for 7 generations ago approximately 207, 74 and 67, respectively. These authors included a downward trend in Ne for all breeds due to selection. Liu et al. (2017) estimated the range of Ne on each chromosome in Chinese Merino sheep from 140.36 to 183.33 for five generations in the past. Prieur et al. (2017) with analysis of New Zealand sheep population obtained the value of Ne for five previous generations ranged from 71 to 237 for Texel and Romney breeds, respectively. Burren et al. (2014) reported that the Ne estimates for five generations ago varied from 18 to 31 in seven local Swiss sheep breeds.
Table 3 The estimated values of weighted mean of excess or deficiency of heterozygote excess (D index) and contemporary effective population size with two different methods for each chromosome in Baluchi sheep
Neb: effective number of breeders and LD: linkage disequilibrium.
Figure 1 shows the trend of effective population size changes in Baluchi sheep across last generations. The historical Ne showed a decreasing trend over the previous 1000 generations, with an increasingly steeper slope since about 550 generations ago. In this study, the historical Ne estimate for 4 generations ago was obtained and ranged from 48 to 75.
Figure 1 Trend of effective population size changes in Baluchi sheep across last generations
Our finding of the historical Ne estimate for Baluchi sheep was relatively consistent with the results obtained by Ghafouri-Kesbi (2012), Mohammadi et al. (2017), Moradi et al. (2017), Zhao et al. (2014) (for German Mutton Merino and Dorper breeds) and Prieur et al. (2017). Also, the contemporary Ne estimate in this study was relatively consistent with the results obtained by Mokhtari et al. (2014) and Burren et al. (2014). However, the value of Ne estimated in the present study was lower than the results obtained by Tahmoorespur and Sheikhloo (2011), Leroy et al. (2013), Al-Mamun et al. (2015), Zhao et al. (2014) (for Sunite breed) and Liu et al. (2017). The difference in results could related to the method used in these studies. Using the SNeP software, the user can apply several corrections about account of sample size, mutation, phasing and recombination rate. Thus, it seems that the Ne estimate of Baluchi sheep for 4 generations ago is more reliable, logical and consistent with the estimates of previous research. Zhdanova and Pudovkin (2008) recommended to combine the methods of Ne estimation altogether and then integrate the obtained estimates in order to set a reliable result. Frankham et al. (2014) suggested that the Ne over five generations should be at least 100 for preventing inbreeding depression and Ne in the long-term should be at least 1000 for retaining evolutionary potential.
Table 4 Mean, range (minimum and maximum values), standard deviation (SD) and coefficient of variation (CV) for different estimates
1 Inbreeding coefficient was calculated for each individual and then averaged across the population.
In this study, the estimate of Ne for Baluchi sheep was small reflecting events occurred during the population’s history. Because the Baluchi breed has desirable characteristics for meat production, the selection programs have been more intensively performed over the years.
The average inbreeding coefficients estimated in Baluchi sheep using the SNP genotypes are presented in Table 4. More results of inbreeding coefficient estimates in sheep have been reported based on pedigree information. For instance, Rashidi et al. (2018) estimated the average inbreeding coefficients as 0.85%, 4.54% and 1.22% in the whole population of Baluchi, Iran-Black and Zandi sheep, respectively. Gholizadeh and Ghafouri-Kesbi (2016) reported the mean inbreeding in the whole population of Baluchi sheep 1.60%. Few studies reported the inbreeding coefficients estimates using SNPs chip data. Al-Mamun et al. (2015) reported the average FGRM, FHOM and FUNI ranged from -0.085 to -0.006, -0.084 to 0.002 and -0.084 to 0.002, respectively in five sheep breeds in Australia. Mastrangelo et al. (2018) estimated the mean FGRM and FHOM in local dairy cattle breeds ranged from 0.036 to 0.098 and -0.015 to 0.025, respectively. Zhang et al. (2015) obtained the average FGRM, FHOM and FUNI ranged from -0.062 to 0.345, -0.234 to -0.001 and -0.031 to 0.057, respectively in three cattle breeds. Our estimates of inbreeding coefficients were relatively consistent with the results obtained by Al-Mamun et al. (2015) and Zhang et al. (2015). No significant inbreeding was observed in the studied population. The results of inbreeding coefficients estimated using genomic data are usually different with pedigree-estimated inbreeding coefficients. The assumption in the inbreeding coefficient calculated on the pedigree is that all alleles of founder are different, consequently, the value of inbreeding for all time is positive. Therefore, it is expected that the genomic estimates to be more accurate than pedigree-based estimates.
Our results indicate that the effective population size for Baluchi sheep demonstrates a declining trend. The small Ne causes a reduction in genetic diversity and production performance due to inbreeding depression. Therefore, there should be efforts to avoid further reduction in the Ne. The differences observed between the results can be due to the different nature of calculation for each method. Because SNeP software utilizes sample size, mutation, phasing and recombination rate, Ne estimate of this software (48 to 75 for 4 generations ago) is more reliable in Baluchi sheep.
This study was supported by Iran National Science Foundation (INSF).