Investigation of GDF9 and BMP15 Polymorphisms in Mehraban Sheep to Find the Missenses as Impact on Protein

INTRODUCTION

Reproduction is a complex process and fecundity traits such as ovulation rate (OR) and litter size (LS) are genetically affected by many minor genes and also some major genes, called fecundity (Fec) genes (Drouilhet et al. 2009). Utilization of major fecundity genes in sheep production can cause genetic improvement in breeding programs, and consequently high performance in reproductive traits (Notter, 2008). So far, some major genes relevant to OR and LS traits such as bone morphogenetic protein recep-tor-1B(BMPR1B), bone morphogenetic protein 15 (BMP15), growth differentiation factor 9 (GDF9) and beta-1, 4-N-acetyl-galactosaminyl transferase 2 or lacaune gene(B4GALNT2) have been reported and all of which have hyper prolificacy-associated mutations (Drouilhet et al. 2013; Bodin et al. 2007). BMPR1B, BMP15 and GDF9 genes are all part of the ovary-derived transforming growth factor-β (TGFβ) superfamily. The essential growth factors and receptors in ovarian follicular development are coded by these genes (Pramod et al. 2013). Causative mutations in sheep such as FecB^B mutation in BMPR1B in Booroola Merino (Wilson et al. 2001) and Javanese (Bradford et al. 1986); mutations in BMP15 in Romney (FecX^I, FecX^H) (Galloway et al. 2000), Cambridge and Belclare (FecX^B, FecX^G) (Hanrahan et al. 2004), Lacaune(FecX^L, FecX^I) (Bodin et al. 2007) and mutations in GDF9 in Belclare and Cambridge (FecG^H) (Hanrahan et al. 2004), Norwegian White Sheep (FecG^NW) (Vage et al. 2013), Icelandic(FecG^T) (Nicol et al. 2009) and Santa Ineˆs sheep (FecG^E) (Silva et al. 2011), have been clearly identified and also their significant effect on fecundity traits have been showed in different sheep breeds. Several methods based utilization of the physico-chemical properties of amino acids, as well as information about the role of amino acid side chains in protein structure, have been developed for assessing the effects of mutation on protein function (Reva et al. 2011). Mutations can influence protein folding, protein stability, protein function, protein-protein interactions, protein expression and subcellular localization (Reva et al. 2011). This work examined the polymorphic sites of GDF9 and BMP15 genes in Mehraban sheep to identify missense mutations. In addition, protein folding and physical interactions were assessed through protein modeling to interpret missense effects on protein conformation and functionality.

MATERIALS AND METHODS

Sampling

To perform the experiment, twelve Mehraban ewes, located in Hamedan province, western Iran, were selected. Blood samples (5 mL per animal) were collected from jugular vein using venoject containing EDTA. Genomic DNA was extracted from whole blood using DNPTM Kit (SinaClon BioScience Co.) and kept at -20 ˚C. Quantity and quality of gDNA were determined using NanoDrop® ND-1000 based on absorbance at A260/A280 ratio.

PCR amplification and purification of products

Specific primers were designed by Primer-BLAST tool (www.ncbi.nlm.nih.gov/tools/primer-blast/). The primers were used to amplify both exons 1 and 2 of BMP15 (GenBank No. NC_019484.1) and GDF9 (GenBank No. NC_019462.1) genes (Table 1). Polymerase chain reaction (PCR) was carried out on 50 ng of genomic DNA in a final 20 μL reaction containing 4 μL of 10 × PCR buffer, 2 μL of 200 µM dNTPs, 1 μL of each primer at a concentration of 0.05 μM and 0.5 unit of Taq polymerase (Promega, Wisconsin, USA). The PCR protocol used and initial of 5 min denaturation step at 94 ˚C followed by 35 cycles at 94 ˚C for 30 s, 55 ˚C for 30s 72 ˚C for 30 s and a final 7 min extension step at 72 ˚C. Purification of PCR products was carried out in a final 15 μL volume consisting of mixture of T-SAP (1 U/µL) and Exonuclease I (20 U/µL) with Programs PCR as follows: digestion 37 ˚C for 45 minutes, inactivation of enzymes 80 ˚C for 30 minutes in 35 cycles.

DNA sequencing

