Genetic variants of BMP15 and GDF9 genes and their association with litter size in Tibetan sheep

1 Introduction

The Qinghai–Tibet Plateau is the highest plateau in the world, with an average altitude of more than 4,000 m, and is often referred to as the “Roof of the World” and the “Third Pole of the Earth.” This region is characterized by harsh environmental conditions, including hypoxia, low temperature, and strong ultraviolet radiation (Sun et al., 2022; He et al., 2023). Tibetan sheep (Ovis aries) have long inhabited the Qinghai–Tibet Plateau and represent an important indigenous sheep genetic resource adapted to cold and hypoxic high-altitude environments. Through long-term natural selection, Tibetan sheep have evolved into a unique population with stable adaptive characteristics to high-altitude hypoxia (Liu et al., 2020; Lu et al., 2022). They also play important roles in supporting the livelihoods of local herders, livestock production, and the stability of grassland ecosystems on the Qinghai–Tibet Plateau (Jing et al., 2020). However, Tibetan sheep are seasonal breeders and generally produce single lambs, with an average litter size of approximately 1.05 lambs per lambing, which limits the production efficiency of local animal husbandry to some extent (Zhang, 2019; Yang et al., 2022). Previous studies have shown that high-altitude hypobaric hypoxia and oxidative stress can affect luteal development and function in sheep, thereby potentially reducing reproductive efficiency (Parraguez et al., 2013). Furthermore, ovarian transcriptomic studies in Tibetan sheep have shown that high-altitude stress can suppress the expression of follicular development-related genes and interfere with LH/FSH signaling pathways, suggesting that plateau environmental stress may regulate reproductive performance in Tibetan sheep by affecting ovarian function and follicular development (Li et al., 2022b). Therefore, the relatively low litter size of Tibetan sheep may not be determined by a single factor, but may be associated with multiple factors formed during long-term adaptation to the plateau environment, including hypoxia adaptation, energy metabolism, endocrine regulation, and changes in ovarian function.

BMP15 and GDF9 are oocyte-secreted members of the transforming growth factor-β (TGF-β) superfamily and are recognized as key intraovarian regulators of female fertility in mammals (Dong et al., 1996; Nilsson and Skinner, 2002; Belli and Shimasaki, 2018; Liu et al., 2021; Juengel et al., 2004). During folliculogenesis, these two factors are mainly expressed in oocytes of growing follicles and act in a paracrine manner on surrounding granulosa and cumulus cells (Juengel et al., 2004; Knight and Glister, 2006). BMP15 regulates granulosa cell proliferation, cumulus expansion, and follicular development, thereby supporting normal follicular growth and oocyte maturation (Fabre et al., 2006; Knight and Glister, 2006). Similarly, GDF9, which is located on ovine chromosome 5, is predominantly expressed in oocytes and ovaries and is essential for early follicular development, cumulus–oocyte complex function, normal ovulation, fertilization, and subsequent reproductive performance (Dong et al., 1996; Sadighi et al., 2002; Hanrahan et al., 2004; Tang et al., 2018).

Through their effects on granulosa cell activity, follicular development, oocyte maturation, and ovulation, BMP15 and GDF9 are closely associated with ovulation rate and litter size in sheep. Previous studies have reported that sequence variations in these genes can influence prolificacy, and BMP15 and GDF9 have been investigated as candidate genes for prolificacy in Chinese sheep breeds such as Small-tailed Han sheep, as well as in other sheep populations (Chu et al., 2005; Najafabadi et al., 2021; Ji et al., 2023). In addition, BMP/Smad pathway-related molecules have been investigated in Hu sheep with different fecundity levels (Xu et al., 2010). However, the genetic characteristics of BMP15 and GDF9 and their potential associations with litter size in Tibetan sheep remain unclear. Therefore, the present study aimed to characterize the molecular features, expression patterns, and single nucleotide polymorphisms of BMP15 and GDF9 and to evaluate their associations with litter size in Tibetan sheep.

2 Materials and methods

2.1 Animals and sample collection

All experimental procedures were conducted in accordance with the guidelines for animal experiments issued by Qinghai University and were approved by the Animal Research Ethics Committee of Qinghai University (No. PJ202401-75). A total of 200 healthy adult Tibetan ewes were selected for genomic DNA extraction and association analysis between the selected SNP loci and litter size. All animals were obtained from the same commercial farm in Gonghe County, Hainan Tibetan Autonomous Prefecture, Qinghai Province, China, were raised in the same flock, and were maintained under a unified production-management system throughout the study. The experimental ewes were 3–4 years old, weighed 40–55 kg, and had complete and consecutive litter-size records for 2–3 parities. All animals were naturally mated during the same natural breeding season, namely autumn, and the lambing records were registered and verified by farm technicians. In the present study, litter size was defined as the number of lambs born per ewe per lambing. The animals were managed under a grazing system with supplementary feeding, and the same nutritional-management program was applied during the mating and gestation periods. Housing conditions, grazing schedules, feeding standards, and health-management practices were kept consistent among all animals. The inclusion criteria were as follows: clinically healthy adult Tibetan ewes with a normal reproductive history, complete reproductive records for 2–3 consecutive parities, and complete and reliable litter-size records. Ewes with a history of abortion, reproductive-tract disease, severe malnutrition, serious systemic disease, or incomplete reproductive records were excluded. Body condition scores were not systematically recorded at the beginning of the study and could not be reliably reconstructed retrospectively; therefore, body condition score was not included as a covariate in the statistical analysis. Nevertheless, ewes showing severe malnutrition or visibly abnormal body condition were excluded during clinical screening. No exogenous reproductive hormones, including estrus-synchronization or ovulation-induction treatments, were administered during the mating and litter-recording periods.

