Polymorphism of the BMPR1B Variants for Prolific Traits in the Indonesian Local Ettawah Goat

TL;DR

This study identifies genetic variations in the BMPR1B gene linked to higher fertility in Indonesian Ettawah goats, suggesting potential for breeding improvements.

Contribution

The study introduces a dual-primer PCR method to detect BMPR1B polymorphisms and links the G allele to increased prolificacy in Ettawah goats.

Findings

The BMPR1B gene's G allele is associated with higher fecundity in Indonesian Local Ettawah Goats.

Homozygous GG genotype showed significantly higher prolificacy compared to AA and GA genotypes.

A dual-primer PCR approach effectively identifies BMPR1B polymorphisms for breeding applications.

Abstract

A polymerase chain reaction (PCR) assay, employing a novel combination of A and G allele primers of the Bovine Morphogenetic Protein Receptor type 1B (BMPR1B) gene, revealed significant polymorphism within the Indonesian Local Ettawah Goat (ILEG). Polymorphism analysis classified the goats into three distinct genotypes: homozygous allele A (AA), heterozygous (AG), and homozygous allele G (GG). These findings also indicated a correlation between the presence of the G allele and increased prolificacy and highlighted the potential of BMPR1B variants as genetic markers for enhanced reproductive performance and overall livestock improvement.

The Indonesian Local Ettawah Goat (ILEG) exhibits substantial genetic variation, suggesting its potential for high productivity and promote sustainable practices in farm animal breeding. This study aimed to investigate the molecular characteristics of prolific ILEG by identifying potential candidate genes through polymerase chain reaction (PCR) analysis of the bone morphogenetic protein receptor type 1B (BMPR1B) gene with two variants: alleles G and A. The research involved PCR amplification and sequencing of the BMPR1B A allele, followed by a combined PCR approach integrating both A and G alleles for genotyping. Blood samples were collected from 73 does with documented prolificacy history and 358 does without prolificacy histories, sourced from seven village breeding operations in East Java. PCR amplification yielded fragments of 556–1181 base pairs in all samples. Haplotype analysis revealed 15 unique haplotypes with a diversity of 0.94 and a mutation frequency of 27.15%. Integration of the BMPR1B alleles G and A revealed polymorphic prolific traits. Polymorphism analysis of 385 ILEGs demonstrated allele frequencies of 0.55 for allele A and 0.45 for the allele G. Average fecundity rates associated with the BMPR1B polymorphism were 1.49 offspring for the homozygous AA, 1.60 for the heterozygous GA, and 1.89 for the homozygous GG. While overall differences among genetic groups were approached statistically significantly (Kruskal–Wallis, p = 0.056), pairwise comparison (Mann–Whitney test) revealed that homozygous GG was significantly associated with higher prolificacy compare to the heterozygous GA (p = 0.029) and homozygous AA (p = 0.040). Similar results were also obtained from data without documented history. These findings suggest that the GG polymorphism of BMPR1B may increase prolificacy in ILEG. Furthermore, the higher frequency of allele G highlights the importance of considering prolificacy traits in breeding selection strategies to enhance sustainable genetic improvement and increase litter size in ILEG. It is recommended to apply dual-primer specific amplification and fragment size differentiation as key molecular approaches for polymorphism of the BMPR1B gene and prolificacy, since these methods can highlight genetic variation and provide valuable markers for breeding programs of the Indonesian Local Etawah Goat.

Figures

Full text

1. Introduction

The increasing impact of climate change and disasters on global food systems reinforces the need to ensure food security from animal protein. Goats, due to their widespread distribution across diverse geographical regions including Asia, Africa, Europe, the Americas, and Australia, and a substantial global population [1,2,3,4,5,6], represent a significant livestock resource. In Indonesia, a Southeast Asian country, the goat population is considerable, with approximately 1.8 million goats spread throughout its 34 provinces. East Java is particularly noteworthy for its high concentration of goats, especially the Indonesian Local Ettawah Goat (ILEG), also known as Peranakan Ettawah, which constitutes the second-largest goat population in rural areas. These goats are primarily raised for dairy and meat using both artificial insemination and natural breeding methods [7].

