Validation study of c.746A>G SNP of BMPR1B gene associated with litter size in sheep genetic resources of Bangladesh

Indigenous sheep is one of the important livestock species in Bangladesh primarily used for mutton production. The total heads of sheep in Bangladesh are about 3.61 million [1] which are mainly concentrated in some agro-ecological zones namely Coastal regions (COR), Barind tracts (BAT) and Jamuna River basin (JRB) areas [2]. Garole (GAR) is a world-famous dwarf type sheep breed available in and around the Sundarban delta region (SDR). Sheep genetic resources of Bangladesh are comprised of mostly non-descript indigenous types along with a small proportion of exotic pure breed and crossbreds [3,4]. Indigenous sheep have greater adaptability under hot and humid climatic conditions with better disease resistance capability. Moreover, the capability of bi-annual lambing with multiple birth and withstand ability to utilize low quality roughage make them a suitable ruminant for small and marginal farmers [5,6]. On the other hand, notable differences were observed among the individuals of within and between sheep populations of Bangladesh for production, reproduction, and morphometric traits. For example, litter size of COR sheep was significantly lower than the sheep populations of SDR and JRB [7] .
Litter size is one of the important traits of multiparous livestock species that is directly related to the economic benefits of farms. Reproductive traits such as ovulation rate and litter size are genetically influenced by major genes known as fecundity genes [8,9]. The BMPR1B or Booroola gene (FecB) was identified as the first major gene that has additive effects on sheep fertility and prolificacy [10,11]. It is located on sheep chromosome 6 and belongs to transforming growth factor beta (TGF-β) superfamily gene [12]. Davis [13] reported that one copy of FecB gene increased ovulation rate by 1.5 and the two homozygous copies by 3.0 in Booroola Merino sheep. Earlier studies reported that mutation(s) of this gene were associated with prolificacy traits in different Indo-Pacific sheep populations. More particularly, the c.746A>G mutation in exon 6 of BMPR1B gene causing non-conservative substitution (Q249R) eventually converts the wild type allele (+) to mutant allele (B) that had significant association with higher prolificacy of ewes [9,14,15]. However, a recent study depicted the prevalence of the aforesaid mutation and its association with litter size in indigenous sheep populations of Bangladesh [16]. Moreover, significant association had also been reported among the locations where JRB population possessed higher prolificacy compared to other regions. Hence, validation study is needed for further confirmation of c.746A>G polymorphism of BMPR1B gene and its association with litter size trait considering all available sheep populations of Bangladesh before applying at the field level as a molecular marker.

Blood sampling and phenotypic data collection
An approval was taken from the university ethical committee (No. 838/BAURES/2020 ESRC/AH/01) for conducting research as per institutional guidelines. Multiparous unrelated mature ewes those who had at least two parity records were included in this study. Blood sampling and litter size information were collected from several sheep populated regions of Bangladesh namely southern coastal region, Sundarban delta region, Jamuna River basin areas and Barind tract. Nomenclature of the respective sheep population was accomplished region basis. Sampling for MUZ and GAR populations were performed from Bangladesh Agricultural University managed and farmers’ flocks as well as from Sundarban delta regions, respectively. Blood samples from two or three individuals were collected from each flock in order to avoid related animals. On the other hand, Bangladesh Livestock Research institute (BLRI), Savar, Dhaka has been conserved all aforementioned sheep populations except Garole. Accordingly, blood sampling and phenotypic data was also collected from the BLRI managed flocks. In total, 192 blood samples (3 to 5 ml) were collected from Jugular vein using vacutainer containing EDTA and were transferred to the Animal Genetics and Genomics lab of Bangladesh Agricultural University, Mymensingh.

DNA extraction and PCR amplification
DNA extraction was performed using AddPrep genomic DNA extraction kit (ADDBIO Inc., Daejeon, South Korea) according to the manufacturer’s instruction. The amount and purity of the extracted DNA was determined using spectrophotometer (Nanodrop 2000, Thermo Fisher Scientific, CA, USA). The primer pair was selected from the published paper by El-Seedy et al.[17] for amplifying the amplicon that harbored the target c.746A>G SNP. FecB fragment of BMPR1B gene was amplified by using following primers with the amplicon size of 190 bp; FecB-forward: 5′-ccagaggacaatagcaaagcaaa-3′ and FecB-reverse: 5′-caagatgttttcatgcctcatcaacaggtc-3′. PCR was carried out using Swift MiniPro thermocycler (ESCO Healthcare, Singapore) in 16µl reaction volume that contains ~100 ng genomic DNA, 10 pmol of each primer and 1x PCR Add Taq master mix (ADDBIO Inc., Daejeon, South Korea). The thermal profile comprised of initial denaturation at 95°C for 10 min followed by 37 cycles of 95℃ for 30 sec., 58℃ for 30 sec., 72℃ for 45 sec. and final extraction at 72℃ for 10 min.

