The study was conducted after ethical approval from the institution Animal Care and Use Committee (ACUC) of the University of Veterinary Sciences, Pakistan. Phenotypic data was collected from the lambing records of three sheep farms i.e., Small Ruminants Training and Research Centre (SRT & RC), UVAS, Livestock Experiment Station Allah Dad, Khanewal, and Livestock Production Research Institute, Bahadur Nagar, Okara, Pakistan. All the three flocks were on semi-intensive production system and homogeneous age and parity sample was selected. Across the farms housing 680 Lohi sheep, comprehensive data on lambing, ewe age at lambing, production systems, and birth type (single or multiple) were recorded. Out of this population, 200 ewes were specifically chosen for molecular analysis, with blood samples collected accordingly. These ewes were subsequently categorized into two distinct groups based on their lambing history: Group A comprised ewes that had previously given birth to singles (n=100), while Group B consisted of ewes that had previously delivered twins (n=100). The selection aimed for a homogeneous sample in terms of age and parity, which is easier to achieve within a single breed like Lohi.

The genetic variants in B4GALNT2 (rs410464838) gene were analyzed and were checked for their association with litter size. A total of 5 ml of blood sample was collected out of each ewe in 0.5 M EDTA (pH: 8) and were stored at -20°C. DNA was extracted using a modified phenol-chloroform protocol [13,14]. Primer pairs for the amplification of the genomic region were designed using Primer 3 software [14] and B4GALNT2 was genotyped by PCR-RFLP method. Forward primer 5’ – TCTTAAATTGCCCTCAGTCC -3’and reverse primer 5’- AAGAACTCTGCCAAGAGAGG -3’; were used for the study. The PCR was done in 25 μl reaction volume containing 50 ng DNA template, 10 pmol forward and reverse primers, 0.4 mM each dNTP, 2.5 mM MgCl₂, 1X Taq buffer and 1 U of Taq DNA polymerase (Thermo Scientific). It was then subjected to thermocycling conditions consisting of 4 minutes of denaturing at 94°C followed by 30 cycles of amplification each containing three steps. The denaturation was carried out at 94°C for 30 seconds, annealing at 58°C for 30 seconds and an extension at 72°C for 45 seconds and lastly final extension of 5 minutes at 72°C. The amplicons formed as a result of PCR were analyzed on 1.5% agarose gel.

PCR-RFLP method was used for genotyping of B4GALNT2 SNPs. After amplification, the 240-base pair (bp) PCR product was digested with restriction enzyme AvaII (Thermo Scientific). Samples were then incubated at 37°C for 16 hours for complete digestion. Digestion of the B4GALNT2 allele yielded fragments of 240 bp for GG genotype, 240, 163 and 77 bp for AG genotype, and 163, and 77 bp for AA genotype. The digested products were separated by 4% agarose gel electrophoresis. Wild type GG and variant alleles were differentiated from each other based upon the different migration profiles of digested DNA fragments.

Selected DNA samples were sequenced to validate the results of RFLP. They were later used as controls in analyzing the RFLP data. DNA sequencing was done commercially by automated ABI PRISM® 3730 Genetic Analyzer. The sequencing results were later analyzed using Chromas software. The intronic variant was analyzed using Ensembl Variant Effect Predictor (VEP) and Global minor allele frequency (GMAF).

The data was entered and analyzed using SAS software. Genotypic and allelic frequencies were calculated. Having estimated the allele frequency in each group, the expected frequencies of the AA, AG and GG genotypes were calculated under the assumption of Hardy-Weinberg Equilibrium (HWE). The data were entered and analyzed using SAS software, where genotypic and allelic frequencies were calculated to assess the genetic diversity within the populations. After estimating the allele frequency in each group, the expected frequencies of the AA, AG, and GG genotypes were computed under the assumption of Hardy-Weinberg Equilibrium (HWE). This assumption is fundamental in population genetics as it provides a baseline to evaluate whether observed genotype frequencies deviate from expected distributions due to evolutionary forces.