The Indonesian Local Ettawah Goat (ILEG) is an indigenous Indonesian dual-purpose breed (meat and milk production) exhibiting promising genetic potential for reproductive traits, particularly prolificacy, with the capacity to produce multiple offspring per parturition; however, this inherent potential remains underexploited. Surveys conducted across 19 breeding villages revealed that the litter size of ILEG ranges from 1.42 to 2.12 offspring, averaging 1.72 (<2 offspring per parturition), with a corresponding phenotypic variation of 0.53 and an estimated breeding value (EBV) ranging from 1.48 to 1.74 [8]. These observations corroborate that the on-farm productivity of ILEG has yet to attain its full genetic capacity, indicative of low reproductive efficiency and a significant disparity (productivity gap) between theoretical potential and realized performance. Consequently, the present investigation is highly pertinent for informing practical strategies to enhance the performance of ILEG meat and milk production systems on smallholder farms in East Java, with a particular focus on improving prolificacy.

Prolificacy, commonly measured by litter size, is a key reproductive indicator that directly influences small ruminant productivity and the profitability of farming systems. Defined as the number of offspring per parturition, litter size reflects a goat’s reproductive potential and overall contribution to livestock productivity. This trait exhibits moderate heritability, enabling the transmission of fecundity-related genes across generations and making genetic improvement through marker-assisted selection a promising strategy. Advances in molecular biology have facilitated the identification of candidate genes associated with fertility in small ruminants, including bone morphogenetic protein receptor type 1B (BMPR1B), BMPR15, GDF9, INHA, and INHB [9]. Among these, BMPR1B is of particular interest due to its established role in ovulation mechanisms and litter size, with genetic polymorphisms shown to enhance both ovulation rate and fertility [10,11].

Molecular studies have examined the genetic basis of prolificacy in various goat breeds, revealing diverse responses to BMPR1B variants [12,13,14,15,16,17,18]. For example, positive associations between BMPR1B and prolificacy have been reported in Black Bengal [19,20], Beetal, and Teddy goats [21], whereas inconsistent results have been observed in tropical breeds such as Malabary, Attappady Black [22], and Tibetan Cashmere [23]. For instance, the genetic basis of prolificacy in Malabari and Attappady goats, using polymerase chain reaction (PCR) and sequencing methodologies on limited sample size (n < 50), identified the presence of the heterozygous AG genotype; however, no statistically significant association was observed between this genotype and prolificacy [22]. Similarly, a PCR-RFLP analysis of 216 Tibetan Cashmere goats identified single nucleotide variations (SNVs) in BMPR1B but no polymorphic restriction fragment length patterns. The absence of detectable polymorphism may reflect limited genetic diversity within these populations, breed-specific differences in BMPR1B allele frequencies, and/or methodological differences in polymorphism detection.

Variation in allele distribution, particularly the G and A alleles of BMPR1B, is a key factor influencing these outcomes. These variations likely stem from structural changes in the BMPR1B receptor, which affect the TGF-B signaling pathway involved in folliculogenesis [8]., although this relationship is breed-dependent. Therefore, investigating genetic polymorphism is crucial for identifying associations with prolific traits and understanding the frequency distribution of genes within a population.

Limited molecular research has explored the genetic basis of prolificacy in the Indonesian Local Ettawah Goat (ILEG), with previous studies reporting associations with the GDF9 gene [7] as well as BMP15, BMPR1B, and KISS1 [24,25,26]. For example, two genotypes were identified within exon 1 of the BMPR1B gene using a PCR-RFLP approach, although a definitive causal relationship with prolificacy was not established [25]. Building on these findings, the present study aimed to characterize haplotype diversity and genotype frequency distributions of the BMPR1B gene in ILEG populations in East Java and to evaluate their association with prolificacy. The outcomes are expected to provide useful molecular markers for ILEG breeding programs, supporting improvements in reproductive efficiency and overall productivity.

2. Materials and Methods

2.1. Data Collection

The descriptive study used samples collected from ILEG parents at various village breeding centers of East Java, including Malang I (Ampel Gading), Malang II (Lawang and Singosari), Blitar, Gresik, Lumajang, Nganjuk, Trenggalek, and Banyuwangi. The sample comprised 385 goats without recorded prolificacy data and 75 goats with documented prolificacy records. Prolificacy data were obtained through direct and individual interviews with breeders at their goat housing facilities. Breeders with a herd size of 3–5 animals were selected to enhance the accuracy of their recall regarding individual litter sizes. The recorders employed in the study consisted of 2–4 undergraduate animal husbandry students who received prior training from the researchers. To avoid biased data, an internal cross-validation was implemented. This involved comparing farmers’ reported data with directly observable physiological characteristics of the animals at the time of data collection, including the number of kids suckling, the age of the offspring, and the reproductive status of the does.