All twelve samples were selected for DNA sequencing. The primers used for sequencing were the same primers used for the PCR. The sequencing reaction was carried out in a final 20 μL volume containing 2 µL Tampon ABI (5X), 0.5 µL Terminators Dichlororhodamine V3.1 (Life technologies Co. California, United States), 1 µL primer (10 µM), 1.5 µL H₂O. The program of PCR sequencing as follows: 95 ˚C for 5 min followed by 1 cycle, and then 25 cycles followed by 95 ˚C for 30 s, 55 ˚C for 30 s, 60 ˚C for 4 min. Sequencing was carried out at the National Institute of Agronomic Research (INRA), France and GenPhySE Centre on the Applied Biosystems 3730 DNA Analyzer platform.

Bioinformatics analysis of polymorphism effects on DNA and protein structures

CLC Genomics Workbench Version 7.6.4 (www.clcbio.com) was used to map all obtained read counts of GDF9 and BMP15 genes to the reference genomes of NC_019462.1 and NC_019484.1, respectively. Option of annotation was used with reference genome in order to find each conflict on the genome level. According to annotation analysis we carefully manifested every known polymorphisms based on other breeds and new polymorphisms based on identified conflicts. Analysis of protein characterizations was accomplished using Protean (DNASTAR Inc., Madison, WI. USA). The mutation scrutinizes in the secondary structure of proteins were envisaged using Protean (DNASTAR Inc., Madison, WI. USA). The Expasy website (http://us.expasy.org/) was used to constitute 3D structure of proteins based on homology modeling. Furthermore, the predicted structures were compared with 3D models of Protein Data Bank (www.pdb.org) TGF-ß 4YCI (Mi et al. 2015), then modeling of protein was carried out using modeler tools of CLC Genomics Workbench Version 7.6.4 based on protein structures of PDB and predicted structures. The sequences of amino acids for different mono-ovulatory and poly-ovulatory species were obtained from biological database such as Uniprot (http://www.uniprot.org/uniprot/). Polymorphism Phenotyping v2 (PolyPhen-2) was used to predict the possible effects of the observed missense mutations on the structure and function of proteins based on an iterative greedy algorithm (Adzhubei et al. 2010). To consider GDF9 and BMP15 in different species, at first, we decided to focus on the sequences of amino acids to achieve a general understanding about substitutions rate of amino acids. A maximum-likelihood method based on the JTT matrix-based model was applied to infer evolutionary history by CLC Genomics Workbench Version 7.6.4 (www.clcbio.com) (Jones et al. 1992).

Table 1 Primers used for exon 1 and 2 of GDF9 (NC_019462.1^*) and BMP15 (NC_019484.1^*) genes amplification

^* GenBank accession numbers at NCBI.

^a amplified fragment length.

The bootstrap consensus tree inferred from 1000 replicates was taken to represent the evolutionary history of the analyzed taxa.

RESULTS AND DISCUSSION

According to Figure 1, the results obtained by DNA sequencing showed six point mutations in exons 1 and 2 of GDF9 gene in comparison to the reference sequence (NC_019462.1). These mutations known as G1 (g.2118 G>A, NC_019462.1), G2 (g.3451 T>C, NC_019462.1), G3 (g.3457 A>G, NC_019462.1), G4 (g.3701 A>G, NC_019462.1), G5 (g.3958 A>G, NC_019462.1) and G6 (g.3974 G>A, NC_019462.1) were previously reported by Hanrahan et al. (2004) (Table 2). Furthermore, sequencing results of the BMP15 gene showed no mutations in comparison to the reference sequence (NC_019484.1). The individuals with simultaneous mutations are surely in complete linkage disequilibrium that it was called as a haplotype situation. According to notable mutations identified in GDF9 of Mehraban sheep, the individuals were classified in three haplotypes as, Wild haplotype (without mutation), haplotype A (simultaneous mutations of G1, G2, G3 and G4) and haplotype B (G5 and G6). It is clear that changing one nucleotide may change the amino acid that the nucleotides are coding for. As it is indicated in Table 2, only G1, G4 and G6 mutations resulted in amino acid substitution at p.Arg87His, p.Glu241Lys and p.Val332Ile positions, respectively. With changing amino acids after mutation occurrence, it impresses on the second structure of TGF-ß proteins. The antigenic index was applied to predict the topological features of GDF9. Only a significant antigenic index was observed in G1 (Arg:0.45>His:-0.3). It was not seen any changes in antigenic index of G4 (Glu:0.3>Lys:0.3) and G6 (Val:-0.3>Ile:-0.3). The degree of hydrophobicity or hydrophilicity of amino acids of GDF9 protein was quantified using hydrophilicity plot. The observed changes were belonged to G1 (Arg:1.63>His:1.49), G4 (Glu:0.422>Lys:0.467) and G6 (Val:-0.244>Ile:-0.278). It is assumed that with occurrence of missense mutations in GDF9, some chemical features such as aliphatic index, molecular weight, 1 microgram, isoelectric point, molar extinction coefficient, absorbance at 280^nm and charged at pH= 7 were affected in GDF9 protein. As illustrated in Figure 2, only changed amino acid position of G1/p.Arg87His was occurred on alpha helix region of GDF9. However, none of change position of G4/p.Glu241Lys (Figure 2A) and G6/p.Val332Ile (Figure 2B) did not occurred on regions of alpha helix, Beta sheets and so on. Because neither the BMP15 nor the GDF9 3-D structures were known (Monestier et al. 2014), the three types of haplotype (wild, A and B) were aligned to predict a model for GDF9 protein the three haplotypes (Figure 3). As illustrated in Figure 3, the predicted GFP9 protein structures was aligned (Figure 3A) with mouse BMP9 (PDB 4YCI by Mi et al. (2015)) for four accession numbers (GI:764091320, GI:764091321, GI:764091322 and GI:764091323) to manifest the turning and folding at protein structure after mutation occurrence. The results of alignment illustrated that the structure of GDF9 protein in haplotype A (G1/p.Arg87His and G4/p.Glu241Lys, simultaneously) and hplotype B (heterozygote for G6/p.Val332Ile) has been rotated 90˚ and 45˚ than wild haplotype, respectively. According to the structural alignment of Mehraban GDF9 with 4YCI chain-A and B (Figure 3B), only chain B was entirely aligned with predicted structures of GDF9 protein and a partial part of 4YCI-D was also aligned with predicted structures of GDF9 protein (Figure 3C). The perfect alignment of predicted structures of GDF9 protein of Mehraban ewes with mouse BMP9 (PDB 4YCI) (Figure 3D) can be inferred to responsible for physical and functional interaction between GDF9 and BMP9 proteins or other TGF-β members.

Figure 1 Identified polymorphisms of GDF9 exon 1 and exon 2 in twelve Mehraban sheep

In comparison of reference genome NC_019462.1, four ewes revealed haplotype A (simultaneous mutations of G1, G2, G3 and G4), one ewe has haplotype B (simultaneous mutations of G5 and G6) and seven animals showed wild haplotype

Table 2 Identified mutations in exon 1 and exon 2 of GDF9 gene

Figure 2 The predicted structures of secondary and three-dimensional of GDF9protein in Mehraban breed

A) Three dimensional (up) and secondary structure (down) in haplotype A (simultaneous mutations of G1, G2, G3 and G4

B) Three dimensional (up) and secondary structure (down) in haplotype B (simultaneous mutations of G5 and G6)