Chi-square tests were used to find the association between genotypes and litter size, allowing for a clear determination of whether specific genotypes are linked to reproductive outcomes. Observed and expected frequencies were then compared using an exact goodness of fit test using the Xmulti function in the XNomial package in R. This approach is particularly useful for small sample sizes or when expected frequencies are low, ensuring robust results. This was done separately for twin and single lambing data, and this allowed assessment if each group was in HWE. In addition, observed genotype frequencies between twin and single samples were compared using a Fisher's exact test, through the Fisher test function in R. This method is advantageous for analyzing categorical data with small sample sizes, providing precise p-values for evaluating associations between genotype distribution and lambing outcomes.

To check for false positives in the association between the B4GALNT2 genotypes and lambing outcomes, a Bayesian logistic regression model was employed. This approach allowed for the incorporation of prior knowledge about the relationship between alleles and reproductive traits, enhancing the robustness of the results. The model utilized Markov Chain Monte Carlo (MCMC) methods to estimate posterior distributions for the parameters, providing a probabilistic framework to assess the likelihood of true associations while controlling for false positives. Together, these statistical methods offer a comprehensive framework for understanding the genetic influences on litter size in Lohi sheep. 

Results

Figures & Tables

Discussion

Ovulation rate is the major source of variation governing multiple births in sheep. Genetic variants in B4GALNT2 affect the ovulation rate [10] and identification of such markers can help in devising a successful selection program for litter size in sheep. Although there are genetic variations reported in different genes with major effects on reproductive performance traits, such as ovulation rate and litter size in different sheep breeds around the world [6,15-16], the effect of environmental factors like age, parity, nutrition, body weight and season cannot be denied. Adequate nutrition is essential for optimal reproductive health, as deficiencies can lead to poor body condition and reduced fertility, while a balanced diet enhances ovulation and litter sizes. Environmental factors, such as temperature and housing conditions, can also affect reproductive performance by causing stress that disrupts hormonal balance. Estimates of heritability regarding prolificacy in sheep varies across different sheep breeds but it is commonly regarded as a complex trait with many genetic loci and environmental factors influencing the litter size [17]. In studies where large heritability was found for the mean prolificacy, with high percentage of twins at lambing environmental variability was greatly reduced [18]. So, we reduced the non-genetic factors as possible to investigate the association of genetic variants in the B4GALNT2 gene with litter size in Lohi sheep. Lohi is a local mutton breed, as for a successful improvement program, it is essential to evaluate indigenous breeds as well [19].

Our study reports a mean twinning percentage of 17.05% while the highest twinning was 24.5% in one of the farms. Previously, the average twinning percentage in Lohi sheep was reported to be 33% by Ahmad and colleagues. This might be due to the genetic causes that remain to be validated. In order to correlate this phenotype with molecular cause, allelic and genotypic frequencies of B4GALNT2 were estimated and the association of the litter size with B4GALNT2 was determined.

A difference of genotype and allele frequencies was observed between ewes with twin births and those with single births. B4GALNT2 genetic variant (rs410464838) resides on chromosome 11 sheep and is biallelic (A/G). This variation is present in the flanking region of the last exon of B4GALNT2. This is a transcript variant occurring within an intron of a coding transcript near the splice site. As mutations in B4GALNT2 enhance gene expression that can result in multiple ovulation, this intronic variant can affect transcription levels, and transcript stability [20,21]. Moreover, intronic variants can also increase the efficiency of mRNA translation [22,23]. However, functional studies can assess its importance in gene regulation. Allele ‘A’ was markedly frequent in animals that gave birth to twins while ‘G’ allele was common in the other group. In animals that gave birth to single animals, GG genotype was the most frequent with a frequency of 0.55 and a corresponding G allele frequency of 0.68 (and A allele frequency of 0.32). In contrast, in the group of ewes that gave birth to twins, AA genotype was the most common (genotype frequency of 0.89) with a corresponding allele frequency of 0.94; there were no homozygous GG sheep observed. The results showed that the potential association of A allele with twining in Lohi sheep. In addition to this, observed genotype frequencies for twin and single samples were compared using a Fisher exact test in R. This comprehensive analytical approach aimed to discern potential associations between the B4GALNT2 variant genotypes and litter size, ensuring a robust examination of genetic influences in both twin and single datasets. The Value of chi-square =31.403 and p<0.0001 showed signification association between genotypes and litter size.