Blood samples were collected from the jugular vein of each Tibetan ewe; 5 mL of blood per individual was transferred to tubes containing 1 mg/mL ethylene diamine tetraacetic acid (EDTA) as an anticoagulant. All samples were transported to the laboratory and stored at −20 °C.

Six 6-month-old, clinically healthy Tibetan ewes from the same flock were selected for exploratory tissue-expression profiling. The animals had similar body weights (35.32 ± 1.96 kg) and were maintained under the same feeding and management conditions. Gross morphological examination at the time of sampling revealed no apparent abnormalities in their reproductive tracts. Six biological replicates were included to obtain preliminary information on the tissue-expression distribution of BMP15 and GDF9, rather than to compare expression among different reproductive stages. For quantitative real-time PCR analysis, each biological sample was analyzed in triplicate as technical replicates. Immediately after slaughter, the ovaries were examined macroscopically, and developing follicles and corpora lutea were visible on the ovarian surfaces of all ewes. However, these morphological features were used only to describe ovarian status and were not used to assign the animals to specific estrous-cycle stages. Because serum reproductive hormone concentrations were not measured and detailed estrous-cycle records were unavailable, the precise estrous-cycle stage of each ewe could not be determined. Therefore, the results should be interpreted as exploratory tissue-expression profiles obtained under consistent feeding and management conditions rather than as expression data generated under strictly controlled estrous-cycle stages. Tissue samples, including the hypothalamus, heart, lung, kidney, liver, spleen, muscle, uterus, ovary, oviduct, rumen, and duodenum, were rapidly collected. All tissue samples were immediately snap-frozen in liquid nitrogen, transported to the laboratory, and stored at −80 °C until total RNA extraction.

2.2 Extraction of blood DNA and quality inspection

Total DNA extraction from the blood of the test Tibetan sheep was performed using the TIANamp Genomic DNA Kit. After testing the quality, DNA was uniformly diluted to 10 ng/μL until further application.

2.3 RNA extraction and cDNA synthesis

Total RNA was extracted from the hypothalamus, heart, lung, kidney, liver, spleen, muscle, uterus, ovary, oviduct, rumen, and duodenum of Tibetan sheep using the TRIzol method. RNA integrity was assessed by 1.5% agarose gel electrophoresis, and RNA concentration was measured using a nucleic acid/protein analyzer. The extracted RNA was stored at −20 °C until further use. Total RNA (1 μg) was reverse-transcribed into cDNA using FastKing gDNA Dispelling RT SuperMix KR118 (TIANGEN Biotechnology Co., Ltd., Beijing, China). Each cDNA sample was then diluted four-fold with RNase-free ddH₂O and stored at −20 °C until further use.

2.4 Gene cloning

The primers for BMP15 and GDF9 were designed based on the ovine BMP15 and GDF9 sequences using Primer Premier 5.0 software (Premier Biosoft International, Palo Alto, CA, USA) (Table 1). These primers were used to amplify cDNA fragments of Tibetan sheep BMP15 and GDF9. PCR amplification was performed using cDNA as the template, and the target PCR products were separated by 1% agarose gel electrophoresis. The PCR products were purified using a Gel Extraction Kit (TIANGEN Biotechnology Co., Ltd., Beijing, China) according to the manufacturer’s instructions. The purified fragments were then ligated into the pMD19-T vector (Takara, Japan) and transformed into competent Escherichia coli cells. Positive recombinant clones were identified by PCR and subsequently sent to Shanghai Genesky Biotechnologies Inc. (Shanghai, China) for sequencing.

Table 1. The primer sequence information used for the cDNA cloning and expression analysis of BMP15 and GDF9.

2.5 Characteristics and bioinformatics analyses of BMP15 and GDF9 cDNA