2.2. Isolation of DNA

Blood samples were collected from dams using Terumo needles and syringes and immediately transferred into EDTA-coated tubes (Onemed, Surabaya, Indonesia). During transport to the laboratory, samples were stored in an icebox and subsequently either processed immediately or preserved at −20 °C until further use. Genomic DNA was extracted using a modified salting-out method based on [27]. The extraction procedure employed standard laboratory equipment, including a water bath (Yamato (Harumi Triton Square Office Tower Y 36F, 1-8-11 Harumi, Chuo-ku, Tokyo 104-6136, Japan), vortex mixer ((Genie= scientific industries inc. 80 Orville Drive, Suite 102, Bohemia, New York, 11716, USA), centrifuge (Andreas Hettich GmbH Föhrenstraße 12, 78532 Tuttlingen, Germany)),) spectrofotometer UV vis Genesys 10 (Thermo Electron Scientific Instruments LLC, Madison, WI, USA), and refrigerated microcentrifuge Micro 22R (Andreas Hettich GmbH Föhrenstraße 12, 78532 Tuttlingen, Germany)).

The extraction process involved repeated red blood cell (RBC) lysis, cell lysis, organic solvent extraction, DNA precipitation, and resuspension. Briefly, 0.5 mL thawed blood was transferred into a sterile tube and mixed with 1 mL RBC lysis buffer, vortexed, and incubated for 5 min at room temperature, followed by centrifugation at 3000 rpm for 10 min at 4 °C. The supernatant was discarded, and the pellet was subjected to 3–5 additional cycles of RBC lysis until a white pellet was obtained. The pellet was then resuspended in 300 µL cell lysis buffer, vortexed for 5 min, and incubated at 37 °C for 1 h. After cooling to room temperature, 600 µL of phenol:chloroform:isoamyl alcohol (25:24:1, v/v) was added, vortexed, and centrifuged at 10,000 rpm for 10 min. The aqueous phase was transferred to a new sterile tube and re-extracted with chloroform:isoamyl alcohol (24:1, v/v), followed by centrifugation at 10,000 rpm for 10 min.

The resulting supernatant was mixed with 20 µL of 10 M ammonium acetate and 1 mL cold absolute ethanol to precipitate DNA, which was then collected by centrifugation at 10,000 rpm for 10 min at 4 °C. The DNA pellet was washed with 70% ethanol, centrifuged, air-dried at 37 °C, and resuspended in 50 µL TE buffer (pH 7.6). DNA samples were stored at −20 °C for long-term preservation.

DNA quality and concentration were assessed using a UV–Vis spectrophotometer by measuring absorbance at 260 and 280 nm. The A260/A280 ratio ranged from 1.53 to 1.92, and the final DNA yield was 60–150 ng/µL.

2.3. PCR Amplification

Primer was order from Integrated DNA Technology Pte. Ltd.-Singapore. Polymerase chain reaction (PCR) was performed using the following reaction mixture: 2 μL ddH2O, 1 μL forward primer, 1 μL reverse primer, 5 μL Promega GoTaq Green Master Mix (Promega cat number M7122), and 1 μL DNA template. The PCR cycling conditions consisted of an initial denaturation step at 94 °C for 3 min, followed by 35 cycles of denaturation at 94 °C for 30 s, annealing at 52 °C for 30 s, and extension at 72 °C for 30 s. A final extension step was performed at 72 °C for 10 min. Two primer sets were used in this study to amplify the BMPR1B gene: Allele G primer set and Allele A primer set [19], as detailed in Table 1.

Table 1 details the nucleotide sequences of BMPR1B primers targeting the A and G alleles used in this study, along with the corresponding A/G notation for band identification on PCR results. BMPR1B, a serine/threonine protein kinase receptor, belongs to the protein kinase superfamily. STRING database analysis, integrating data from database searches, text mining, experimental evidence, and homology predictions, suggests several functional partners of BMPR1B, with confidence scores ranging from 0.983 to 0.999. The predicted interactions between BMPR1B and these potential partners are illustrated in Figure 1.

2.4. DNA Sequencing