FecB genotyping
FecB genotyping was carried out in order to validate the association between c.746A>G polymorphism of BMPR1B gene with litter size in different sheep populations of Bangladesh. PCR-RFLP method was employed to identify and categorize the sheep populations whether it is a Booroola carrier (BB) or non-carrier (++). The PCR product digestion was performed using AvaII restriction enzyme (New England Biolabs Inc., CA, USA). Digestion was done by mixing 10 µl of PCR products with 5 U (10 U/µl) restriction enzyme and then incubated for 10 minutes at 37℃. Visualization of the digested product was done by gel electrophoresis in 2.5% agarose gel stained with 5 µl (10 µg/µl) of safe gel stain dye (green) and the images were captured by GDS 200 digital gel documentation system (Sunil Bio Inc., Seoul, South Korea).

Statistical analysis
Genotypic and allelic frequencies were calculated according to Falconer and Mackay [18]. Single marker association analysis was carried out to evaluate the relationships between the resultant genotypes of BMPR1B gene polymorphism and litter size trait using agricolae package in R [19]. Mean separation was tested by pastecs package in R [20]. The management system and population were considered as fixed effect on litter size where the effects were calculated according to the following model.
Y_ijkl = µ+ M_i +P_j + G_k +I_l +e_ijkl
Where, Y_ijkl = the dependent variable (litter size); µ = the overall mean; M_i = the fixed effect of i^th management system; P_j = the fixed effect of j^th population, G_k = the effect of k^th genotype, I_l = interaction effects and e_ijkl = the residual error.

PCR-RFLP genotyping of BMPR1B gene fragment
PCR–RFLP genotyping was performed for 192 individuals from five different sheep populations of Bangladesh to categorize them as wild (++), carrier (+B) or mutant (BB) types based on c.746A>G polymorphism of BMPR1B gene. AvaII cleavage site (G|GACC) created three different banding patterns ++ (190 bp), +B (190, 160 and 30 bp) and BB (160 and 30 bp) in agarose gel electrophoresis. The PCR-RFLP genotyping of BMPR1B (FecB) gene fragment detected all three genotypes in sheep populations of Bangladesh (Figure 1). The high prolific Booroola mutant digested to yield a 160 bp whereas non carrier products remain uncut (190 bp).

Genotype and allele frequencies of FecB mutation of BMPR1B gene
Table 1 describes the genotypic and allelic frequencies of FecB genotypes in different sheep populations of Bangladesh. Wild type genotype (FecB⁺⁺) was dominant in COR (60%) and MUZ (95%) ewes whereas the proportion of mutant BB ewes were found negligible (<2%). Frequency of the FecB⁺⁺, FecB^+B and FecB^BB genotypes in prolific JRB were 0.09, 0.25 and 0.66 respectively. Similarly, 50% of BAT ewes detected as FecB^BB genotype and it was 55% in case of GAR sheep population. Among the investigated sheep samples, the prolificacy attributed BB genotype was highest in JRB (66%) while it was absent in MUZ sheep population. Accordingly, the FecB^B allele frequency was predominant in JRB ewes (78%) and was lowest in MUZ individuals (2%). Considering all sheep populations, 35% (67) individuals were non-carrier having wild type (++) genotype and the frequency of heterozygous (FecB^+B) and homozygous carrier (FecB^BB) were 0.33 and 0.32, respectively.