Of course in reason of the low sequence identities between BMP15 and GDF9 on the one hand and TGFB1 on the other hand, these results need to be interpreted with caution. Moreover, the potential impact of the identified amino acid substitution, in this study, on the straucture and functional of GDF9 was assessed using PolyPhen2. As shown in Figure 4 the amino acide substitutions of p.Arg87His, p.Glu241Lys and p.Val332Ile were benign with scores of 0.00, 0.00 and 0.031, respectively (Figure 4A and 4B). With attention to possible impact of the missense on protein function, it can be imagined that the six point mutations of GDF9 containing G1/p.Arg87His, G2/p.Val157Val, G3/p.Lue159Lue, G4/p.Glu241Lys, G5/p.Glu326Glu, G6/p.Val332Ile are benign. However, both amino acid substitutions of G7/p.Val371Met and G8/p.Ser315Phe were probably damaging with the scores of 0.99 and 1.00 respectively (Figure 4C and 4D). They were belonged to locus of FecG^NW and FecG^H, respectively (Vage et al. 2013; Hanrahan et al. 2004). Phylogenetic tree of GDF9 demonstrated that wild hapl type of Mehraban and Texel breed (O77681), human (O60383) and chimpanzee (H2QRG9), rat (Q9QYW4) and mouse (Q07105), dog (D0F1P5) and giant panda (C5HV43) were grouped together (Figure 5A). However, haplotypes of A and B distanced from Texel breed (O77681) with bootstrapping values of 63 and 76, respectively.