The PCR amplicons were submitted to 1st BASE (PT Genetika Science, Jakarta, Indonesia) for Sanger sequencing. The resulting sequence data was visualized as chromatograms, with each nucleotide base identified by a distinct color: guanine (G) in black, cytosine (C) in blue, adenine (A) in green, and thymine (T) in red.

2.5. Polymorphism

Polymorphism was performed on dams with both recorded and unrecorded prolificacy data. The frequencies of phenotypes, genetic group, and alleles were then analyzed based on the polymorphism results for both groups of dams. To assess differences in the frequencies of the AA, GA, and GG genetic groups.

2.6. Statistical Analysis

PCR data for the BMPR1B alleles (A and G) were classified into three genotypes: GG, GA, and AA. These genotypes were then combined with the mean prolificacy values of each ewe to assess associations. Data normality was evaluated using the Shapiro–Wilk test to determine the suitability of analysis of variance (ANOVA). For data that failed to satisfy normality assumptions, nonparametric tests (Kruskal–Wallis and Mann–Whitney) with Bonferroni correction were applied. Descriptive statistics, including mean, median, standard deviation, and interquartile range (IQR), were used to summarize the data. The IQR, defined as the difference between the 75th percentile (Q3) and the 25th percentile (Q1), was reported as a measure of data dispersion. All statistical analyses were performed in Python version 3.10 using the SciPy library.

3. Results

3.1. PCR Results

PCR analysis was conducted on goat genomic DNA to determine the presence of A and G alleles of the BMPR1B gene. Two primer sets were utilized: one specific to the A allele, generating a 1100 bp amplicon (Figure 2A), and the other specific to the G allele, generating a 100 bp amplicon (Figure 2B). Genetic groups were assigned based on the presence or absence of these amplicons; for instance, samples exhibiting amplification only at 1100 bp were designated as homozygous for A allele (AA) whereas, samples exhibiting amplification only at 100 bp were designated as homozygous for G allele (GG). Samples demonstrating amplification at both 1100 bp and 100 bp were identified as heterozygous (AG).

3.2. Haplotype and Sequencing

Sequencing of the BMPR1B allele A in 19 samples revealed no consistent polymorphisms. Any observed variations were interpreted as individual polymorphisms specific to a single sample. Individual polymorphisms identified in the sequencing results are summarized in Table 2 and Table 3.

Sequence analysis of the BMPR1B amplicons (n = 19) revealed nucleotide frequency variations at specific positions, as summarized in Table 2. Specifically, at positions 133 and 141, adenine (A) was the most frequent nucleotide, observed in 63% and 68% of the sequences, respectively, while cytosine (C) was present in 37% and 32% of sequences. Position 135 showed a high frequency of guanine (G) at 95%, with thymine (T) detected in only 5% of the sequences. Similarly, positions 137 and 143 exhibited a high frequency of guanine (G), at 95% and 63%, respectively, with adenine (A) detected in only 5% and 37%. At position 136, adenine (A) was highly frequent (95%), with guanine (G) detected in the remaining 5%. Finally, position 148 exhibited a thymine (T) frequency of 63%, with guanine (G) observed in 37% of the sequences.

Table 3 presents the results of sequencing analysis aimed at identifying nucleotide variation and corresponding amino acid residue changes within the BMPR1B gene across a 133–150 nucleotide sequence. At nucleotide positions 133 and 135, variation from adenine (A) to cytosine (C) and guanine (G) to thymine (T), respectively, resulted in a change in the methionine to leucine at location 45, due to the codon change from AUG to CUU. Similarly, at positions 136 and 137, variation from adenine (A) to guanine (G) and guanine (G) to adenine (A) led to a change in the serine to aspartic acid at amino acid location 46, resulting from the codon change AGC to GAG. At amino acid position 47, nucleotide changed from A to C at nucleotide sequence 141, with the corresponding codon change from CCA to CCC, did not result in an amino acid change (proline). A nucleotide changed from G to A at nucleotide position 141, leading to the codon change from AGA to AAA at amino acid location 48, resulting in a change from arginine to lysine. Finally, a variation from thymine (T) to guanine (G) was observed at nucleotide position 148 and amino acid location 50, resulting in a change from the amino acid leucine to valine, reflected by the codon change from TTA to GTA.