Effect of c.746A>G SNP genotypes on litter size trait
Table 2 represents the effects of fixed factors and their interactions with litter size trait. There were highly significant (P<0.001) association observed between litter size trait with population, genotype and farming systems. On the contrary, none of their interactions had significant effects on litter size trait (P>0.05) in the studied populations. The effects of the genotypes on mean litter size performance among the studied sheep populations are presented in Table 3. Litter size performance among the individuals of three genotypes (FecB⁺⁺, FecB^+B, FecB^BB) were insignificant (P>0.05) in COR, MUZ and GAR sheep populations. However, in JRB sheep population, litter size differed significantly (P<0.05) among the homozygous BB (2.50±0.17), heterozygous +B (1.5±0.19) and homozygous non-carrier ++ (1.67±0.33) genotypes. Similarly, significant associations of genotypes with litter size trait were also observed in BAT ewes where the mean litter sizes of FecB⁺⁺, FecB^+B and FecB^BB genotypes were 1.33±0.33, 1.79±0.10 and 2.06±0.06, respectively. Altogether, the FecB^BB genotype had 0.67 and 1.03 more lambs than FecB^+B and FecB⁺⁺ genotypes. Among the FecB^BB genotypes, the highest litter size was observed in JRB population. Population wise mean litter size is presented in Table 4. Comparatively better litter size was observed in JRB (2.17±0.15), GAR (1.90±0.10) and BAT (1.88±0.07) sheep populations whereas COR and MUZ ewes had lower litter size (1.12±0.03 and 1.02±0.02). The overall litter size was found 1.59±0.05 that differed significantly among those populations (P<0.001). JRB population had the highest litter size and MUZ occupied the lowest position among the aforementioned five populations.

Effect of management system on litter size trait
The mean litter size in the studied sheep populations on the basis of management system (farmer’s and on station level) is given in Table 5. Even though management system had insignificant (P>0.05) effects in COR and MUZ sheep populations, litter size differed significantly between two production systems (P<0.05) in BAT (1.58±0.13 versus 2.05±0.05) and JRB (1.68±0.19 versus 2.43±0.18) populations. Moreover, litter size was comparatively better in farmer’s condition than on-station management system for these two populations. On the other hand, unavailability of litter size data for GAR sheep at institutional flocks limits the comparison between two management systems.

