Molecular Characterization of TGF-Beta Gene Family in Buffalo to Identify Gene Duplication and Functional Mutations

The TGF-β superfamily is ubiquitously distributed from invertebrates to vertebrates with diverse cellular functioning such as cell adhesion, motility, proliferation, apoptosis, and differentiation. The present study aimed to characterize the TGF-β gene superfamily in buffalo through evolutionary, structural, and single nucleotide polymorphism (SNPs) analyses to find the functional effect of SNPs in selected genes. We detected 32 TGF-β genes in buffalo genome and all TGF-β proteins exhibited basic nature except INHA, INHBC, MSTN, BMP10, and GDF2, which showed acidic properties. According to aliphatic index, TGF-β proteins were thermostable but unstable in nature. Except for GDF1 and AMH, TGF-β proteins depicted hydrophilic nature. Moreover, all the detected buffalo TGF-β genes showed evolutionary conserved nature. We also identified eight segmental and one tandem duplication event TGF-β gene family in buffalo, and the ratio of Ka/Ks demonstrated that all the duplicated gene pairs were under selective pressure. Comparative amino acid analysis demonstrated higher variation in buffalo TGF-β gene family, as a total of 160 amino acid variations in all the buffalo TGF-β proteins were detected. Mutation analysis revealed that 13 mutations had an overall damaging effect that might have functional consequences on buffalo growth, folliculogenesis, or embryogenesis.


Introduction
Transforming growth factor beta (TGF-β) is a diverse gene family that contains a variety of growth factors, and all the members of this family have a set of three isoforms, TGF-1, TGF-2, and TGF-3, which are made up of interconnected dimeric polypeptide chains [1,2]. The TGF-β superfamily is ubiquitously present in both invertebrates and vertebrate species and plays a vital role in dorsoventral modeling, mesodermal initiation and patterning, and also in limb bud development, neuronal differentiation, and bone and cartilage formation [3][4][5]. The functional diversity of the TGF-β superfamily is crucial for the development of different body tissues and organs, particularly in vertebrates [6]. Furthermore, these proteins have a critical role in regulating and mediating basic cellular processes including cell motility, adhesion proliferation, apoptosis, and differentiation, as well as processes at the tissue or organism level, such as angiogenesis, growth, propagation, wound healing, and fibrosis [7].
The TGF-β superfamily members play a significant role in controlling the gene expression [8], and they can also regulate the expression of noncoding RNAs, such as microRNAs (miRNAs) and long noncoding RNAs (lncRNAs) [9]. In vertebrates, more than 30 genes belonging to the TGF-β superfamily have been reported so far, including TGF-β isoforms (TGF-β1, TGF-β2, and TGF-β3), bone morphogenetic proteins (BMPs), growth differentiation factors (GDFs), activins, inhibins, nodals, and Müllerian inhibitory factor (MIF) [7]. The TGF superfamily functions via Smads, which are well-defined downstream mediators. Smads control the gene expression either by activating or repressing a gene by interacting with high-affinity DNA-binding transcription factors and transcription co-regulators [8].
Members of the TGF-gene family have a key role in bovine physiology, for example, bone morphogenetic protein 1 (BMP1) contributes in the selection and dominance of follicles by regulating the proliferation and apoptosis of granulosa cells [10]. Similarly, bone morphogenetic protein/suppressor against decapentaplegic (BMP4/SMAD) signaling pathway has been suggested to play a role in regulation of follicular development by regulation granulosa cells in bovines [11]. Further, the bone morphogenetic protein 4 (BMP4) knockdown in bovines cumulus cells inhibited the proliferation of cumulus cells, apoptosis, and cell cycle arrest [12]. Growth differentiation factor-9 (GDF9) and bone morphogenetic protein 15 (BMP15) are also important candidate genes for oocyte maturation and embryo development and they were found to be highly expressed in the oocytes and embryos of buffalo [13]. Similarly, anti-Müllerian hormone (AMH) is also important for development and proper functioning of corpus luteum in buffalo [14]. The Inhibin Subunit Beta A (INHBA) gene was characterized in buffalo bulls and single amino acid variations detected at cleaving sites with potential association with growth, maintenance, and reproduction [15]. Moreover, myostatin (MSTN) is involved in skeletal muscle growth, and additionally was also found to be involved in folliculogenesis [16]. Water buffalo is an economically important genetic asset that contributes more than 15% of the world's total milk supply [17]. Despite having huge commercial importance, buffaloes are being neglected for their genetic resources, especially the breeding-and physiology-related regulators. No comparative genomic studies are available on buffalo and cattle to investigate poor estrus expression and lower reproductive efficiency in buffalo compared to cattle. Therefore, there is a dire need to shed light on the biological entities that would help to develop our understanding with different regulatory gene families which ultimately benefit buffalo development. The present study aimed to explore the evolutionary, physicochemical characterization, and gene structure analyses of the TGF-β superfamily in buffalo. Further, we also conducted a single nucleotide polymorphism (SNPs) analysis to find the functional effect of SNPs in selected genes of the TGF-β superfamily in buffalo.