Haplotype diversity was high (0.9415), with 15 distinct haplotypes identified (Table 4). Haplotypes were identified using PCR sequencing of the BMPR1B gene allele A, which is known to play a crucial role in reproductive traits. This high level of genetic diversity was attributed to geographical differences among the breeding villages. An exception was haplotype 3, which exhibited genetic similarity across the Malang I (animal code AG), Banyuwangi (BW), and Gresik (GR) breeding villages. Additionally, genetic diversity was observed even among different individuals within the same village. For example, individuals from Malang I (animal code AG) displayed four different haplotypes (2, 4, 6, and 12), while individuals from Banyuwangi (BW) showed haplotypes 5, 8, 10, and 13, and individuals from Gresik (GR) had haplotypes 1 and 9. The specific positions of each haplotype are detailed in Table 4.

3.3. BMPR1B Polymorphism of Does with Documented Prolificacy Records

PCR analysis of 73 goats for the BMPR1B gene variants identified three distinct genetic groups: homozygous GG, heterozygous AG, and homozygous AA. Each genetic group presented a specific mean prolificacy of does (Table 5).

Table 5 presents the results of polymerase chain reaction (PCR), performed in duplicate using primer sets specific to both the A and G alleles of the BMPR1B gene. It showed that homozygous A individuals (AA), identified by amplification with only the BMPR1B Allele A primer set, comprising 11 individuals (genotype frequency = 0.15) and exhibited an average prolificacy of 1.49 offspring per kidding per goat.

Table 5 presents the results of polymerase chain reaction (PCR), performed in duplicate using primer sets specific to both the G and A alleles of the BMPR1B gene. The homozygous A allele individuals (AA), identified by amplification with only the BMPR1B allele A primer set, comprising 11 individuals (genotype frequency = 0.15) and exhibited an average prolificacy of 1.49 offspring per kidding per goat.

Moreover, heterozygous individuals (GA), identified by amplification with both BMPR1B Alleles A and G primer sets, comprised 51 individuals (genetic group frequency = 0.70) and exhibited an average prolificacy of 1.60 offspring per kidding per goat. In addition, homozygous variants individuals (GG), identified by amplification with only the BMPR1B allele G primer set, comprised 11 individuals (genotype frequency = 0.15) and exhibited an average prolificacy of 1.89 offspring per kidding per goat.

The results of the descriptive analysis showed that the average prolificacy of GA, GG, and AA individuals is summarized in Table 6. The average number of offspring per kidding per goat for homozygous GG, heterozygous GA, and homozygous AA was 1.89 ± 0.35, 1.60 ± 0.52, and 1.49 ± 0.58, with medians of 2.00, 1.50, and 1.40, and IQRs of 0.00, 1.00, and 0.84, respectively. Meanwhile, ewes with heterozygous GA and homozygous AA had average prolificacy values of 1.60 and 1.49 with medians of 1.50 and 1.40, standard deviations of 0.52 and 0.58, and IQRs of 1.00 and 0.84 offspring per kidding per goat.

The Shapiro–Wilk test indicated that the GG, GA, and AA genotypic groups deviated significantly from normality (p ≤ 0.05). As the data did not satisfy the assumptions required for analysis of variance (ANOVA), non-parametric methods were applied. Overall differences among groups were evaluated using the Kruskal–Wallis test, with pairwise comparisons performed using the Mann–Whitney test. This analytical strategy is appropriate for non-normally distributed data and provides more robust statistical inferences.

Boxplots and dot plots of prolificacy across the GG, GA, and AA genotypic groups (Figure 3) further corroborated these findings. The GG group displayed a notably homogeneous distribution, with most values tightly clustered around 2.0, yielding an interquartile range (IQR) approaching zero. In contrast, the GA and AA groups exhibited greater variability, as evidenced by their higher IQR values.

The results of the median comparison and pairwise tests among the three genotypic groups (GG homozygous, GA heterozygous, and AA homozygous) for prolificacy are summarized in Table 7. Pairwise Mann–Whitney analyses revealed significant differences between GG and GA (p = 0.029) and between GG and AA (p = 0.040), indicating that GG homozygotes generally exhibit higher prolificacy than the other groups. In contrast, no significant difference was observed between GA and AA (p = 0.508), suggesting comparable prolificacy levels in these genotypes.

The GG genotypic group exhibited the highest median prolificacy, with significant increases of 33.3% and 42.9% relative to the GA and AA groups, respectively. By contrast, the difference between GA and AA was minimal (7.1%), indicating comparable prolificacy patterns. Collectively, these findings demonstrate that the GG genotype is consistently associated with higher prolificacy, whereas the GA and AA genotypes exhibit similar reproductive performance.