In our study, we noticed that the single birth frequencies were not in Hardy-Weinberg equilibrium (HWE) that can be attributed to several potential factors. Firstly, HWE assumes a large, randomly mating population with no evolutionary influences such as natural selection, genetic drift, mutation, or gene flow. If any of these conditions are violated, deviations from HWE can occur. For example, if there is non-random mating—such as inbreeding or assortative mating, this can lead to an increase in homozygosity and affect allele frequencies. However, due to the difference of genotypic and allelic frequencies in both groups ‘A allele’ of B4GALNT2 gene can be considered as a potential marker for fecundity in Lohi sheep. However, functional studies and testing of this polymorphism in large samples are required to validate the results. Previously mutation in B4GALNT2 has been found to be strongly associated with litter size in ten breeds of sheep [24,26] However, this study was conducted on 200 animals (100 of each group) and therefore studies on larger groups of animals are required to validate the results of the study and investigate the effect across different parities. It is also important to validate this association on different populations of Lohi sheep before implementation in an MAS program. From an operational point of view, there are many challenges that need to be overcome to incorporate this marker as part of a MAS program. Likely, it would be part of a selection program that would select for more than one gene and more than one trait. Further, any selection would need to be made to ensure that the resultant ‘improved breed’ is appropriate for the environment in which it is going to be raised.

Future research should focus on conducting functional studies to confirm the specific biological mechanisms by which B4GALNT2 influences twinning rates in sheep. This could involve gene knockout or overexpression experiments to observe the effects on fertility. Additionally, expanding the investigation across diverse sheep populations will help assess the prevalence and impact of different alleles in various environmental contexts. Longitudinal studies tracking reproductive performance over multiple breeding seasons could further elucidate the long-term effects of B4GALNT2 variants on litter size, ultimately contributing to more effective genetic selection strategies. 

Author Contributions

Haiba Kaul conceptualized and supervised the study, Asad Ali and Afzal Rashid were also involved in supervision, Fouzia Razaq, collected data and performed experiments, Bilal Bin Majeed helped in data analysis and DNA sequencing and Shagufta Naz helped in statistical analysis and was involved in proofreading the manuscript.

Acknowledgement

The authors declare that there is no conflict of interest regarding the publication of this paper.

References