Reproductive efficiency of multiparous animal is primarily determined by prolificacy attributing traits where sheep is considered as a model animal for analyzing fecundity genes. Polymorphisms in the several major genes affecting litter size traits have been investigated in different sheep populations worldwide. The FecB was the first identified major gene affecting prolificacy and also had positive effects on ovulation rate and litter size in Booroola Merino ewes [21]. Fabre et al. [22] reported that both ovulation rate and litter size was increased with each copy of FecB mutation in BMPR1B gene. Earlier studies reported the FecB mutation of BMPR1B gene in Australian Booroola Merino [9,13], Indian Garole and Nilagiri sheep [11,23], Indonesian Javanese [24], Chinese Small Tailed Han and Hu sheep [24,25] and Iranian Baluchi and Kalehkoohi sheep [26,27]. Taken together, FecB mutation is predominant in different Asian sheep breeds and following this prolificacy attributed mutation was also detected in the sheep populations of Bangladesh and thus supports those previous findings.
The resultant genotypes of c.746A>G SNP of BMPR1B gene are almost equally distributed in in Bangladeshi sheep populations where 32% of ewes were prolificacy attributed mutant (FecB^BB) type and is quite lower than the findings of Chu et al. [15] who reported the frequencies of FecB^BB genotype were 56 and 68% in Chinese Small tail Han and Hu sheep, respectively. However, FecB^BB genotype frequencies varied from 6 to 13% in Indonesian fat-tailed [28], Iranian Kalehkoohi [27] and Luzhong mutton sheep [29] and is lower than the present findings. No wild type ewes (FecB⁺⁺) were detected by Polley et al. [11] in prolific Indian GAR sheep and likewise only 7% ewes with wild type genotypes were also found in Bangladeshi GAR population. In Luzhong mutton sheep only 5.2% ewes were mutant (FecB^BB) [29] which is comparable with the  frequency of mutant ewes in COR sheep of this study (2%). Altogether, genotype and allele frequency are population or breed specific that varies largely due to genomic structure of the population, sample size and random genetic drift.
Like Bangladeshi sheep populations, there are different sheep breeds or populations worldwide reported significant association with FecB gene mutation. More particularly, Chinese Small-tailed Han sheep and Hu sheep [15], Luzhong mutton ewes [29], Iranian Baluchi and Kalehkoohi sheep [26,27], Indian Garole × Malpura crossbred sheep [30] as well as Mongolian sheep breeds [31] had significant association with FecB mutation of BMPR1B gene and are in agreement with the present study. GAR is a highly prolific sheep breed both in India and Bangladesh having coexisting mutations in two different prolificacy attributed genes BMPR1B and GDF9 [11] that support the present findings. According to Davis et al. [32] FecB mutation has been fixed in some GAR populations and homozygous carrier genotype (FecB^BB) is the original genotype of this breed. It was also reported that Booroola gene was introduced into Australian Booroola Merino through GAR sheep [32]. However, in Indian Bonpala sheep showed non-significant association between genotypes and litter size [33] that contradicts to this study. In addition, the FecB mutation was absent in Thoka, Woodlands, Olkuska, Lacaune, Belclare, and Cambridge sheep breeds [32] and five other Egyptian sheep breeds (Rahmani, Ossimi, Awassi, Barki and Awassi x Barki crossbred) [34] that support the current data of Coastal and Muzaffarnagari sheep populations where incidence of FecB mutation was negligible.
Due to additive effect of FecB mutation on sheep litter size, carrier genotype FecB^BB and FecB^B+ of BMPR1B gene in Small Tailed Han and Hu sheep gave 0.78 and 0.58 more lambs (P<0.01) than those with wild type FecB⁺⁺ genotype [15]. Similar results were also found in three regional sheep populations JRB, BAT and GAR sheep of this study where mutant ewes (FecB^BB) resulted in 0.83, 0.73 and 0.67 more lambs than their wild type counterparts (FecB⁺⁺). Likewise, Mishra et al. [30] found 65.6% higher litter size in ewes with mutant FecB^BB genotype compared to non-carriers (FecB⁺⁺) in Indian Garole × Malpura crossbred sheep. Mahdavi et al. [27] reported that the homozygous FecB^BB and carrier FecB^B+ genotypes possessed 0.52 and 0.35 more lambs compared to homozygous wild-type (FecB⁺⁺) in Iranian Kalehkoohi sheep. Altogether, the above stated findings are in agreement with this study having various degrees effects of mutant genotypes on litter size trait in sheep.
This FecB mutation exert effects not only on litter size trait but also reported to have effects on other traits related with reproductive efficiency of ewes such as fertility and embryo survival, lamb survival, lamb growth, carcass quality, wool production that supports the present findings [35]. Moreover, the g.29380965A>G locus of BMPR1B gene found to had negative effect on the litter size of STH sheep of China, but combination with FecB genotype produce significantly (P<0.05) higher litter size [36]. On the other hand, some non-genetic factors were reported to have significant effect on litter size. Chu et al. [37] observed significant (P<0.05) effect of lambing season and parity as non-genetic factors and it was also stated that improved ovulation rate may accelerate by nutrition and management system [35]. Among the Bangladeshi sheep populations, the litter size was significantly better in farmers level (P<0.01) for JRB and BAT populations that revealed precise selection of ewes at farmers’ level for prolificacy traits. Taken together, the present findings revealed strong evidence for FecB mutation and its association with litter size trait in the regional sheep populations of Bangladesh. However, the effects of this mutation were close to zero in Coastal and Muzaffarnagari sheep due to low frequency of mutant alleles.

The present study validated that c.746A>G mutation of BMPR1B gene had significant association with litter size in the sheep populations of Bangladesh. Population wise analysis depicted JRB sheep as the most prolific sheep followed by GAR and BAT sheep populations due to the additive effect of FecB mutation. This mutation could be used as a molecular marker for selecting prolificacy attributed FecBBB genotype even at earlier stage of sheep and thereby possible to employ molecular information guided breeding against the long term conventional selective breeding to accelerate genetic gain particularly for low heritable trait like litter size.

This research was supported by a grant from BAS-UDSA (Project no. BAS-USDA PALS LS-19). The authors are thankful to Goat and Sheep Production Research Division, Bangladesh Livestock Research Institute (BLRI) for their cordial support in blood sampling and phenotypic data collection. The authors are also thankful to the sheep owners for their cooperation in sheep identification and blood sampling.

SIM and MSAB designed of the experiments. TMM, SIM and MHP performed the experiments. NHD and SA were involved in phenotypic data recording. MMH and MHP analyzed data. SIM and TMM contributed to drafting the article. MMH and MSAB contributed to revising it critically for important intellectual content and made the final approval of the version to be published.

There is no conflict of interest among the authors.