3.4. Polymorphism Results Without Recorded Prolificacy Data

The distribution of genotype frequencies across seven breeding villages is summarized in Table 8. Polymorphism of 385 goats with unrecorded prolificacy data for the BMPR1B gene revealed three distinct genotypes: homozygous Allele A (AA), heterozygous (AG), and homozygous Allele G (GG). The observed genotype frequencies and corresponding goat numbers were as follows: AA = 0.04 (53 goats), AG = 0.82 (316 goats), and GG = 0.14 (16 goats). Based on these genotype frequencies, the estimated frequency of the A allele was 0.55, while the estimated allele frequency of the G allele was 0.45.

Table 9 presents the distribution of genotypes, including the number of heterozygous, homozygous allele G, and homozygous allele A individuals, and their corresponding genotype frequencies per breeding village across several breeding village locations in East Java Province, Indonesia. Polymerase chain reaction (PCR) analysis revealed that eight villages (Malang I, Malang II, Blitar, Gresik, Lumajang, Nganjuk, Trenggalek, and Banyuwangi) exhibited a high frequency of the heterozygous genotype (GA), which ranged from 73.91% to 94.74%.

The frequency of the homozygous for allele G (GG) ranged from 5.26% to 20.63% in these villages, with the exception of Lumajang. In addition, the homozygous AA exhibited the lowest frequency, ranging from 0% to 6.38%, with the exception of Lumajang. In Lumajang breeding village, the homozygous GG had the lowest frequency (10.87%), while the frequency of the homozygous AA was 15.22%.

Moreover, the results of the pairwise comparison test using Mann–Whitney with Bonferroni correction showed that the comparison between the GA (Heterozygous) and GG (Homozygous G) groups was significantly different, with a p-value of 0.00093 (<0.0167). Similarly, the comparison between GA and AA (Homozygous A) was also significant with a p-value of 0.00086 (<0.0167). In contrast, the comparison between GG and AA did not show a significant difference, with a p-value of 0.038 (p < 0.0167 after Bonferroni correction).

Figure 4 illustrates the distribution of these genetic group frequencies among the different breeding villages. The AA genetic group exhibited the lowest frequency across all breeding villages, with Malang I, Blitar, Gresik, and Nganjuk displaying a frequency of 0. Conversely, the heterozygous genotype displayed the highest frequency across all breeding villages.

4. Discussion

Figure 2 shows that PCR amplification yielded 100 bp fragments for allele G and 1100 bp fragments for allele A, corresponding to BMPR1B gene variants. These results confirm that the primers effectively discriminated between allelic forms, enabling reliable identification of polymorphisms within the population. Samples exhibiting only the 1100 bp band with the allele A primer were classified as AA, those showing only the 100 bp band with the allele G primer as GG, and those displaying both bands as GA.

Table 2 and Table 3 summarize the significant variations in the nucleotide sequences at positions 133, 135, 136, 137, 141, 143, and 148. These variations include transitions of Adenine (A) to Cytosine (C) (A-C), Guanine (G) to Thymine (T) (G-T), Adenine (A) to Guanine (G) (A-G), Guanine (G) to Adenine (A) (G-A), Thymine (T) to Guanine (G) (T-G). Notably, A to C transition (A-C) occurred at positions 133 and 141. G to T transitions (G-T) occurred at positions 135, an Adenine (A) to Guanine (G) transition was observed at position 136, Guanine (G) to Adenine (A) was observed at positions 137 and 143, A to C transition occurred at 141, and a Thymine (T) to Guanine (G) transition (T-G) was observed at position 148. The frequencies of these nucleotide variations varied sharply, from 5% to a high 37%, with an average incidence of 22.57% within the analyzed (ILEG) population. However, variations at positions 133, 141, 143, and 148 exceeded this average, exhibiting 37% and 32% frequencies for the A to C, G to A, T to G and A to C variations, respectively, which highlights the prevalence of these variations and their potential impact. Sequencing analyses revealed that such variation patterns in goats are prevalent and corroborates findings in other small ruminant species [28,29]. Correspondingly, distinct variations are reported [30] at position 773 (G to C) in Assam hill goats utilizing sequencing with the BMPR1B primer. Furthermore, a mutation at position 109 (A to G) was documented in sheep due to crossbreeding between Merino and Garut breeds [10]. A BMPR1B variations at position 159 (A to G) has also been reported in Mehraban sheep [31]. These findings underscore the significance of these nucleotide variations in understanding genetic variations and potential breeding strategies in small ruminants.