Identification of TGF-β Genes in Buffalo
Whole-genome, proteome, and annotation data of Mediterranean river buffalo (UOA_WB_1), cattle (ARS-UCD1.2), sheep (Oar_rambouillet_v1.0), goat (ARS1), human (GRCh38.p12), and horse (EquCab3.0) were downloaded from National Center for Biotechnology Information (NCBI) Genome database (https://www.ncbi.nlm.nih.gov/) (accessed on 2 March 2022). Both the Basic Local Alignment Search Tool (BLAST) and hidden Markov model (HMM) searches were performed to identify all TGF-β protein isoforms in buffalo at the genome level. The non-redundant TGF-β gene sequences of cattle (Bos taurus), human (Homo sapiens), goat (Capra hircus), horse (Equus caballus), and sheep (Ovis aries), were retrieved from the UniProt (https://www.uniprot.org/) (accessed on 2 March 2022), and subjected as a query via BLASTp with a threshold of e-value = 10 −5 by using BLOSUM62 matrix with a six word size, eleven gap cost with an extension of 1, and a conditional composition score matrix adjustment. Additionally, the buffalo dataset was also searched with HMMER [18,19] (http://hmmer.org/) (accessed on 3 March 2022) software using HMM profile of the TGF-β domain (PF00019) from the Pfam online database [20] with an E-value 1.0 × e −5 . To avoid ambiguity, duplicate sequences were deleted after retrieving the relevant protein sequences. To confirm the TGF-β domains in protein sequences, these non-redundant sequences were analyzed in Simple Modular Architecture Research Tool (SMART) (http://smart.embl-heidelberg.de/) (accessed on 4 March 2022), and NCBI-CDD database was used for buffalo TGF-β proteins conserved domains searches.

Multiple Sequence Alignment
To identify the indels and visualize sequence variations, all the TGF-β protein sequences were aligned using Multiple Align Show (https://www.bioinformatics.org/sms/ multi_align.html) (accessed on 7 March 2022).

Structural Features Analysis
The conserved motifs were evaluated in the MEME suite (https://meme-suite.org/ meme/tools/meme) (accessed on 8 March 2022) [22]. All the buffalo TGF-β protein sequences were given in FASTA format as query. Site distribution was selected as one occurrence per sequence to find 10 MEME motifs with the minimum and maximum motif widths ranging between 6 and 50, respectively. To examine the pattern of introns and exons in TGFβ genes, CDs and genomic sequences were loaded in Gene Structure Display Server (GSDS) (http://gsds.gao-lab.org/) (accessed on 10 March 2022) and then the Tbtools (v1.098721) software (https://github.com/CJ-Chen/TBtools) (accessed on 10 March 2022), which used in-house scripts general feature format (GFF) file to depict the final gene structure [23].

Phylogenetic Analysis
All the TGF-β amino acid sequences of cattle, buffalo, goat, sheep, horse, and human were aligned in ClustalW, and the neighbor-joining (NJ) molecular phylogenetic tree using MEGA7 v.7.0 software (https://megasoftware.net/) (accessed on 10 March 2022) was constructed with a bootstrap value of 1000 replicates, adopting the Poisson model with pairwise deletion [24].

Synteny and Gene Duplications Analysis of TGF-β Superfamily Genes
Chromosomal locations of buffalo TGF-β genes were acquired from their genome resources, and a genome annotation (GFF) file was given as an input to the MCScanX program to map the physical locations of genes on chromosomes and then visualized in TBtools. Furthermore, the buffalo and cattle dual synteny plots were aligned for TGF-β genes collinearity [25]. Additionally, pairwise alignment of homologous gene pairs of TGF-β genes using MEGA7 v.7.0 with the MUSCLE algorithm was used to assess the occurrences of duplications for the buffalo TGF-β gene family [24]. DnaSP v6.0 software (http:// www.ub.edu/dnasp/) (accessed on 10 March 2022) was also used to estimate pairwise synonymous substitutions per synonymous site (Ks) and nonsynonymous substitutions per nonsynonymous site (Ka) adjusted for multi hits [26].

Genomic Identification of Buffalo TGF-β Genes
By using the BLAST and HMMER software, 193 non-redundant protein sequences encoded by 32 TGF-β genes were detected from the whole genomes of five mammalian species (cattle, buffalo, sheep, goat, and horse) while the human sequences were encoded by 33 TGF-β genes. Similarly, about 32 genes of the TGF-β gene family were also identified from the buffalo genome ( Figure 1).

Gene Structure Analysis of Buffalo TGF-β Genes
The phylogenetic relationship, motif distribution, gene structure, and conserved region analyses were conducted to discover the structural features of TGF-β gene family in buffalo, as shown in Figure 2A-D. These analyses were executed while taking their phylogenetic evolutionary relationships ( Figure 2A). From buffalo TGF-β genes, a total of 10 MEME conserved motifs were identified, of which motifs MEME-8 and MEME-9 have a higher number of amino acids, 49 and 45, respectively, and both the motifs were annotated as TGF-β domain after the Pfams search ( Figure 2B and Table 1). Additionally, the results were also checked against the NCBI-CDD database for verification ( Figure 2C). Additionally, along with the domain of TGF beta, the domain of AMH-N and PTZ00449 superfamily were also identified. Moreover, the gene structure analysis showed that the introns and the upstream and downstream untranslated regions (UTRs) structure varied greatly and all the buffalo TGF-β genes had different numbers of exon and intron patterns ( Figure 2D). higher number of amino acids, 49 and 45, respectively, and both the motifs were annotated as TGF-β domain after the Pfams search ( Figure 2B and Table 1). Additionally, the results were also checked against the NCBI-CDD database for verification ( Figure 2C). Additionally, along with the domain of TGF beta, the domain of AMH-N and PTZ00449 superfamily were also identified. Moreover, the gene structure analysis showed that the introns and the upstream and downstream untranslated regions (UTRs) structure varied greatly and all the buffalo TGF-β genes had different numbers of exon and intron patterns ( Figure 2D).

Characterization of Physicochemical Properties of Buffalo TGF-β Genes
Physiochemical attributes of the buffalo TGF-β gene family, including TGF-like and BMP-like, were analyzed for their chromosomal distribution, exon count, number of the amino acids (A.A) in each peptide, molecular weight (Da), isoelectric point (pI), instability index (II), aliphatic index (AI), and Grand Average of hydropathicity Index (GRAVY), as presented in Tables 2 and 3 (Tables 2 and 3), and the aliphatic index values were found to be >65, which exhibited thermostable characteristics of all the buffalo TGF-β proteins. According to the instability index, all the proteins appeared unstable as values were >40 (Tables 2 and 3). Owing to lower GRAVY values, all the buffalo TGF-β proteins had hydrophilic nature except for GDF1 and AMH, which depicted hydrophobic nature (Tables 2 and 3).

Comparative Amino Acid and Functional Mutation Effect Analysis of Buffalo TGF-β Gene Family
Comparative amino acid analysis for all the TGF-β genes of buffalo was conducted using cattle as a reference (Figures S2-S33). Four buffalo TGF-β genes, including TGFB2, TGFB3, INHBA, and BMP7, shared 100% homology with cattle genes (Figures S2-S5), while the indels were also observed in buffalo and cattle INHBE, GDF6, BMP6, TGFB1, and GDF11 genes, where four indels in each of INHBE and GDF6, three in BMP6, and one in each of TGFB1 and GDF11 were assessed (Figures S6-S10). Furthermore, a total of 160 amino acid alterations were calculated in all the buffalo TGF-β genes, where a single amino acid variation was detected in each BMP2, BMP4, and BMP10 gene (Table 5). Similarly, TGFB1, INHBB, BMP5, GDF1, GDF5, and GDF11 had two amino acid changes but two genes (GDF6 and GDF15) and four genes (MSTN, BMP6, BMP8B, and GDF3) with three and four amino acids differences, respectively, were also found (Table 5). Moreover, GDF7 with six, INHBC, INHBE, NODAL, and GDF10 with seven, only BMP3 with eight, INHA, GDF2, and GDF9 with nine, and BMP8A and BMP15 with ten amino acid substitutions were also evaluated in buffalo (Table 5). In addition, higher ratio of non-synonymous nucleotide substitution was perceived in LEFTY2 and AMH genes, with 16 and 18 amino acids replacements in buffalo (Table 5 and Figures S6-S33). Additionally, the functional effect of these mutations was assessed using different software (Table 5) and a total of 13 detected amino acid mutations in different TGF-β genes of buffalo, including MSTN (E116 > D), , and AMH (A334 > T), due to nucleotides non-synonymous alterations were supposed to have overall damaging effect on protein structure and functions, while the other amino acid substitutions have overall neutral affect and caused no serious influence on protein structure and functions (Table 5).

Discussion
The advancements in high-throughput genome sequencing technologies, typically the next-generation sequence data availability, make it easy to scan the genetic variability of genes, such as SNPs ordering with their functional effect which control a specific phenotypic trait that allow to understand animals' genetics at molecular level [27][28][29][30][31][32]. For farm animals, the candidate gene studies analyze available genetic resources to identify functional genes and their potential association with productivity traits, such as disease resistance, production ability, and adaptation [29]. In buffaloes, comparative genomics offers an opportunity to explore the genetics of economically important physiological traits though discovering novel genes and their regulatory mechanisms, which significantly contributes to the development of the buffalo industry [32].
The TGF-gene family has a set of coding genes with conserved structure and biologically diverse functioning that release signaling molecules which regulate fundamental cellular pathways, and dynamic biological processes such as apoptosis, communication, proliferation, differentiation, and tissue remodeling throughout growth, repair, and organo-genesis [33,34]. In the present study, 32 TGF-β coding genes were identified from buffalo genome that were categorized into two major groups (TGF-β-like and BMP-like), where the TGF-β-like had six while BMP-like had only five sub-groups or set of genes sharing higher sequence homology with Bos taurus than Capra hircus and Ovis aries (Figure 1). The phylogenetic pattern observed in our study is in agreement with earlier findings reporting the evolution of the TGF-gene family in 9 invertebrates and 15 chordates species [35]. Furthermore, TGF-β peptide was the only conserved motif found in the buffalo TGF-gene family ( Table 1).
All the members of the TGF-β superfamily are predominantly involved in a variety of cellular activities by binding to a particular receptor to generate a signal for a specific cellular function. Except for inhibin and GDNF subfamilies, all the other TGF-superfamily ligands bind to a unique set of double transmembrane Ser/Thr kinase receptors, type I and type II receptors [36]. Firstly, the dimeric ligand binds to cognate type II receptors and then this complex recruits or activates type I receptors, eventually causing R-SMAD phosphorylation and activation (receptor-regulated SMADs) [37]. Generally, there are two models for R-SMAD activation, firstly, SMAD1/5/8 interact with AMH, BMPs, and some GDFs through ALK-2, -3, and -6; secondly, SMAD2/3 respond to activins, NODAL TGF-βs, and some GDFs via ALK-4, -5, and -7. Finally, the activated R-SMADs form a complex with SMAD (Co-SMAD) or SMAD4, and then this complex could translocate into the nucleus and influence or regulate the target gene expression through interacting with other transcription factors in different cell types or tissues [37,38].
The GDF1 and GDF3 could interact with nodal and initiated signals for the long term, which are important in left-right patterning [39], while GDF3 is substantially involved in robust nodal signaling during germ layer formation [40]. Further, AMH which is distributed in Sertoli cell of testes, have functional involvement in gonads development and sex differentiation through inhibiting the Müllerian duct by smad1 pathway [41]. Moreover, the GDF8 (myostatin) and GDF11 act through type II receptors (ActRII or ActRIIB) or type I receptors (ALK-4 or -5 type) that activate the Smad2/3 pathway and regulate the muscle mass by inhibiting the muscle differentiation and regeneration [42][43][44][45]. It is unlikely that GDF5 is involved in muscle hypertrophy through BMP signaling [46]. The BMP9 and BMP10 are recognized as ligands of ALK-1 receptor which crucially regulate vascular system homeostasis and heart development [47][48][49], while BMP6 controls the iron metabolism and expression of hepcidin by acting as an endogenous ligand [50,51].
Similarly, the GDF9 and BMP15 are functionally expressed in oocytes and work synergistically through the smad3 pathway via ERK1/2 and SRC kinase-dependent signaling [52]. Although GDF10 and BMP3 are closely related and have positive role in endochondral bone formation in adult animals, they negatively affect the bone morphogenesis at embryonic stage [53]. Additionally, BMP8A is involved in the regulation of spermatogenesis through smad2/3 and smad1/5/8 pathways and BMP8B prevents the male germ cells apoptosis [54,55], whereas GDF5/6/7 contributes to the formation of normal bones, limb joints, skull, and axial skeleton [56].
Organisms use gene duplication mechanisms including retroposition, genome, or chromosomal duplication, and crossing over to acquire novel gene or genetic variants, which tremendously contributes to the evolution of functional processes [57]. Discovering the dynamics that create duplicate gene copies, as well as the subsequent trajectories among duplicated genes, is critical because these studies shed light on localized and genome-wide attributes of evolutionary factors that influence intra-and inter-specific genome contents, evolutionary interactions, and relationships [57]. It is difficult to measure the rate of duplications that occur but selective pressure and mutations with functional effects are vital for developing redundant genetic variants [58]. Over generations of an organism, the duplicated gene could accumulate mutations faster than in a single functional gene copy and possibly develop a novel function [59]. Earlier, it has been reported that in ice fish, the apparent mutations in a duplicated digestive gene transform to the antifreeze gene while duplication leads to a distinct snake venom gene [59], and in pigs, synthesis of 1-β-hydroxytestosterone [60]. In the present study, we identified eight segmental and one tandem duplication event in the buffalo TGF-β gene family. The Ka/Ks ratio for duplicated gene pairs demonstrated that all these gene pairs were under selective pressure, among which four gene pairs were under positive selection with Ka/Ks ratio > 1, while five gene pairs had Ka/Ks ratio < 1, exhibiting that these are purifying selection pressure [61]. Our findings are consistent with earlier studies and suggest that observed gene duplication in the TGF-β gene family resulted in an increased buffalo genome size and complexity [62].
In mammals, seven receptors belonging to type I called activin receptor (kinase-like 1-7) and five type II receptors (such as ACVR2A, ACVR2B, AMHR2, BMPR2, and TβR2) have already been identified that can induce the heterotetrameric complex [38]. In this reverence, the sequence and structural variations in TGF-β superfamily ligands should be crucial for their differential binding affinities for distinct type I and type II receptors. Furthermore, the type II receptor could also influence the binding affinity of ligands with type I receptors [63]. The comparative amino acid analysis of our study presented a total of 160 amino acid substitutions in buffalo TGF-β genes, of which 13 exhibited an overall damaging effect on respective protein structure and functions. Except for GDF3, GDF9, and BMP15, all TGF-β superfamily ligands use an extra conserved cysteine for intermolecular disulphide bond formation, which stabilizes the dimer [64]. Even though the majority of ligands tend to be homodimers, ligand heterodimerization has also been reported [64], such as heterodimerization between GDF9 and BMP15 [65], GDF1 and NODAL [39], and activin βA-βB [64].
Likewise, the members of the TGF-β superfamily play a dynamic role in controlling different biological processes such as folliculogenesis [66], skeleton development [67], nodal signaling [40], and fat yield [68], while absence or altered expression of these genes can induce impaired development. For example, GDF3 is important for folliculogenesis [66], skeleton development [67], and nodal signaling [40], but zebrafish mutant embryos showed early embryonic lethality [40]. Similarly, in humans it has been reported that mutations in GDF3 cause ophthalmic and skeletal abnormalities, and the patients with ocular and/or skeletal (Klippel-Feil) anomalies were found to have several missense variations, including one with heterozygous changes in GDF3 and GDF6 [69]. Additionally, GDF6, which is important for proper development of bones and joints in the limbs, head, and axial skeleton but abnormal development of joints, ligaments, and cartilage, could occur in the absence of the GDF6 gene [56], and GDF6 mutations in zebrafish have been reported associated with reduced eye size and different skeletal defects [70].
On the other hand, GDF9 is primarily involved in folliculogenesis and controls the development of theca and granulosa cells, as well as plays an important role in oocyte differentiation and maturation [71,72]. In sheep, GDF9 natural mutations resulted in higher ovulation rate and twin or triplet births in heterozygotes, while in homozygotes a complete primary ovarian failure leading to complete sterility was also reported [73]. Furthermore, the polymorphism in the GDF9 gene in sheep was also associated with litter size, milk production, and prolificacy [74,75]. In buffalo, GDF9 plays a significant role in buffalo oogenesis as it is present throughout folliculogenesis and embryogenesis stages, while during follicular development, the GDF9 gene surged at first, then reduced [76]. Moreover, alterations in the maternal mRNA transcript of GDF9 gene in buffalo oocytes resulted in the deleterious seasonal effects on oocytes development [77]. In the present study, we identified a total of nine amino acid substitutions, of which one belongs to GDF3, three to GDF6 along with four amino acid deletion, and five associated with GDF9. Further, the functional mutation analysis revealed that only two mutations of GDF9 exhibited damaging effect, which might have functional influence on folliculogenesis or embryogenesis of buffalo. The present study provides evolutionary insights into the TGF-β gene superfamily in buffalo, revealing the functional importance of SNPs in different genes with developmental and physiological consequences on buffalo performance. The mutations identified in the present study provide a basis for further studies to investigate their potential utility for genomic selection for targeted improved breeding and utilization of buffalo genetic resources.

Conclusions
In the present study, we comprehensively described the molecular structure and functional role of mutations in the TGF-β superfamily in buffalo. All the TGF-β genes in buffalo showed evolutionary conserved nature. Eight segmental and one tandem duplication events in the buffalo TGF-β gene family were identified, and all the duplicated gene pairs were under selective pressure. These duplications might have played a role in adaptation and speciation to specific functions and ecological niches [78]. We identified 13 mutations in different TGF-β genes with an overall damaging effect that might have functional consequences on buffalo growth and development, folliculogenesis, and/or embryogenesis. The present study is the first report on comparative genome-wide characterization of the TGF-β superfamily in buffalo and it provides insight into the evolutionary importance of gene duplications and mutations in the TGF-β superfamily with respect to adaptation and speciation. Furthermore, these findings would be helpful in understanding the crucial role of these genes and their potential utility for selective breeding in buffalo for economically important traits such as reproduction, growth, and development. Further studies are warranted to elucidate functional effects of mutations identified in different TGF-β genes in the present study to confirm their physiological manifestation in buffalo and potential effects on growth, development, and reproductive performance in buffalo.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/genes13081302/s1, Figure S1: Molecular phylogenetic analysis of all the buffalo TGF-β genes for putative duplication detection; Figure S2: Comparative amino acid analysis of buffalo and cattle TGFB2 gene; Figure S3: Comparative amino acid analysis of buffalo and cattle TGFB3 gene; Figure S4: Comparative amino acid analysis of buffalo and cattle INHBA gene; Figure S5: Comparative amino acid analysis of buffalo and cattle BMP7 gene; Figure S6: Comparative amino acid analysis of buffalo and cattle INHBE gene; Figure S7: Comparative amino acid analysis of buffalo and cattle GDF6 gene; Figure S8: Comparative amino acid analysis of buffalo and cattle BMP6 gene; Figure S9: Comparative amino acid analysis of buffalo and cattle TGFB1 gene; Figure S10: Comparative amino acid analysis of buffalo and cattle GDF11 gene; Figure S11: Comparative amino acid analysis of buffalo and cattle BMP2 gene; Figure S12: Comparative amino acid analysis of buffalo and cattle BMP4 gene; Figure S13: Comparative amino acid analysis of buffalo and cattle BMP10 gene; Figure S14: Comparative amino acid analysis of buffalo and cattle INHBB gene; Figure S15: Comparative amino acid analysis of buffalo and cattle BMP5 gene; Figure S16: Comparative amino acid analysis of buffalo and cattle GDF1 gene; Figure S17: Comparative amino acid analysis of buffalo and cattle GDF5 gene; Figure S18: Comparative amino acid analysis of buffalo and cattle GDF15 gene; Figure S19: Comparative amino acid analysis of buffalo and cattle MSTN gene; Figure S20: Comparative amino acid analysis of buffalo and cattle BMP8B gene; Figure S21: Comparative amino acid analysis of buffalo and cattle GDF3 gene; Figure S22: Comparative amino acid analysis of buffalo and cattle GDF7 gene; Figure S23: Comparative amino acid analysis of buffalo and cattle INHBC gene; Figure S24: Comparative amino acid analysis of buffalo and cattle NODAL gene; Figure S25: Comparative amino acid analysis of buffalo and cattle GDF10 gene; Figure S26: Comparative amino acid analysis of buffalo and cattle BMP3 gene; Figure S27: Comparative amino acid analysis of buffalo and cattle INHA gene; Figure S28: Comparative amino acid analysis of buffalo and cattle GDF2 gene; Figure S29: Comparative amino acid analysis of buffalo and cattle GDF9 gene; Figure S30: Comparative amino acid analysis of buffalo and cattle BMP8A gene; Figure S31: Comparative amino acid analysis of buffalo and cattle BMP15 gene; Figure S32: Comparative amino acid analysis of buffalo and cattle LEFTY2 gene; Figure S33: Comparative amino acid analysis of buffalo and cattle AMH gene.