Indigenous sheep stands third in position among the ruminants in Bangladesh. They are sparsely distributed throughout the country with a relatively higher concentration in three agro-ecological zones (AEZs) namely Jamuna river basin (JRB), Sundarban delta region (SDR) and Coastal region (COR) (Hassan and Talukdar, 2011). Sheep of Bangladesh are mostly non-descript indigenous type along with a small proportion of crossbreds (Bhuiyan, 2006) and the total population is estimated to be 3.47 million (Department of Livestock Services, 2017). In Bangladesh, the major attention towards sheep rearing mainly for mutton production under small scale mixed farming system. They have the withstand ability of adaptations in hot and humid agro-climatic conditions, capable of bi-annual lambing with multiple births and better resistance to diseases and parasites (Hassan and Talukder, 2011; Hailemariam et al. 2018). Notably, sheep are the only domestic animals can utilize stubbles of cultivated crops, farm waste, small vegetation in harvested fields and fellow lands, and even graze on aquatic weeds which make them suitable for traditional low input production system (Khan et al. 2009; Banerjee et al. 2010). Locally adapted sheep population possesses a lot of genetic variation compared to highly specialized breeds. Despite of many positive attributes possess by this farm animal genetic resource, but very little attention has been paid so far in Bangladesh for their conservation and genetic improvement. Mitochondrial DNA (mtDNA) polymorphisms have proven to be useful for elucidating the domestication history and phylogenetic relationships in many species due to uniparental (maternal) inheritance, no recombination and higher nucleotide substitution rate than nuclear genome (Brown et al. 1992; Xu et al. 2007; Di Lorenzo et al. 2016). More particularly, displacement-loop (D-loop) is the main non-coding regulatory region of mitochondria, also known as control region that controls transcription and replication of mtDNA. The higher variability in the D-loop region compared to other parts of mtDNA facilitates for genetic diversity measures and phylogeographic origin of different species. The study of within and between populations genetic diversity is essential for effective management practices and to develop sustainable conservation strategies (Tanaka et al. 1996; Gaouar et al. 2005; Pariset et al. 2011). The mtDNA sequence polymorphisms have been investigated to address the questions of genetic diversity, maternal origin, domestication history, and genetic distance and/or relationships among sheep breeds or populations across the world (Hiendleder et al. 2002; Tapio et al. 2006; Liu et al. 2016; Ghernouti et al. 2017). Hence, it would be worthwhile to investigate the maternal genetic pool of indigenous sheep of Bangladesh for conservation and effective management plan of this genetic resource. However, mtDNA sequence based study has yet been conducted in sheep populations of Bangladesh. The objectives of this study were the assessment of genetic diversity, phylogenetic relationships, genetic distance and maternal origin of Bangladeshi sheep through utilization of mtDNA sequence information.
MATERIALS AND METHODS
Description of study area and sampling
Blood samples were collected from four different AEZs namely COR (n=26), SDR (n=22), BAT (n=18) and JRB areas (n=24) those possessed relatively higher concentration of indigenous sheep. For sampling, the following Upazilas from the above stated AEZs were included in this study; Subarnachar and Companiganj Upazilas of Noakhali district (COR), Mohadebpur Upazila of Naogaon district (BAT), Shamnagar Upazila of Sathkhira district (SDR), and Bhuapur, Tangail Sadar and Gobindaganj Upazilas of Tangail and Gaibandha districts, respectively (JRB). Information on age, sex, flock size and breed/type of animals were recorded during sampling. In each flock, one or two blood samples were collected from unrelated mature individuals. Blood samples (3-5 mL) were taken aseptically from the jugular vein using venoject tubes containing ethylenediaminetetraacetic acid (EDTA) as anticoagulant. The collected blood samples were kept in icebox and transferred to Animal Genomics and Breeding Laboratory, Department of Animal Breeding and Genetics, Bangladesh Agricultural University (BAU) for further use. Genomic DNA was extracted from the whole blood sample using PrimePrepTM genomic DNA isolation kit (GeNet Bio, South Korea). Concentration and purity of the extracted DNA samples were assessed using NanoDrop spectrophotometer (model ND2000).
PCR and sequencing of mtDNA D-loop fragment
One primer pair was designed based on mtDNA genome sequence information (GenBank accession number: AF039578, Hiendleder et al. 1998) to amplify 482 bp of mtDNA D-loop region using Primer-BLAST tool available in NCBI database. The primer sequence information is as follows; forward: 5’-GTGAAACCAACAACCCGCTT-3’ and Reverse: 5’-TTGAGTATTGAGGGCGGGAT-3’. Polymerase chain reaction (PCR) was carried out in 20µl volume mixtures containing 10× buffer, 10mM dNTP mixtures, 3 U Prime Taq polymerase (GeNet Bio, Korea), 10 µM of each primer and approximately 50 ng of genomic DNA. Biometra T-gradient thermocycler was used to amplify mtDNA D-loop sequence using the following cyclic conditions: initial denaturation at 95 ˚C for 10 min, followed by 37 cycles at 95 ˚C for 30 sec, at 58 ˚C for 30 sec and at 72 ˚C for 45 sec for denaturation, annealing and extension steps, respectively and a final extension at 72 ˚C for 10 min. Visualization of PCR products were performed under ultraviolet light using 2.0% agarose gel stained with green gel and the gel images were documented by digital gel documentation system (GDS-200, Sunil Bio. Inc., South Korea). PCR products were purified using the PrimePrep PCR purification kit (GeNet Bio, South Korea). A total of 62 purified PCR products were finally selected for sequencing of mtDNA D-loop fragment both forward and reverse directions from commercial sequence service provider (Solgent Co. Ltd., South Korea). Sequencing was performed using DNA sequencer 3100 (Applied Biosystems, USA).
The generated raw sequences were retrieved and edited using Bioedit software (Hall, 1999). The intra-population diversity measures such as number of segregating sites, nucleotide diversity, haplotype number and diversity were calculated by using DnaSP v.5.10.01 (Librado and Rozas, 2009). The corrected average pairwise differences, population pairwise difference (FST) and analysis of molecular variance (AMOVA) were calculated using Arlequin v.3.5 (Excoffier and Lischer, 2010) to detect the population genetic differentiation among the Bangladeshi sheep populations. Phylogenetic tree was constructed to determine the evolutionary relationship among the Bangladeshi indigenous sheep including sheep breeds and species of different countries by using MEGA v. 6.0 (Tamura et al. 2013) under UPGMA statistical method and Kimura-2 parameter model. The generated sequences of this study were added with the representative sequences those retrieved from the NCBI database for phylogenic analysis to assess the possible matrilineage of indigenous sheep of Bangladesh. The accession numbers of mtDNA D-loop sequences of standard breeds and indigenous sheep populations (Ovis aries) from different countries are as follows; KP228656.1, KP228655.1, KP228654.1, DQ903236.1, DQ903238.1, DQ903240.1, DQ073050.1, AY829406.1, AY829402.1, AY829403.1, AY829404.1, KU899996.1, KU899997.1, AY829388.1, AY829387.1, AY829386.1, DQ903299.1, DQ903300.1, DQ903302.1, AY829389.1, AY829390.1, AY829391.1, DQ073052.1, DQ073051.1, DQ097452.1, DQ097460.1, DQ097466.1, DQ097450.1, DQ097456.1, DQ097453.1, DQ097462.1, DQ097461.1, DQ097459.1, DQ097451.1, DQ097449.1 and DQ097454.1. The accession numbers of mtDNA D-loop sequences of Ovis sub-species are as follows; KR011779.1 and KR011780.1 (Ovis orientalis) and KR011782.1 and KR011776.1 (Ovis aries musimon).
RESULTS AND DISCUSSION
The mtDNA sequence polymorphism analysis and genetic diversity measures
The PCR amplification generated 482 bp mtDNA D-loop fragments from 62 sheep DNA samples (Figure 1). However, after trimming of 5′ and 3′ end and excluding alignment gap a total of 444 bp were remained for further analysis. The genetic diversity parameters like number of segregating sites, nucleotide diversity, haplotypes number and haplotypes diversity were assessed from 4 different sheep populations. The DNA polymorphisms analysis identified 10 segregating sites that defined 11 haplotypes among the indigenous sheep population of Bangladesh (Figure 2). More than one haplotypes were found in one or more regional Bangladeshi sheep population. However, none of the haplotypes distributed in all the four regional sheep populations. One major haplotype (Hap1) shared by SDR as well as in COR sheep. Higher number of segregating sites generated the highest number of haplotypes in JRB sheep population. In contrast, the lowest numbers of haplotypes were found in COR and SDR sheep samples. The number of haplotypes identified in this study is comparable with the previous findings of Moradi et al. (2017) who reported low haplotype diversity with only 9 haplotypes from 257 individuals of Zel and Lori-Bakhtiari sheep breeds of Iran. However, the number of identified haplotypes in this study are relatively lower than the previous findings reported in Asian, European or African sheep breeds (Chen et al. 2006; Tapio et al. 2006; Gorkhali et al. 2015; Liu et al. 2016; Rafiaand Tarang, 2016; Ghernouti et al. 2017; Ilmira et al. 2018). This differentiation might be attributed with a weakness of selection precision applied on this population in related with artisanal breeding system practices by the breeders. In this study, the haplotype diversity ranged from 0.1052 to 0.6818 and the nucleotide diversity ranged between 0.00044 and 0.0029 (Table 1). Region-specific estimated results showed that JRB population contained the highest segregating sites and nucleotide diversity (0.003) while the lowest nucleotide diversity (0.00044) was found in COR. Othman et al. (2015) stated that three Italian sheep breeds Sarda, Laticauda and Muflon possessed low haplotype and nucleotide diversity that varied from 0.678 to 0.988 and 0.004 to 0.006, respectively and is consistent to this study. The haplotype diversity and nucleotide diversity were estimated to be 0.992 ± 0.010 and 0.019 ± 0.001, respectively in Tibetan sheep breeds (Liu et al. 2016); 0.996 ± 0.003 and 0.0372 ± 0.0001, respectively in Iranian indigenous sheep breeds (Rafiaand Tarang, 2016) and 0.993 ± 0.005 and 0.039 ± 0.001, respectively in Kazakh native sheep breeds (Ilmira et al. 2018). All of the aforementioned values were higher than the present findings that indicate relatively high level of diversity existed in those populations. The haplotypes and nucleotide diversity ranged between 0.883 and 0.970, and varied from 0.009 to 0.029, respectively in indigenous sheep breeds of Nepal (Gorkhali et al. 2015). Singh et al. (2013) reported nucleotide diversity was 0.014 ± 0.007 in Indian sheep breeds and is quite higher than this study. The existence of large number of singletons in Bangladeshi sheep mtDNA sequences suggesting divergent mtDNA sequences possessed by females.
The calculated average pairwise nucleotide difference ranged from 0.044 to 0.549 (Table 2). The highest corrected average pairwise difference was found between COR and JRB sheep (0.549) and the lowest were observed between SDR and JRB (0.052) sheep populations. Population pairwise estimates (FST) showed (above diagonal) COR sheep population had highly significant (P<0.01) distance with three other sheep populations. The pairwise FST values ranged from 0.018 to 0.352 (Table 2).
Figure 1 Image of PCR product (left picture) of mtDNA D-loop region (M: 100 bp DNA marker; S1, S2, S3 and S4: respective samples; NC: negative control) and chromatogram of indigenous sheep mtDNA D-loop sequence of Bangladesh (right picture)
Figure 2 Sequence variations of 11 haplotypes derived from 4 regional indigenous sheep population of Bangladesh
SDR, COR, JRB and BAT denote Sundarban delta region, Coastal region, Jamuna river basin and Barind tract area of Bangladesh
Table 1 Measures of genetic diversity in Bangladeshi indigenous sheep population
N: sample size; S: segregating sites; H: no. of haplotypes; Hd: haplotype diversity and π: nucleotide diversity.
Table 2 Corrected average pairwise difference (below diagonal) and population pairwise FST (above diagonal) estimates among different regional indigenous sheep populations of Bangladesh
* (P<0.05) and ** (P<0.01).
The highest distance was found between JRB and COR populations while the lowest value was observed between SDR and BAT sheep populations. The genetic differentiation (FST) between Iranian sheep breeds Zel and Lori-Bakhtiari was 0.052 ± 0.003 (Moradi et al. 2017) and is comparable to this study. Liu et al. (2016) reported that the FST values ranged from -0.046 to 0.237 among the fifteen Tibetan sheep populations where none of the value varied significantly among the studied populations and almost all FST values of the present study fall within the specified range except between JRB and COR sheep populations. The low genetic distances among JRB, BAT and SDR (0.018 to 0.037) sheep populations of Bangladesh are supported by the findings of Khan et al. (2009) and Deb et al. (2019) who reported the genetic distances varied from 0.0664 to 0.165 in those sheep populations based on microsatellite marker information. Agaviezor et al.(2012) reported low FST ranged between 0.0013 and 0.0033 in Nigerian sheep breeds which are inconsistent with the present results. Luo et al. (2005) reported quite higher genetic variability than this study among Chinese (24.22%) and Mogolian sheep (26.85%) populations. In another investigation, a higher population differentiation was observed among the breeds (29.78%) in Ethiopian indigenous sheep (Adhena, 2018). On the other hand, the mean number of pairwise differences obtained in this study was lower than the reported values of Pereira et al. (2006) and Othman et al. (2015) that ranged between 2.612 and 12.359 in Egyptian, Italian and European sheep breeds. Based on the pairwise nucleotide divergence and population differentiation (FST), coastal sheep can be regarded as separate sheep population from the other three regions sheep of Bangladesh. Analysis of molecular variance (AMOVA) revealed that only 11.16% genetic variation was accounted among the studied four sheep populations of Bangladesh and the other 88.84% variability was associated with within population differentiation (Table 3). Comparing to the present results, even higher genetic variability was found within populations (95.04 to 97.0%) and low variability was distributed among populations (3.0 to 4.94%) and different geographical regions in Iranian and European sheep breeds (Pariset et al. 2011; Rafiaand Tarang, 2016; Moradi et al. 2017). Liu et al. (2016) reported that most of the genetic diversity occurred intra-population variation (98.238%) whereas the remaining came from differences among the studied Tibetan sheep breeds (1.762%). In domesticated animal species, genetic variation attributed to within population is much higher than that of genetic differentiation between populations (Gizaw et al. 2011). Altogether, genetic differentiation among the present studied populations was found low as compared to other standard breeds or varieties. Low differentiation between populations could be attributed to migration of animals to be one of the important factors that lead to natural gene flow from one population to others and thus reducing the differentiation between populations (Kantanen et al. 1995).
To detect the evolutionary relationship of Bangladeshi indigenous sheep with other sheep breeds or populations, two different approaches were employed for phylogenetic tree construction based on reconstructed haplotypes of indigenous sheep of Bangladesh along with Ovis aries as reference (Figure 3) as well as considering haplotypes information from Bangladeshi indigenous sheep, different sheep breeds or populations around the world and also including different sheep species (Figure 4). Phylogenetic analysis revealed the close relationships among the sheep of four different regions without forming any specific clusters or groups for a particular population. More particularly, one major cluster included almost 60% Bangladeshi sheep population data from all considered regions with the Ovis aries reference data and thus pinpointed indigenous sheep are the sole progenitors of Ovis aries (Figure 3). Another cluster showed that the sheep of JRB region is closely related to the BAT or SDR regional sheep population than COR sheep population. Based on microsatellite data, the UPGMA dendogram showed four distinct clusters where JRB and BAT clustered together, and sheep populations of COR, SDR (Garole) and Nagpuri sheep formed separate clusters without significant differences (Khan et al. 2009; Deb et al. 2019) and is supported by the present findings. The absence of any distinct genetic structure among the studied populations suggested strong gene flow due to intermixing among the sheep populations of Bangladesh or common maternal lineages among the regions (Rafia and Tarang, 2016). Moreover, Figure 3 shows indigenous sheep of Bangladesh formed cluster with only Ovis aries sheep species that also depicted their single ancestor. More particularly, it showed a close phylogenetic relationship between most of the Bangladeshi sheep samples and sheep breeds from India (Garole and Bonpala), China (Tibetan) and Turkey (Morkaraman and Tuj). This finding suggested common maternal ancestor of those sheep populations. This is also supported from the facts that domestic sheep (Ovis aries) rapidly spread from the domestication center (Southwestern Asia currently known as Iran, Turkey and Cyprus) to the west following human migrations, to the east particularly towards South Asia and China (Chessa et al. 2009) and supports the matrilineage history of indigenous sheep of Bangladesh.
Table 3 AMOVA analysis of indigenous sheep of Bangladesh among four different regions
Figure 3 Phylogenetic tree of Bangladeshi sheep based on mtDNA D-loop hypervariable region
SDR: Sundarban delta region; JRB: Jomuna river basin; COR: Coastal region and BAT: Barind tract
Figure 4 Phylogenetic tree based on mtDNA D-loop region of Bangladeshi indigenous (SDR: Sundarban delta region; JRB: Jamuna river basin; COR: Coastal region and BAT: Barind tract) sheep breeds along with indigenous sheep populations of different countries and other sub-species (Ovis orientalis and Ovis aries musimon) of sheep
As expected, two sheep species (Ovis orientalis and Ovis aries musimon) along with some other Algerian (Ouled djellal), Mongolian (Texel), Chinese indigenous (Hu, Kazakh and Tashikuergan) and Turkish (Hemsin, Karayaka and Akkaraman) sheep population formed separate clusters that revealed distant relationship with Ovis aries as well as Bangladeshi indigenous sheep.
The analysis of mtDNA D-loop sequence variations depicted low genetic diversity among the indigenous sheep of Bangladesh. AMOVA analysis revealed that 11.16% variability was accounted for population differentiation and 88.84% variation was distributed within the population. The low genetic differentiation among the studied sheep populations might be due to common maternal lineage and strong gene flow among themselves through inter se breeding. Phylogenetic analysis showed that Bangladeshi indigenous sheep are the progenitors of a single ancestor Ovis aries. This is the first attempt on mtDNA D-loop sequence variability analysis in indigenous sheep of Bangladesh. Altogether, this finding provides important information on genetic diversity aspect of Bangladeshi sheep that could be utilized in conservation and future breeding program.
This research was supported by a grant from Bangladesh Academy of Sciences, Bangladesh (Project no.: BAU-USDA PALS LS-19). The authors are thankful to the sheep owners for their cooperation in sheep identification and blood sampling.