Table 2 and Table 3 show three variants in the BMPR1B gene in ILEG, located at 133–148, of 45–50 amino acid locations. At the 45th location, there was a change in the methionine amino acid to leucine (AUG to CUU); at the 46th location, there was a change in the serine amino acid to aspartic Acid (AGC to GAC); at the 47th residue, there was no change in the amino acid proline, but there was a change in the nitrogen base compound (CCA to CCC). At the 48th location, there was a change in amino acid residue from arginine to lysine (AGA to AAA); and at the 50th location, there was a change in amino acid residue from leucine to valine (TTA to GTA).

These amino acid variations resulted in different amino acids, yet failed to affect the classification of the amino acid group (essential or non-essential) except at 48th amino acid location there was change from non-essential amino acid to essential amino acid (arginine to lysine). Specifically, variations at location 45 and 50 also resulted in amino acids that retained the properties of the original essential group. The protein resulting from variations at location 46 produced a non-essential group, while location 48 produced an essential group. Reportedly [32], the essential amino acids were methionine, lysine, threonine, leucine, isoleucine, valine, phenylalanine, tryptophan, cysteine, arginine, and histidine, while non-essential amino acids were alanine, aspartate, asparagine, glutamine, glutamate, cysteine, tyrosine, serine, proline, and glycine.

Table 4 show the haplotype diversity of ILEG in the East Java region of Indonesia had a high category of 0.9415, and 94.5% haplotype was recorded [33] between lineage A and B based on analysis of the mitochondrial D-loop region in goats in the Yangtze River region, China. Moreover, high haplotype diversity in seven goat breeds in Russia was similarly reported [34] at 0.843–1; in addition, a study [35] claims that haplotype diversity in seven goat breeds in China ranged from 0.943–1000 with a high category, which corroborated the finding of high genetic diversity in East Asia goats [36].

Table 5, Table 6 and Table 7 and Figure 3 highlight the genetic group frequencies within the population, indicating that AA, AG, and GG were present at 0.15, 0.70, and 0.15, respectively. The BMPR1B gene is connected to two distinct alleles: Allele A, which does not correlate with high fecundity, and the potent Allele G, regarded as a prolific carrier gene. Our PCR analysis, utilizing primers specific to both alleles, confirmed the existence of AA, AG, and GG genotypes in the sample.

Based on statistic descriptive and analysis statistics of prolificacy results AA genetic group yielded an average of 1.49 ± 0.48 offspring; GA produced 1.60 ± 0.52 whereas GG averaged 1.89 ± 0.58, median 33.30–42.90% higher than GA and AA. Significant differences were observed between GG and GA, as well as between GG and AA, whereas no difference was detected between GA and AA. This was supported by the IQR of 0.00 in the GG group, indicating a homogeneous distribution of prolificacy, in contrast to the higher IQR values (0.84–1.00) observed in GA and AA, which reflected greater variability. These findings suggest that individuals with the GG genotype exhibit consistently high reproductive performance, a desirable trait in productivity-related contexts. By comparison, the phenotypic diversity observed in GA and AA may reflect both genetic heterogeneity and environmental influences. For breeding programs, the stability of the GG genotype underscores its potential value in improving prolificacy, while the variability in GA and AA highlights the need for targeted management strategies to optimize their productivity. Collectively, these results indicate that the GG genotype is strongly associated with enhanced prolificacy, whereas GA and AA genotypes display comparable reproductive performance.

This evidence is supported by findings that Black Bengal goats had average litter sizes of 2.7, 3.04, and 3.11 for AA, AG, and GG genotypes, respectively [20]. The prolificacy of ILEG is lower than that of Black Bengal goats, which may be attributed to genetic background differences. Specifically, the frequency of the BMPR1B variant in ILEG (27.15%) is lower than that reported for Black Bengal goats (40%) [20]). Although different genotyping methods were applied (dual-primer PCR in the present study and PCR-RFLP in Black Bengal goats) both approaches targeted the same BMPR1B SNP associated with prolificacy. Thus, the genotype frequencies obtained here are suitable for direct comparison and highlight genetic differences between the two populations. Additionally, research on local Bangladeshi sheep that employed the same genomic primers reported averages of 1.19 ± 0.07, 1.44 ± 0.09, and 2.13 ± 0.09 for the ++, +B, and BB genotypes, respectively [37]. Notably, numerous studies consistently reveal that heterozygous genotypes often outperform AA homozygous genotypes in terms of prolificacy [38].