Figure 3 Alignment of predicted 3D structures of GDF9 protein via homology modeling

A) Alignment of predicted models of GDF9 protein for three type haplotype wild, A and B

B) Alignment of predicted models of GDF9 protein and mouse BMP9 (PDB 4YCI chain-A,B)

C) Alignment of predicted models of GDF9 protein and mouse BMP9 (PDB 4YCI chain-C,D)

D) Perfect structural alignment between predicted models of Mehraban GDF9 protein and mouse BMP9 (PDB 4YCI)

Figure 4 Impact of amino acid substitution on the structure and function of proteins

A) Impact of amino acid substitution (G1/p.Arg87His and G4/p.Glu241Lys) on GDF9

B) Impact of amino acid substitution (G6/p.Val332Ile) on GDF9

C) Impact of amino acid substitution (G7/p.Val371Met) on GDF9

D) Impact of amino acid substitution (G8/p.Ser315Phe) on GDF9

Figure 5 Molecular phylogenetic by maximum likelihood method

The tree is drawn to scale, with branch lengths (gray notes) and bootstrap values (bold and italic notes) measured in the number of substitutions per site

A) Circular cladogram of GDF9 protein in mono and poly ovulation species

B) Circular cladogram of BMP15 protein in mono and poly ovulation species

As well as, phylogenetic tree of BMP15 in different species showed that Mehraban breed and Texel breed (Q9MZE2), human (O95972) and western lowland gorilla (G3QUI0), rat (E9PTU9) and mouse (Q9Z0L4), dog (E2QX74) and giant panda (D215P3), olive baboon (A0A096MY95) and green monkey (A0A0D9RJW4), chicken (Q645R0) and common turkey (G1NBP2) were grouped together (Figure 5B). Therefore, evolutionary origin of GDF9 and BMP15 proteins represents mono- or poly ovulatory species are not affected by protein sequence. Functional effect of protein mutations is specified by the diversity of their impact on molecular function (Reva et al. 2011). GDF9 and BMP15 have the main role in the process of follicular development and oocyte maturation (Våge et al. 2013). It is reported eight point mutations (G1, G2, G3, G4, G5, G6, G7 and G8) occurring in exons 1 and 2 of GDF9 gene in Cambridge/Belclare sheep, while only five mutations including G1/p.Arg87His, G4/p.Glu241Lys, G6/p.Val332Ile, G7/p.Val371Met and G8/p.Ser315Phe, led to amino acid substitution (Hanrahan et al. 2004). Following to study of Hanrahan et al. (2004), G1 was never associated to prolificacy, neither in the Cambridge/Belclare breeds nor in other breeds (Hanrahan et al. 2004). Anyway the role of G1 polymorphism may be retained unknown G1 polymorphism was roughly illustrated in some of the Iranian native breeds such as Mehraban (present study and Abdoli et al. (2013)), Baluchi (Moradband et al. 2011), Moghani and Farahani (Potki et al. 2015), Lori (Zamani et al. 2015), Afshari and Shal (Eghbalsaied et al. 2014), Moghani and Ghezel (Barzegari et al. 2010). But there are conflicting results on significant association of G1 with prolificacy (Barzegari et al. 2010; Javanmard et al. 2011; Moradband et al. 2011; Eghbalsaied et al. 2014; Zamani et al. 2015; Abdoli et al. 2013), Eghbalsaied et al. (2014) indicated the existence of simultaneous/dual mutation (G1 and G4) in Afshari and Shal sheep breeds. Similar observations have been reported in Davisdale flock (Juengel et al. 2004). Among all the polymorphisms identified in GDF9, only five mutations that are FecG^H (also known as G8 in Hanrahan et al. 2004), FecG^T in Icelandic Thoka sheep (Nicol et al. 2009), FecG^E in Brazilian Santa Ines (Silva et al. 2011), FecG^NW (also known as G7) in White Norwegian sheep (Vage et al. 2013) and FecG^V in Brazilian Île-de-France (Souza et al. 2014) are proven to be associated with an altered activity of GDF9 that affect the ovarian function in European sheep breeds. Moreover, the BMP15 gene regulates proliferation and granulosa cells differentiation and promoting ligand expression, which plays an important role in female fertility in mammals. Also as an important role of BMP15 in early follicle growth is relating to mono- and poly-ovulatory animals (Moore and Shimasaki, 2005). Although any point mutations were not identified in BMP15 gene for Mehraban breed in the present study, but several missenses in BMP15 gene with positive effect on reproduction traits, have been identified in different sheep breeds (Galloway et al. 2000; Bodin et al. 2002; Hanrahan et al. 2004; Monteagudo et al. 2009). Phenotypically changes in ovarian function due to alterations in BMP15 and GDF9 genes emerge to differ among species and may be relevant to their mono- or poly ovulating status (Monestier et al. 2014). Screening of G8 mutation (FecG^H) in GDF9 and Booroola SNP in BMPR-1B in various Iranian sheep breeds including Shal, Ghezel, Baluchi and Sangsari, did not indicated the presence of these SNPs (Akbarpour et al. 2008; Ghaffari et al. 2009; Moradband et al. 2011; Kasiriyan et al. 2011). The ovulation rate observed in double heterozygotes (FecX^G/FecX⁺ and FecG^H/FecG⁺) reflects an essentially additive effect of these mutations on prolificacy. However, homozygous carriers of FecG^H or of either of the mutations in BMP15 and ewes with the genotype FecX^G/FecX^B, are sterile due to arrested follicle development (Mullen et al. 2013). Without any functional testing, we cannot conclude that the observed mutations of GDF9 in Mehraban sheep are associated with ovarian function. The direct prediction of a mutation’s impact on molecular function based on first principles is currently impossible for a number of reasons: e.g. lack of data (3D structures and complexes) and lack of accurate and efficient approaches for de novo modeling of protein structure and function on the molecular level. However, evolutionary analysis does provide a powerful tool, as natural selection of a particular sequence variant by definition reflects the aggregate effect of molecular changes on cell, tissue and organ physiology (Reva et al. 2011). In spite impact of mutation on performance is proved the based on association analysis with traits, but it must be proven through biological models as well. Demars et al. (2013), have been done a genome-wide association studies (GWAS) to identify genetic variants responsible for the highly prolific phenotype in two sheep flocks of Grivette and Olkuska, also they have done functional analyses based on altered the protein signaling activity in vitro. They identified two novel non-conservative BMP15 mutations, that both mutations altered the BMP15 signaling activity, suggesting a novel kind of BMP15 variant responsible for an atypical high prolificacy, in contrast to all other BMP15 variants described so far (Demars et al. 2013). Therefore, to acknowledge effect of mutations on protein functionality, it should be proven as statistical model with followed up by genotyping a large number of ewes with known reproductive performance, bioinformatics model (mathematical algorithms) and biological model (in vitro functional assay). Despite the identified mutations in current article were benign at alone. But the based on protein modeling, simultaneous mutations in genetic situations of Mehraban sheep will be able to create a particular conformation for protein structure and synergistic interactions with TGFβ superfamily such as BMP families.