The bioinformatic characteristics of the target protein were systematically analyzed using a series of professional online bioinformatic analysis tools: Open reading frame analysis was conducted via the ORF finder tool on the NCBI platform (https://www.ncbi.nlm.nih.gov/orffinder); the physical and chemical properties of the protein were calculated using ProtParam (https://web.expasy.org/protparam/); ProtScale (https://web.expasy.org/protscale/) was employed to predict the hydrophilic and hydrophobic properties of the protein sequence. Signal peptide prediction was performed with the SignalP 5.0 Server (http://www.cbs.dtu.dk/services/SignalP/), and the TMHMM 2.0 Server (http://www.cbs.dtu.dk/services/TMHMM/) was used to identify transmembrane structural domains of the protein. NetPhos 3.1 Server (http://www.cbs.dtu.dk/services/NetPhos/) was applied for the prediction of protein phosphorylation sites; SOPMA (https://npsa.lyon.inserm.fr/) was utilized to analyze the secondary structure of the protein, and the tertiary structure of the protein was predicted through the SWISS-MODEL online server (http://swissmodel.expasy.org/).

2.6 Real-time quantitative PCR

For quantitative real-time PCR (qPCR) analysis of BMP15 and GDF9, specific primers were designed using Primer Premier 5.0 software (Premier Biosoft International, Palo Alto, CA, USA). The primer design was based on the ovine BMP15, GDF9, and GAPDH sequences downloaded from the National Center for Biotechnology Information (NCBI) database. The specificity of the qPCR primers was further rechecked by sequence alignment against the ovine reference genome to reduce the possibility of nonspecific amplification. Primer sequences are listed in Table 1, and all oligonucleotides were synthesized by Sangon Biotech (Shanghai) Co., Ltd.

Using GAPDH as the internal reference gene, the mRNA expression levels of BMP15 and GDF9 in different Tibetan sheep tissues were quantified using an ArchimedX real-time PCR instrument. The total qPCR reaction volume was 20 μL, consisting of 10 μL 2× SuperReal Color PreMix, 0.6 μL forward primer (10 μM), 0.6 μL reverse primer (10 μM), 2 μL cDNA template, 6.3 μL RNase-free ddH₂O, and 0.5 μL Dilution Buffer. The thermal cycling conditions were as follows: initial denaturation at 95 °C for 15 min, followed by 40 cycles of denaturation at 95 °C for 10 s, annealing at 60 °C for 32 s, and extension at 72 °C for 20 s. No-template controls were included in each qPCR run to exclude contamination. For each tissue, samples from six ewes were used as biological replicates, and each sample was analyzed in triplicate as technical replicates. Only reactions showing stable amplification among technical replicates were used for subsequent expression analysis. The relative expression levels of BMP15 and GDF9 were calculated using the 2−ΔΔCt method.

2.7 SNP selection and genotyping

The four SNPs analyzed in this study were selected from whole-genome resequencing data previously generated for Tibetan sheep by our research group. Functional annotation of the four SNPs was performed using Ensembl Variant Effect Predictor (VEP, release 116) based on the ovine (Ovis aries) ARS-UI_Ramb_v3.0 reference genome assembly (RefSeq accession: GCF_016772045.2). The genomic coordinates of BMP15 and GDF9 were referenced to the X-chromosome sequence NC_056080.1 and chromosome 5 sequence NC_056058.1, respectively, whereas transcript-level annotation was based on the RefSeq transcripts NM_001114767.2 and NM_001142888.2. The exon number, HGVSc, codon change, amino acid position and substitution, HGVSp, dbSNP identifier, and predicted variant consequence were determined from the VEP output. The potential functional effects of missense variants were further evaluated using SIFT and BLOSUM62. For BMP15 g.54284258G>A and GDF9 g.42115010C>T, the allele order reported in the original genotyping report was opposite to the reference allele–alternative allele orientation on the forward strand of the reference genome. Therefore, to maintain consistency with the original genotyping dataset, the original locus labels were retained in the subsequent analyses. BMP15 and GDF9 encode oocyte-derived growth factors involved in follicular development, ovulation, and the regulation of litter size. Accordingly, exonic variants, particularly missense variants, may affect protein physicochemical properties, spatial conformation, receptor-binding capacity, or downstream biological functions. Although synonymous variants do not alter the amino acid sequence, they may still influence gene function through codon-usage bias, mRNA stability, or translational efficiency. On this basis, the four SNPs were selected for population genotyping and association analysis with litter size. Detailed information on their genomic positions, variant types, and amino acid substitutions is provided in the Results section and Table 2.

Table 2. Functional annotation of SNPs in the BMP15 and GDF9 genes.

The variants BMP15 g.54284258G>A and GDF9 g.42115010C>T designated in the original genotyping data correspond to the standardized variants on the forward reference strand NC_056080.1:g.54284258A>G and NC_056058.1:g.42115010T>C, respectively.

Genomic DNA was extracted from jugular venous blood samples. Genotyping of the four selected SNP loci within the BMP15 and GDF9 genes was performed using improved multiplex ligase detection reaction (iMLDR) following PCR amplification on an ABI 3730XL Genetic Analyzer (Applied Biosystems, USA), with technical support provided by Shanghai Genesky Biotechnologies Inc. The primer sequences used for PCR amplification and iMLDR genotyping are summarized in Table 3.

Table 3. iMLDR genotyping primer.

Degenerate bases Y, R, K and M are used in the primers, corresponding to C/T, G/A, G/T and A/C respectively.

PCR reactions were conducted in a total volume of 20 μL, containing 1× GC-I buffer (Takara), 3.0 mM Mg²⁺, 0.3 mM dNTPs, 1 U HotStarTaq polymerase (Qiagen Inc., Germantown, MD, USA), 1 μL of genomic DNA template, and 1 μL each of forward and reverse primers. The amplification program was set as follows: initial denaturation at 95 °C for 2 min; 11 cycles of denaturation at 94 °C for 20 s, annealing at 65 °C for 40 s, and extension at 72 °C for 1.5 min; followed by 24 cycles of denaturation at 94 °C for 20 s, annealing at 59 °C for 30 s, and extension at 72 °C for 1.5 min; with a final extension step at 72 °C for 2 min. The PCR products were purified via incubation with 5 U of shrimp alkaline phosphatase and 2 U of exonuclease I at 37 °C for 1 h, followed by heat inactivation at 75 °C for 15 min to remove residual dNTPs and primers.

For the iMLDR assay, two allele-specific probes and one fluorescently labeled probe were used for each SNP locus (Table 3). The ligation reaction mixture contained 1 μL of 10× ligation buffer, 0.25 μL of thermostable Taq DNA ligase, 0.4 μL of 1 mM 5′-end ligation primer mix, 0.4 μL of 2 mM 3′-end ligation primer mix, 5 μL of purified multiplex PCR product, and nuclease-free water. The iMLDR cycling conditions consisted of 38 cycles of 94 °C for 1 min and 56 °C for 4 min, and the reaction products were stored at 4 °C after completion. Subsequently, 0.5 μL of the iMLDR reaction product was mixed with 9 μL of Hi-Di formamide and 0.5 μL of LIZ-600 size standard, denatured at 95 °C for 5 min, and separated on an ABI 3730XL Genetic Analyzer. Genotyping data were processed and analyzed using GeneMapper Software v4.1 (Applied Biosystems, USA).

2.8 Statistical analysis

Genotype frequency, gene frequency, homozygosity (Ho), heterozygosity (He), effective allele numbers (Ne) and polymorphic information content (PIC) were calculated for each genotype of the BMP15 and GDF9 genes in Tibetan sheep using Microsoft Excel 2010 software. heterozygosity (He), Effective allele numbers (Ne) and Polymorphic information Content (PIC) were calculated with reference to the method of Li et al (Li et al., 2022a).

The relative expression of BMP15 and GDF9 genes was analyzed using the 2^-ΔΔCT method, where ΔCt=Ct_{(target gene)}- Ct₍_GAPDH₎, and 2^-ΔΔCT represents the fold change in the relative expression of each target gene normalized to GAPDH. Before statistical testing, G*Power 3.1.9.7 software was used to perform post-hoc statistical power analysis to verify the rationality of the sample size and ensure the reliability of the subsequent experimental statistical results. One-way analysis of variance (ANOVA) was performed using SPSS 23.0 software, and multiple comparisons were conducted using Duncan’s multiple range test. A P-value less than 0.05 was considered statistically significant. All experimental results are presented as mean ± standard deviation, and all graphs were plotted using GraphPad Prism 8.0 software.

General Linear Model was used to analyze the correlation between polymorphisms of each gene and litter size, General Linear Model procedure model was as follows:

Yijn= μ + Pi+ Gj+ IPG+ eijn

Where Y_ijn is the phenotypic observation value; μ is the overall population mean; P_i is the fixed effect of the i^th parity (i = 1, 2, or 3); G_j is the effect of the j^th genotype (j = 1, 2, or 3), I_PG is the interactive effect of parity and genotype; and e_ijn represents random error.

3 Results

3.1 Molecular cloning and sequence analysis of Tibetan sheep BMP15 and GDF9 genes

3.1.1 Sequence and open reading frame analysis of cloned BMP15 and GDF9 fragments

Open reading frame prediction was performed on the cloned Tibetan sheep BMP15 and GDF9 cDNA fragments using the NCBI ORF Finder tool, followed by alignment with ovine reference sequences from the NCBI database. The analysis identified predicted ORFs of 468 bp and 405 bp within the cloned BMP15 and GDF9 fragments, encoding 155 and 134 amino acids, respectively. Considering that these predicted amino acid lengths are markedly shorter than those of the reported full-length ovine BMP15 and GDF9 precursor proteins, and that alignment with the NCBI reference sequences confirmed incomplete coverage of the full-length CDSs, the sequences obtained in this study should be regarded as partial coding fragments of BMP15 and GDF9 rather than complete full-length CDSs (Juengel et al., 2004). Therefore, subsequent bioinformatic analyses were performed based on the predicted amino acid sequences encoded by these cloned fragments.

3.1.2 Multiple sequences alignment and phylogenetic analysis

Multiple sequence alignment and homology analysis were performed using the predicted amino acid sequences encoded by the cloned BMP15 and GDF9 fragments from Tibetan sheep and homologous sequences from Ovis aries, Bos taurus, Bos mutus, Homo sapiens, Sus scrofa, Mus musculus, Maylandia zebra, Canis lupus familiaris, Pan troglodytes, Macaca mulatta, and Gallus gallus. The alignment results showed that the predicted amino acid sequences of the cloned Tibetan sheep BMP15 and GDF9 fragments had the highest homology with those of domestic sheep (Ovis aries), followed by yak (Bos mutus). The lowest homology was observed between the Tibetan sheep BMP15 fragment and Canis lupus familiaris, and between the Tibetan sheep GDF9 fragment and Maylandia zebra, respectively (Figure 1).

3.1.3 Prediction of protein physicochemical properties

The physicochemical properties of the protein fragments encoded by the cloned BMP15 and GDF9 sequences in Tibetan sheep were predicted using the ExPASy ProtParam online tool. The results indicated that both predicted protein fragments consisted of five elements: carbon (C), hydrogen (H), nitrogen (N), oxygen (O), and sulfur (S). For the protein fragment encoded by the cloned BMP15 sequence, the molecular formula was C₇₇₉H₁₂₈₁N₂₁₇O₂₀₉S₇, with a total atom number of 2493, an extinction coefficient of 11460 M^-1·cm^-1, an instability index of 35.21, and an aliphatic index of 108.13. For the protein fragment encoded by the cloned GDF9 sequence, the molecular formula was C₆₈₆H₁₀₆₂N₁₉₄O₁₈₅S₆, with a total atom number of 2133, an extinction coefficient of 24200 M^-1·cm^-1, an instability index of 53.62, and an aliphatic index of 82.39. Because an instability index below 40 indicates predicted stability, the BMP15-encoded protein fragment was predicted to be stable, whereas the GDF9-encoded protein fragment was predicted to be unstable (Table 4).

Table 4. Amino acid composition of the protein fragments encoded by the cloned BMP15 and GDF9 sequences in Tibetan sheep.

3.1.4 Protein hydrophilicity prediction

Hydropathy profiles of the protein fragments encoded by the cloned BMP15 and GDF9 sequences in Tibetan sheep were analyzed using the ExPASy ProtScale tool with the Kyte–Doolittle method, and their grand average of hydropathicity (GRAVY) values were calculated using ExPASy ProtParam. As shown in Figure 2, both protein fragments contained alternating hydrophobic and hydrophilic regions, with relatively hydrophobic regions observed near their N-termini. The GRAVY values of the BMP15- and GDF9-encoded protein fragments were 0.026 and −0.323, respectively. Therefore, the BMP15-encoded fragment was nearly neutral in overall hydropathicity, with only slight hydrophobicity, whereas the GDF9-encoded fragment was hydrophilic overall.

3.1.5 Protein signal peptide prediction

Signal peptide prediction was performed for the protein fragments encoded by the cloned BMP15 and GDF9 sequences in Tibetan sheep using the SignalP 5.0 online server with the Eukarya model. As shown in Figure 3, no typical signal peptide was predicted in the BMP15-encoded protein fragment, with predicted probabilities of 0.0812 for a Sec/SPI signal peptide and 0.9188 for other proteins. In contrast, a Sec/SPI-type signal peptide was predicted in the GDF9-encoded protein fragment, with a predicted probability of 0.8842. The potential cleavage site was located between amino acid residues 25 and 26 (LDS|LP), with a cleavage-site probability of 0.3352. These results suggest that the cloned GDF9 fragment retains an N-terminal sequence feature potentially associated with targeting to the classical secretory pathway.

3.1.6 Transmembrane structure prediction

Transmembrane domains in the protein fragments encoded by the cloned BMP15 and GDF9 sequences were predicted using TMHMM Server v.2.0. As shown in Figure 4, one N-terminal transmembrane-like segment was predicted in the BMP15-encoded protein fragment, spanning approximately amino acid residues 10–27, with a predicted topology of o10–27i. In contrast, no transmembrane helix was predicted in the GDF9-encoded protein fragment.

3.1.7 Prediction of protein N-glycosylation sites and phosphorylation sites

Potential N-linked glycosylation sites in the protein fragments encoded by the cloned BMP15 and GDF9 sequences in Tibetan sheep were predicted using the NetNGlyc 1.0 Server. As shown in Figure 5, one potential N-linked glycosylation site was predicted at Asn95 in the BMP15-encoded protein fragment. This site was located within the consensus N-X-S/T sequence motif NRT, and its prediction score exceeded the threshold of 0.5. In contrast, no potential N-linked glycosylation site with a prediction score above the threshold was identified in the GDF9-encoded protein fragment.

Potential phosphorylation sites in the protein fragments encoded by the cloned BMP15 and GDF9 sequences in Tibetan sheep were predicted using the NetPhos 3.1 Server, with a prediction threshold of 0.5. As shown in Figure 6, 13 potential phosphorylation sites were identified in the BMP15-encoded protein fragment, comprising nine serine residues, two threonine residues, and two tyrosine residues. Similarly, 13 potential phosphorylation sites were identified in the GDF9-encoded protein fragment, comprising six serine residues, six threonine residues, and one tyrosine residue.

3.1.8 Protein secondary structure prediction

The secondary structures of the protein fragments encoded by the cloned BMP15 and GDF9 sequences in Tibetan sheep were predicted using the SOPMA online tool. As shown in Figure 7, the predicted secondary structure of the BMP15-encoded protein fragment consisted of 50.97% random coils, 32.26% α-helices, 9.03% extended strands, and 7.74% β-turns. The predicted secondary structure of the GDF9-encoded protein fragment consisted of 52.24% random coils, 32.84% α-helices, 11.94% extended strands, and 2.99% β-turns. Thus, the proportions of the four structural elements in both protein fragments followed the order of random coils, α-helices, extended strands, and β-turns from highest to lowest. Random coils and α-helices were the predominant predicted secondary-structure components of both protein fragments. The preliminary tertiary structure models of protein fragments encoded by cloned BMP15 and GDF9 sequences were further constructed via the SWISS-MODEL online server, and the predicted structures are presented in Figure 8. Given that only partial coding fragments of the two genes were cloned in this study, these tertiary structure models are merely preliminary predictive results.

3.2 Tissue expression analysis of BMP15 and GDF9 genes in Tibetan sheep

BMP15 and GDF9 genes were expressed in 12 tissues including hypothalamus, heart, liver, spleen, lung, kidney, ovary, uterus, oviduct, rumen, duodenum and muscle of Tibetan sheep (Figure 9). The BMP15 gene was expressed highest in oviduct (P<0.05), followed by liver (P<0.05), and lowest in muscle. The GDF9 gene was preferentially expressed in ovary (P<0.05), followed by kidney and uterus (P<0.05), and lowest in heart.

3.3 Analysis of polymorphisms in BMP15 and GDF9 genes and their correlation with litter size in Tibetan sheep

The VEP results classified the four SNPs as two synonymous and two missense variants (Table 2). BMP15 NC_056080.1:g.54284258A>G (rs55628000; originally designated as g.54284258G>A in the genotyping report), located in exon 2, corresponded to NM_001114767.2:c.752T>C. The codon changed from CTG to CCG, resulting in a leucine-to-proline substitution at amino acid position 251 (NP_001108239.1:p.Leu251Pro); this SNP was therefore classified as a missense variant. BMP15 NC_056080.1:g.54283966C>T (rs587958054), also located in exon 2, corresponded to NM_001114767.2:c.1044G>A. The codon changed from GTG to GTA without altering the encoded valine residue (NP_001108239.1:p.Val348=); this SNP was therefore classified as a synonymous variant. GDF9 NC_056058.1:g.42116345C>T (rs410123449), located in exon 1, corresponded to NM_001142888.2:c.260G>A. The codon changed from CGC to CAC, resulting in an arginine-to-histidine substitution at amino acid position 87 (NP_001136360.2:p.Arg87His); this SNP was therefore classified as a missense variant. GDF9 NC_056058.1:g.42115010T>C (rs160076413; originally designated as g.42115010C>T in the genotyping report), located in exon 2, corresponded to NM_001142888.2:c.477A>G. The codon changed from CTA to CTG without altering the encoded leucine residue (NP_001136360.2:p.Leu159=); this SNP was therefore classified as a synonymous variant. SIFT predicted both BMP15 p.Leu251Pro and GDF9 p.Arg87His to be tolerated, with scores of 0.16 and 0.12, respectively, while their corresponding BLOSUM62 scores were −3 and 0. These in silico predictions provide only preliminary evidence of the potential functional effects of the two missense variants, and their biological consequences require further experimental validation.

No significant deviations from Hardy–Weinberg equilibrium were detected for any of the four SNPs (P > 0.05). Genotype and allele frequencies and population genetic parameters are summarized in Tables 5 and 6, respectively.

Table 5. Gene frequencies and genotype frequencies in BMP15 and GDF9.

Table 6. Population genetic structure analysis of BMP15 and GDF9 genes.

df=1, χ²_0.05 = 3.841; df=2, χ²_0.05 = 5.991.

Association analysis between single nucleotide polymorphism (SNP) loci and litter size was conducted in Tibetan sheep, and the results are presented in Table 7. No significant association with litter size was detected for the g.54284258G>A and g.54283966C>T loci in the BMP15 gene, nor for the g.42116345C>T and g.42115010C>T loci in the GDF9 gene. Both g.54283966C>T and g.42116345C>T harbored two genotypes, namely CC and CT. No significant difference in litter size was observed between CC and CT genotype individuals at either locus (P > 0.05). Similarly, no significant differences in mean litter size were detected among genotypes at the g.54284258G>A and g.42115010C>T loci (P > 0.05). The mean litter size, standard deviation (SD), and 95% confidence interval (95% CI) for each genotype are detailed in Table 7. Owing to the limited sample size of some rare genotypes, such as the GG genotype at g.54284258G>A (n = 3) and the CT genotype at g.54283966C>T (n = 7), all individuals within these genotype groups had identical litter size values, resulting in a standard deviation of 0 and zero-width 95% CIs. For other rare genotypes, such as the CT genotype at g.42116345C>T, the 95% CI was relatively wide, indicating limited robustness of the mean estimation for rare genotype groups. To further evaluate the magnitude of genotype effects, effect sizes were estimated based on the association analysis. The partial η² values for the four loci were 0.008, 0.007, <0.001, and 0.006, respectively, indicating very small genotype effects on litter size in the studied population. Together with the non-significant association results and the imbalanced genotype distribution, these findings suggest that the four SNP loci analyzed in this study were not significantly associated with litter size in the current Tibetan sheep population.

Table 7. Correlation analysis of polymorphisms and litter size in BMP15 and GDF9 genes.

3.4 Population genetic analysis of polymorphism in sheep BMP15 and GDF9

Linkage disequilibrium (LD) analysis was performed for the SNPs within the BMP15 and GDF9 genes, and the results are presented in Figure 10 and Table 8. One haplotype block was identified within each gene. For BMP15, the SNP pair g.54284258G>A and g.54283966C>T showed a D′ value of 1.000, a LOD score of 0.24, and an r² value of 0.001. For GDF9, the SNP pair g.42115010C>T and g.42116345C>T showed a D′ value of 1.000, a LOD score of 2.67, and an r² value of 0.039. Although the D′ values were high for both SNP pairs, the very low r² values indicated weak allelic correlation between the two loci within each gene. Therefore, these SNP pairs were not considered to exhibit strong LD in the studied Tibetan sheep population.

Table 8. Pairwise linkage disequilibrium parameters of SNPs within BMP15 and GDF9 genes in Tibetan sheep.

3.5 Haplotype analysis of BMP15 and GDF9 genes

Haplotype reconstruction was performed separately for the BMP15 and GDF9 genes based on their intragenic SNP loci. As summarized in Table 9, three haplotypes (H1–H3) were detected for the BMP15 gene, designated H1 (CA), H2 (CG), and H3 (TA). Haplotype H1 exhibited the highest frequency (0.9075) in the studied population, whereas haplotype H3 presented the lowest frequency (0.0175). For the GDF9 gene, three haplotypes (H1–H3) were also inferred, namely H1 (CC), H2 (TC), and H3 (TT). Haplotype H1 was the most prevalent at a frequency of 0.6667, while haplotype H3 was the least common (0.0167).

Table 9. Haplotype analysis of BMP15 and GDF9 genes.

4 Discussion

BMP15 and GDF9 are oocyte-derived members of the transforming growth factor-β (TGF-β) superfamily and act as important intraovarian regulators of follicular development and female fertility (Juengel et al., 2004; Knight and Glister, 2006). BMP15 is mainly produced by oocytes and exerts its biological functions by acting on surrounding granulosa and cumulus cells, thereby regulating granulosa cell proliferation, differentiation, and follicular development (Knight and Glister, 2006). GDF9 is also essential for normal follicle progression, as follicular development is arrested beyond the primary stage in GDF9-null mice and in ewes homozygous for naturally occurring inactivating mutations in the GDF9 gene (Dong et al., 1996; Hanrahan et al., 2004). Both BMP15 and GDF9 are synthesized as precursor proteins and subsequently processed into biologically active mature forms (Hashimoto et al., 2005; Juengel et al., 2004). Therefore, understanding the molecular structure and functional characteristics of BMP15 and GDF9 is important for elucidating their roles in folliculogenesis, ovulation, fertility, and litter size regulation in sheep.

In the present study, cDNA fragments of BMP15 and GDF9 were cloned from Tibetan sheep and subjected to sequence analysis and bioinformatic prediction. ORF prediction showed that the predicted ORFs within the cloned BMP15 and GDF9 fragments were 468 bp and 405 bp in length, encoding 155 and 134 amino acid residues, respectively. However, these predicted amino acid sequences were markedly shorter than the previously reported full-length ovine BMP15 and GDF9 precursor proteins. Therefore, the sequences obtained in the present study should be regarded as partial coding fragments rather than complete full-length coding sequences (CDSs) (Juengel et al., 2004). Accordingly, all subsequent bioinformatic analyses, including hydrophobicity, signal peptide prediction, transmembrane domain prediction, phosphorylation-site prediction, and secondary-structure prediction, were performed based on the predicted amino acid sequences encoded by these cloned fragments. Previous studies have shown that BMP15 and GDF9 are highly conserved across mammalian species and play important roles in ovarian follicular development and the regulation of female fertility (Galloway et al., 2000; Hanrahan et al., 2004; Juengel et al., 2004). To further evaluate the phylogenetic relationships between Tibetan sheep BMP15 and GDF9 fragments and their homologous sequences from other species, phylogenetic trees were constructed based on the predicted amino acid sequences. The results showed that the Tibetan sheep BMP15 fragment shared high sequence similarity with homologous sequences from other species, with similarity exceeding 90% in most species except Maylandia zebra. Similarly, the predicted Tibetan sheep GDF9 fragment showed high similarity to homologous sequences from most analyzed species, except Gallus gallus and Maylandia zebra, which exhibited relatively lower similarity. These results suggest that the analyzed regions of BMP15 and GDF9 are relatively conserved during evolution. Further bioinformatic prediction showed that no typical signal peptide was detected in the cloned BMP15 fragment. Considering that the BMP15 sequence obtained in this study represents only a partial coding fragment and may not cover the complete N-terminal region of the full-length BMP15 precursor protein, this result should not be interpreted as evidence against the secretory nature of the full-length BMP15 protein. For the cloned GDF9 fragment, signal peptide prediction suggested the presence of a Sec-type signal peptide, with a predicted cleavage site between the 25th and 26th amino acid residues, indicating that this fragment may retain partial sequence features of a secretory protein. Subsequent bioinformatic analyses showed that the protein fragment encoded by the cloned BMP15 sequence had an instability index of 35.21 and was therefore predicted to be stable, whereas the protein fragment encoded by GDF9 had an instability index of 53.62 and was predicted to be unstable. Hydropathy analysis indicated that the BMP15-encoded fragment was nearly neutral in overall hydropathicity, with slight hydrophobicity, whereas the GDF9-encoded fragment exhibited an overall hydrophilic character. TMHMM predicted one transmembrane-like segment at the N-terminus of the BMP15 protein fragment, whereas no transmembrane region was identified in the GDF9 protein fragment. The predicted secondary structures of both protein fragments were mainly composed of α-helices, random coils, extended strands, and β-turns. Overall, this study provides a preliminary characterization of the molecular features of the cloned BMP15 and GDF9 gene fragments in Tibetan sheep. However, because only partial coding fragments rather than complete full-length coding sequences were obtained, the protein structure and function predictions remain subject to certain limitations. Future studies should include full-length CDS cloning, protein-level validation, and in vitro functional assays to comprehensively clarify the complete molecular structures and biological functions of BMP15 and GDF9 proteins in Tibetan sheep.

Previous studies have shown that BMP15 and GDF9 are closely associated with follicular development, oocyte maturation, ovulation, and reproductive performance in sheep and goats (Galloway et al., 2000; Hanrahan et al., 2004; Knight and Glister, 2006; Tang et al., 2018; Yang et al., 2012). In the present study, the mRNA expression of BMP15 and GDF9 was examined in 12 tissues of Tibetan sheep, including the hypothalamus, heart, lung, kidney, liver, spleen, muscle, uterus, ovary, oviduct, rumen, and duodenum. Among these tissues, GDF9 showed the highest expression level in the ovary, which is consistent with its established role as an oocyte-derived factor involved in follicular development and oocyte maturation, suggesting that it may participate in the regulation of ovarian reproductive function in Tibetan sheep. Notably, BMP15 exhibited the highest mRNA expression level in the oviduct, with a relatively high expression signal also detected in the liver. Previous tissue-expression studies in goats have shown that GDF9 can be expressed in multiple tissues, whereas BMP15 generally exhibits a more tissue-restricted expression pattern, suggesting that the expression profiles of these genes may vary among species or breeds (Pan et al., 2015). The oviduct plays important roles in gamete transport, fertilization, early embryonic development, and embryo–maternal communication (Li and Winuthayanon, 2017; Kölle et al., 2020). Therefore, the relatively high expression of BMP15 observed in the oviduct may be related to the local reproductive tract microenvironment and provides preliminary evidence of a potential extra-ovarian expression pattern. However, this observation remains exploratory and cannot establish a specific biological function of BMP15 in the oviduct. The biological significance of the relatively high BMP15 mRNA signal detected in the liver also remains unclear and should be regarded as a preliminary and unexpected observation. Because the present study assessed only relative mRNA expression and did not include qPCR amplicon sequencing, systematic evaluation of reference-gene stability across tissues, or protein-level and spatial-localization validation, the current evidence is insufficient to confirm actual BMP15 protein expression or biological function in the oviduct, liver, or other extra-ovarian tissues. Future studies should include larger sample sizes and validate the mRNA-expression results using independently designed primers, amplicon sequencing, and multiple validated reference genes. Western blotting, immunohistochemistry, and in situ hybridization should also be performed to further characterize BMP15 protein abundance and spatial localization in extra-ovarian tissues.

Previous studies have shown that BMP15 and GDF9 are important candidate genes involved in follicular development, ovulation rate, and litter size regulation in sheep and goats (Galloway et al., 2000; Davis, 2004; Hanrahan et al., 2004; Pramod et al., 2013). In particular, mutations in BMP15 have been reported to increase ovulation rate in a dosage-sensitive manner and to cause infertility in homozygous carriers in sheep (Galloway et al., 2000). In Cambridge and Belclare sheep, mutations in BMP15 and GDF9 have also been associated with increased ovulation rate, and homozygous mutation carriers may exhibit infertility (Hanrahan et al., 2004; Davis, 2004). Therefore, these two genes are widely regarded as important candidate genes affecting reproductive traits in sheep.

Two single-nucleotide polymorphisms (SNPs) were identified in each of the candidate genes BMP15 and GDF9. The BMP15 gene is located on the X chromosome, and both g.54283966C>T and g.54284258G>A are located in exonic regions. Standardized VEP annotation showed that g.54283966C>T is a synonymous variant, changing the codon from GTG to GTA without altering the valine residue at position 348 (p.Val348=). The original locus label g.54284258G>A corresponds to g.54284258A>G on the forward strand of the reference genome and is a missense variant. This substitution changes the codon from CTG to CCG, resulting in a leucine-to-proline substitution at amino acid position 251 (p.Leu251Pro). The GDF9 gene is located on chromosome 5, and both g.42116345C>T and g.42115010C>T are located in exonic regions. The g.42116345C>T variant is a missense variant that changes the codon from CGC to CAC, resulting in an arginine-to-histidine substitution at position 87 (p.Arg87His). The original locus label g.42115010C>T corresponds to g.42115010T>C on the forward strand of the reference genome and is a synonymous variant. This substitution changes the codon from CTA to CTG without altering the leucine residue at position 159 (p.Leu159=). SIFT predicted both p.Leu251Pro and p.Arg87His to be tolerated, with scores of 0.16 and 0.12, respectively. Association analysis showed that none of the four SNPs was significantly associated with litter size in the Tibetan sheep population examined in this study. The absence of significant associations may be attributable to several factors. First, Tibetan sheep generally have relatively low prolificacy, and most ewes produce a single lamb per parity. The resulting limited phenotypic variation in litter size may have reduced the statistical power to detect genotype–phenotype associations. Second, although 200 ewes were included, this sample size may still have been insufficient to detect SNPs with small genetic effects, particularly when some genotypes occurred at low frequencies. In addition, although g.54284258G>A and g.42116345C>T were classified as missense variants after standardized annotation, both were predicted by SIFT to be tolerated. These loci may therefore not be causal variants directly affecting litter size, or they may not be in strong linkage disequilibrium with functional variants influencing prolificacy. Accordingly, the present findings indicate only that these four SNPs were not significantly associated with litter size in the Tibetan sheep population and sample examined in this study.

Linkage disequilibrium (LD) refers to the non-random association of alleles at different genetic loci and is widely used in population genetics and trait association analyses (Oraguzie et al., 2007; Slatkin, 2008). Pairwise LD indices, such as D^/and r², are commonly used to evaluate the strength of LD between markers (Liu et al., 2008). Haplotype analysis can integrate genetic information from multiple linked SNPs and may therefore provide more genetic-effect information than single-marker association analysis (Clark, 2004; Liu et al., 2008). In the present study, haplotype analyses were performed based on the SNP loci within the BMP15 and GDF9 genes. Three haplotypes were inferred for BMP15, namely H1 (CA), H2 (CG), and H3 (TA), and three haplotypes were inferred for GDF9, namely H1 (CC), H2 (TC), and H3 (TT). However, although D′ reached 1.0 for both SNP pairs, the corresponding r² values were very low (0.001 for BMP15 and 0.039 for GDF9), indicating weak allelic correlation between the loci. Therefore, the two SNPs within each gene were not considered to be in strong LD in the studied population.

Overall, the present results only indicate that the four SNP loci and their haplotypes analyzed in this study were not significantly associated with litter size in the experimental Tibetan sheep population. However, the possible effects of other variants within these two genes or genes involved in related reproductive pathways on litter size cannot be excluded. Further studies with larger sample sizes, validation in multiple populations, and in vivo and in vitro functional assays are needed to further identify and validate reliable molecular markers for litter size selection in Tibetan sheep.

5 Conclusion

The predicted open reading frames (ORFs) of the cloned Tibetan sheep BMP15 and GDF9 gene fragments were 468 and 405 bp in length, encoding 155 and 134 amino acids, respectively. Because these sequences did not cover the complete coding regions of ovine BMP15 and GDF9, they were considered partial coding fragments. Bioinformatic analyses predicted that the BMP15-encoded fragment was stable, whereas the GDF9-encoded fragment was unstable. Hydropathy analysis indicated that the BMP15 fragment was slightly hydrophobic, whereas the GDF9 fragment was hydrophilic overall. One transmembrane-like segment was predicted at the N-terminus of the partial BMP15 fragment, whereas a Sec-type signal peptide was predicted in the partial GDF9 fragment. The predicted amino acid sequences of both cloned fragments showed high similarity to the corresponding ovine reference sequences. Because these predictions were based on partial coding fragments, their structural and functional significance requires further validation using full-length coding sequences and protein-level experiments. Exploratory tissue-expression profiling showed that BMP15 transcript abundance was highest in the oviduct, whereas GDF9 expression was highest in the ovary. Three major haplotypes were identified for each gene, and no strong linkage disequilibrium was detected between the two SNPs within either BMP15 or GDF9. Association analysis showed that none of the four SNPs examined was significantly associated with litter size in the Tibetan sheep population included in this study. Overall, the partial sequence characteristics, exploratory tissue-expression profiles, and population-genetic data obtained in this study provide baseline information for further investigation of the potential roles of BMP15 and GDF9 in Tibetan sheep reproduction and for the future identification and experimental validation of reproduction-related functional variants.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fanim.2026.1904014/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Footnotes

Edited by: Yuanfang Li, Harbin Institute of Technology, China

Reviewed by: Linkai Zhu, Cornell University, United States

Hendro Sukoco, Universitas Sulawesi Barat, Indonesia