1. Haile A, Getachew T, Mirkena T, Duguma G, Gizaw S, et al. Community-based sheep breeding programs generated substantial genetic gains and socioeconomic benefits. Animal, (2020); 14(7):1362-1370.
2. Thornton PK. Livestock production: recent trends, future prospects. Philosophical Transactions of the Royal Society B: Biological Sciences, (2010); 365(1554): 2853-2867.
3. Zhou M, Pan Z, Cao X, Guo X, He X, et al. Single nucleotide polymorphisms in the HIRA gene affect litter size in Small Tail Han Sheep. Animals, (2018); 8(5): 71
4. Khazaal NM, Alghetaa HF, Al-Shuhaib MB, Al-Thuwaini TM, Alkhammas AH. A novel deleterious oxytocin variant is associated with the lower twinning ratio in Awassi ewes. Animal Biotechnology, (2023); 34(8): 3404-3415.
5. El Fiky ZA, Hassan GM, Nassar MI. Genetic polymorphism of growth differentiation factor 9 (GDF9) gene related to fecundity in two Egyptian sheep breeds. Journal of Assisted Reproduction and Genetics. (2017); 34(12): 1683-1690.
6. Zongsheng ZH, Heng YA, Yaosheng YU, Yifan XI, Manjun ZH, Liang H. Identification of SNP within the sheep RXRG gene and its relationship with twinning trait in sheep. Kafkas Universitesi Veteriner Fakultesi Dergisi, (2018); 24(1): 39-43.
7. Jia JL, Xie DJ, Zhang YY, Zhang HX, Zhang LP, et al. Association of BMPR-1B gene 3"-UTR region polymorphism with litter size in Tibetan sheep. Kafkas Universitesi Veteriner Fakultesi Dergisi, (2021); 27(3): 285-289.
8. Chantepie L, Bodin L, Sarry J, Woloszyn F, Ruesche J, et al. Presence of causative mutations affecting prolificacy in the Noire du Velay and Mouton Vendéen sheep breeds. Livestock Science, (2018); 216: 44-50.
9. Drouilhet L, Lecerf F, Bodin L, Fabre S, Mulsant P. Fine mapping of the FecL locus influencing prolificacy in Lacaune sheep. Animal Genetics, (2009); 40: 804–812.
10. Drouilhet L, Mansanet C, Sarry J, Tabet K, Bardou P, et al. The highly prolific phenotype of Lacaune sheep is associated with an ectopic expression of the B4GALNT2 gene within the ovary. PLoS Genetics, (2013); 9(9): e1003809.
11. Ahmad Z, Yaqoob M, Younas M. The Lohi sheep: a meat breed of Pakistan review. Pakistan Journal of Agricultural Sciences, (2001); 38(3-4): 69-72.
12. Grimberg J, Nawoschik S, Belluscio L, McKee R, Turck A, Eisenberg A. A simple and efficient non-organic procedure for the isolation of genomic DNA from blood. Nucleic acids research, (1989); 17(20): 8390.
13. Kaul H, Riazuddin SA, Shahid M, Kousar S, Butt NH, et al. Autosomal recessive congenital cataract linked to EPHA2 in a consanguineous Pakistani family. Molecular vision, (2010); 16(58-59): 511-517
14. Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, et al. Primer3–new capabilities and interfaces. Nucleic Acids Res, (2012); 40(15): e115.
15. Xu SS, Gao L, Xie XL, Ren YL, Shen ZQ, et al. Genome-wide association analyses highlight the potential for different genetic mechanisms for litter size among sheep breeds. Frontiers in genetics, (2018); 10(9): 118.
16. Akhatayeva Z, Mao C, Jiang F, Pan C, Lin C, et al. Indel variants within the PRL and GHR genes associated with sheep litter size. Reproduction in Domestic Animals, (2020); 55(11): 1470-1478.
17. Murphy TW, Keele JW, Freking BA. Genetic and nongenetic factors influencing ewe prolificacy and lamb body weight in a closed Romanov flock. Journal of Animal Science, (2020); 98(9): skaa283.
18. SanCristobal-Gaudy M, Bodin L, Elsen JM, Chevalet C. Genetic components of litter size variability in sheep. Genetics Selection Evolution, (2001); 33(3): 249-271.
19. Sharif N, Ali A, Dawood M, Khan MI, Do DN. Environmental effects and genetic parameters for growth traits of Lohi sheep. Animals, (2022); 12(24): 3590.
20. Xu L, Chen Z, Chen S, Chen Y, Guo J, et al. An Identification of Functional Genetic Variants in B4GALNT2 Gene and Their Association with Growth Traits in Goats. Genes, (2024); 15(3): 330.
21. Park JS, Woo SJ, Song CS, Han JY. Modification of surface glycan by expression of beta-1, 4-N-acetyl-galactosaminyltransferase (B4GALNT2) confers resistance to multiple viruses infection in chicken fibroblast cell. Frontiers in Veterinary Science, (2023);10: 1160600.
22. Rogalska ME, Vivori C, Valcárcel J. Regulation of pre-mRNA splicing: roles in physiology and disease, and therapeutic prospects. Nature Reviews Genetics, (2023); 24(4): 251-269.
23. Ji X, Cao Z, Hao Q, He M, Cang M, et al. Effects of New Mutations in BMPRIB, GDF9, BMP15, LEPR, and B4GALNT2 Genes on Litter Size in Sheep. Veterinary Sciences, (2023); 10(4): 258.
24. Guo X, Wang X, Liang B, Di R, Liu Q, et al. Molecular cloning of the B4GALNT2 gene and its single nucleotide polymorphisms association with litter size in small tail han sheep. Animals, (2018); 8(10): 160.
25. Nenova R, Dimitrova I, Stancheva N, Bozhilova-Sakova M, Tzonev T, Minkova T. Genetic markers associated to improving prolificacy of sheep. A review. Bulgarian Journal of Agricultural Science, (2023); 29(2): 371–377.
26. Getaneh M, Taye M, Alemayehu K, Haile A, Getachew T, Ayalew W. A review on candidate genes associated with sheep fertility traits: Implications for genetic improvement of indigenous sheep breeds in developing countries. Ecological Genetics and Genomics, (2024); 14: 100243.