Variations in reproductive performance among ILEG and Black Bengal goats with GG, GA, and AA genotypes reflected differences in prolificacy, largely attributable to breed-specific genetics and the polymorphism detection methods employed. Genetics play a critical role in maximizing productivity in these valuable livestock populations [39,40]. The present findings in ILEG are consistent with those reported for Black Bengal, Beetal, and Teddy goats, where BMPR1B variants show a positive association with prolificacy. By contrast, studies in Malabary, Attappady [22], and Tibetan Cashmere [23] goats have not demonstrated a clear association, likely due to limited genetic variation, small sample sizes, methodological differences, or environmental influences, which were beyond the scope of the present study. These comparisons underscore the breed-specific influence of BMPR1B on prolificacy and highlight the value of molecular screening in optimizing selection strategies.

Genotyping of 385 does of ILEG (Table 8) revealed critical insights into the BMPR1B gene within the breeding villages. The allele frequencies for A and G were 0.55 and 0.45, respectively. Moreover, the genetic group frequencies were 0.04 for the AA, 0.82 for the GA, and 0.14 for the GG. These findings underline the very low frequency of homozygous AA, emphasizing that high productivity traits are desirable and important criteria for selecting ILEG in rural breeding farms. The very low frequency of AA also correlated negatively with productivity, especially the number of kids per litter, meaning that the higher the livestock productivity, the lower the frequency of homozygous AA. Furthermore, this aligns with the research by Zhao et al. [41], which supports the notion that positive selection for highly prolific traits in Guizhou Black goats and Meigu goats has significantly enhanced the proliferation of high-performing individuals while simultaneously reducing the prevalence of lower-performing traits. This evidence strongly advocates for the strategic focus on prolificacy to improve breeding outcomes and overall livestock performance.

In village breeding populations, heterozygous (GA) individuals represented the largest proportion (Table 7; Figure 3). The frequency of the GA genotype was significantly higher than that of either homozygous group (GG or AA). On average, the AA genotype was the least common, occurring at a frequency 231% lower than GG, although this difference was not statistically significant. This distribution pattern suggests effective positive selection for prolificacy traits within the village breeding goat populations. Despite the lack of formal recording practices, farmers relied on observable traits—such as appearance, size, health, and fertility—and their own experiences to select parent candidates for ILEG. This intuitive selection process aligned perfectly with the higher prevalence of heterozygous GA phenotypes, underscoring that ILEG was crucial in promoting long-term genetic diversity within livestock populations. Local customs, such as the exchange of bucks among breeders, supported this diversity of ILEG. These findings support the idea that the advantages conferred by heterozygous GA arise from adaptive mutations, highlighting the potential for sustainable breeding practices that capitalize on genetic diversity [42,43].

5. Conclusions

The study of the BMPR1B gene in Indonesian Local Ettawah Goat (ILEG) revealed significant insights into the genetic basis of prolificacy, with potential applications for enhancing livestock breeding programs. The presence of the allele G correlated with increased fecundity, highlighting BMPR1B’s utility as a genetic marker for reproductive performance. Polymorphism analysis identified three distinct genetic groups (AA, AG, and GG) and their respective frequencies in ILEG populations, both with and without recorded prolificacy data. The heterozygous genotype (GA) was the most prevalent in village breeding goat populations, suggesting successful, positive selection for prolific traits, even in the absence of formal recording practices. The homozygous GG genotype was associated with higher fecundity, and the elevated frequency of allele G underscores the relevance of incorporating prolificacy traits into breeding selection strategies. Accordingly, the breeding program prioritized GG does, while AG bucks were retained to preserve genetic diversity within the herd. Highly prolific individuals were identified through PCR-based polymorphism analysis, enabling the establishment of a core breeding group. Based on these findings, the implementation of sustainable breeding strategies is recommended for ILEG using dual-primer PCR genotyping to efficiently screen for GG and AG individuals, which could increase average prolificacy by an estimated 15–20%.