CONCLUSION

Missense mutation is a kind of point mutation that changes a codon to indicate a different amino acid. This may or may not affect protein function, depending on whether the change is conservative or non-conservative, and what the amino acid actually does. According to sequencing results of Mehraban flock, point mutations of GDF9 are surely in complete linkage disequilibrium. Although neither function wasn’t found for identified point mutations of GDF9 in Mehraban sheep, but we illustrated those four individuals are haplotype A (simultaneous mutations of G1, G2, G3 and G4), one animal was haplotype B (simultaneous mutation of G5 and G6) and seven of them were at wild haplotype (without mutation). Thus we deciphered a perfect correlation between G1/p.Arg87His and G4/p.Glu241Lys in Mehraban sheep. It was verified that G7/p.Val371Met and G8/p.Ser315Phe are probably damaging. However, we haven’t seen any variants of G7 and G8 in Mehraban sheep. Protein modeling of GDF9 at haplotype A and haplotype B individuals had revealed a synergism interaction with TGFβ superfamily. Identified missenses with impressed on turns and helixes of conservative regions of GDF9, caused a synergism to create effects on functions at molecular level.

ACKNOWLEDGEMENT

We would greatly appreciate from Dr Stéphane Fabre, Ph D –HDR, I.N.R.A Centre de recherche de Toulouse, France due to technical assistance, provided reagents or equipment and sequencing of DNA samples without any charges reception. As well as we thank Julien Sarry and Florent Woloszyn for their assistances in the lab.