Bone morphogenetic protein 15 (BMP15) and growth and differentiation factor 9 (GDF9) are TGFbeta-like oocyte-derived growth factors involved in ovarian folliculogenesis as critical regulators of many granulosa cell processes and ovulation rate. Ovarian phenotypic effect caused by alterations in BMP15 and GDF9 genes appears to differ between species and may be relevant to their mono- or polyovulating status. Through phylogenetic analysis we recently showed that these two paralogous genes are strongly divergent and in rapid evolution as compared to other members of the TGFbeta superfamily. Here, we evaluate the amino acid substitution rates of a set of proteins implicated in the ovarian function, including BMP15 and GDF9, with special attention to the mono- or polyovulating status of the species. Among a panel of mono- and polyovulating mammals, we demonstrate a better conservation of some areas in BMP15 and GDF9 within mono-ovulating species. Homology modeling of BMP15 and GDF9 homodimer and heterodimer 3-D structures was suggestive that these areas may be involved in dimer formation and stability. A phylogenetic study of BMP15/GDF9-related proteins reveals that these two genes diverged from the same ancestral gene along with BMP3 and GDF10, two other paralogous genes. A substitution rate analysis based on this phylogenetic tree leads to the hypothesis of an acquisition of BMP15/GDF9-specific functions in ovarian folliculogenesis in mammals. We propose that high variations observed in specific areas of BMP15 and GDF9 in polyovulating species change the equilibrium between homodimers and heterodimers, modifying the biological activity and thus allowing polyovulation to occur.
INTRODUCTION
In female mammals, the biological mechanisms controlling the number of follicles that ovulate at each estrus cycle, i.e., the ovulation rate and to some extent litter size, are still puzzling. Women, cows, and ewes usually have one or two offspring and are considered as mono-ovulatory species, whereas other mammals, such as rodents, dogs, cats, or sows, are highly prolific and produce four or more offspring and are thus considered as polyovulatory species. Numerous studies in sheep have indicated that the ovulation rate and litter size can be genetically regulated by the action of single major genes, named fecundity (Fec) genes [1]. Among these Fec genes, bone morphogenetic protein 15 (BMP15) and growth and differentiation factor 9 (GDF9) are paralogous genes encoding the most important oocyte-derived factors regulating ovarian folliculogenesis [2]. Interestingly, numerous mutations have been found in both genes that are associated with an alteration in ovarian function. In sheep, eight mutations in BMP15 and four mutations in GDF9 are associated with hyperprolificacy or sterility [345678910–11]. In human, numerous missense mutations in BMP15 and in GDF9 have been described associated with premature ovarian insufficiency (POI), dizygotic twinning, or ovarian hyperstimulation syndrome [121314151617181920–21]. In mice, homozygous Bmp15-null females are fertile and exhibit only slight ovarian defects [22], whereas a homozygous null mutation of Gdf9 leads to sterility with a blockade at the primary stage of folliculogenesis resembling the phenotype of BMP15 homozygous sterile ewes [3, 5]. However, heterozygous null mice for Bmp15 or Gdf9 never exhibit increased ovulation rate as observed in BMP15- or GDF9-heterozygous altered sheep. Thus, the roles of BMP15 and GDF9 appear to differ between species and may be relevant to their mono- or polyovulatory status [23, 24].
BMP15 and GDF9 are members of the TGFβ superfamily of cytokines, and, like many other members of this family, they are synthesized as preproproteins. After signal peptide cleavage, dimeric proproteins are processed by protease of the furin family at a proteolytic cleavage site with a conserved RXXR motif [252627–28]. The processing step results in a noncovalent association between the C-terminal mature protein and the N-terminal propeptide. The propeptide or prodomain region appears to be of importance in the BMP15 and GDF9 biological activity, as shown by the deleterious effect of some propeptide mutations or chimeric proprotein in in vitro test [2930–31]. GDF9 and BMP15 are part of the few TGFβ members with only six cysteine residues in their TGFβ domain, or cystine-knot domain, lacking the cysteine usually involved in the covalent dimer formation, but are still able to homodimerize [14, 30]. It has also been shown that BMP15 and GDF9 can heterodimerize [32], cooperate, and synergize to control granulosa cell function [3334–35].
Regarding all above observations, the central hypothesis being examined in this article is to determine if BMP15 and GDF9 have specialized during evolution to be considered as “key genes” controlling ovulation rate in mammals. We recently provided evidence that BMP15 and GDF9 are strongly divergent and in rapid evolution compared to other members of the TGFβ superfamily [36]. To directly compare BMP15 and GDF9 sequence divergences regarding to species ovulating status, we undertook a phylogenetic analysis of these two proteins in mammals. We examined their amino acid sequences to determine if some areas show differential substitution rates, possibly due to specific selective pressure in mono- or polyovulating species. Because BMP15 and GDF9 are paralogous genes, we also determined their evolutionary origins.
MATERIALS AND METHODS
Sequences
Only metazoan sequences were considered for this study. Orthologous BMP15 and GDF9 proteins were identified from all genomic and expressed sequence tag sequences available in general databases such as NCBI and ENSEMBL for most organisms. We also explored specialized databases, as for the amphioxus or elephant shark, using BLASTP, TBLASTN with default parameters (an e-value cutoff at 0.01 was used in all BLAST searches). Accession numbers for each gene are shown in Supplemental Tables S1 to S13 (Supplemental Data are available online at www.biolreprod.org).
Ancestral Character Reconstruction
Ancestral character reconstruction was made by parsimony method using Mesquite software [37] with a user tree of life and the litter size information found in PanTHERIA database [38] as a continuous character. When a species mean litter size was not available, the value for a closely related species was considered. We arbitrarily placed the litter size cutoff value between mono- and polyovulating species at 1.6.
Phylogenetic Analysis
The alignments of amino acid sequences were made using Clustal Omega [39, 40] for the phylogenetic analysis and the substitution rates on the deuterostomes tree (Supplemental Tables S12 and S13) or PRANK software [41] for the substitution rate analysis in mammals (Supplemental Tables S1 to S11).
The evolutionary history was inferred by using the maximum likelihood method. The model of protein evolution that best fit our set of aligned sequences was chosen according to the calculations provided by ProtTest 3.2.1 [42] and MEGA6.06 [43]. Thus, the analysis was based on the JTT matrix-based model [44]. The tree with the highest log likelihood (−16 251.0253) is shown. Branch width was adjusted according to the percentage of trees (out of 500 replicates) in which the associated taxa clustered together (for the branch upstream of the node). Initial tree(s) for the heuristic search were obtained by applying the neighbor-joining method to a matrix of pairwise distances estimated using a JTT model. A discrete gamma distribution was used to model evolutionary rate differences among sites (five categories [+G, parameter = 1.0990]). The rate variation model allowed for some sites to be evolutionarily invariable ([+I], 5.6134% sites). The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. The analysis involved 51 amino acid sequences. All positions with less than 90% site coverage were eliminated. That is, fewer than 10% alignment gaps, missing data, and ambiguous bases were allowed at any position. There were a total of 236 positions in the final dataset. Evolutionary analyses were conducted in MEGA6.06. [43].
Synteny Analysis and Paralogon Detection
Synteny and orthology between human and invertebrate genes were assessed by using the Synteny Database Web server [45] with a window of 50, 100 or 200 genes; Ciona intestinalis was taken as an outgroup for the paralogon research.
Substitution Rates
To evaluate the substitution number per site for the 24 mammal species and from Branchiostoma floridae to eutherians, we used the methodology described in Petit et al. [46] and Monestier et al. [47]. Briefly, a user tree of life (Fig. 1) was used for running the parsimony program Protpars included in the PHYLIP Package [48].
Statistical Analysis of Substitution Rates in Mono-Ovulating and Polyovulating Species
We first estimated for each gene the number of nucleotide nonsynonymous mutations at each amino acid position of the protein sequences within each group of ovulation type: 13 mono-ovulating and 11 polyovulating species. Differences calculated from these raw numbers can be affected by 1) the particular sample of species obtained from grouping by ovulating type and 2) the patterns of substitution at a particular position. Also, substitution rates exhibit large variations even at neighboring sites. To correct for these effects, we transformed raw substitution numbers at each site as explained below.
First, to account for local variation in the number of substitutions, we used smoothed local estimates based on a rolling average over 20 amino acids.
Second, to assess whether the differences in substitution rates could be explained by the ovulation type rather than random sampling of species, we performed 10 000 independent samplings of either 11 (for the mono-ovulating group) or 13 (for the polyovulating group) species out of the 24 total. For each permutation, we estimated the number of substitutions in the resulting phylogenetic tree and computed smoothed local estimates of substitution rate for each permutation and position as described above.
The manual comparison of the distribution of substitution rate obtained by permutation to the substitution rate of our mono- and polyovulating group for each protein permitted us to define, for only BMP15 and GDF9 proteins, “large differential segments” more conserved in mono-ovulating species.
Calculation of standardized substitution rates within ovulation type and global differences between mono- and polyovulating species.
By performing permutations, we were able to sample from the possible distributions of the substitution rates, which we found to be close to a normal distribution at a given position in most cases. Only for positions corresponding to highly conserved regions across all 24 species was this not true. These regions were thus discarded for the rest of the analysis. For other positions, we obtained for each gene, at and for each group, a standardized substitution rate for type (mono- or polyovulating type) and gene at each position. We observed that for all genes considered the substitution rate in mono-ovulating species was globally lower than for polyovulating species (see Results for a discussion of the possible explanations), with polyovulating species rates being slightly higher, and mono-ovulating species strongly lower than expected (Supplemental Fig. S1). To correct for this effect when examining gene regions, we estimated a background substitution rate for each ovulating type and subtracted it.
Differences in substitution rates between ovulating type across protein sequences.
As our final substitution rate estimates were expected to follow a standard normal distribution, the difference between ovulating types was expected to follow a normal distribution of mean 0 and variance 2. We used this distribution to test whether the null hypothesis could be rejected at each position for each gene. We considered as potentially interesting regions where the P value for the null hypothesis was less than 0.01.
We detected three genes harboring regions clearly more conserved in the mono-ovulating species group: BMP15, GDF9, and KIT (see Supplemental Fig. S2A, Fig. 2, Results, and Discussion).
FIG. 2
Substitution rate within various proteins. For each site the substitution rate represented was calculated from a moving average using a window size of 20 amino acids. The site position from start methionine is on the abscissa and the mean substitution rate on the ordinate. The position of the furin cleavage site of BMP15 and GDF9 is represented by a red dashed line. The substitution rate for all 24 species is represented by a green dashed curve. The substitution rate of the 11 polyovulating species is represented in red and the substitution rate of the 13 mono-ovulating species in blue. Dashed box represents large differential segments between mono- and polyovulating species, solid gray box highlights areas significantly more conserved in mono-ovulating species with P value < 0.01 after normalization of the substitution rates.
We also noticed three genes with segments more conserved in the polyovulating species: BMP3, gastric inhibitory polypeptide (GIP), and KIT (see Supplemental Fig. S2B). Going back to the local substitution rates in these regions, we found that they corresponded to regions where mono-ovulating species has substitution rates that were very similar to those of polyovulating species and consistent with the expectation from permutations. They produced significant P values because this is not the average behavior for mono-ovulating genes, but their interpretation is not clear.
Protein Structure Prediction
The 3-D models of GDF9 and BMP15 for the different species were built by homology with the Protein Data Bank TGFB1 3RJR structure [49]. Computation was made by Modeller software [50]. Interaction energies for the different dimers were computed using PISA [51].
Molecular graphics and analyses were performed with the PyMOL Molecular Graphics System, Version 1.5.0.4 (Schrödinger, LLC).
RESULTS
Ancestral Character Reconstruction
For the present study we selected 24 eutherian mammals with complete sequence information for BMP15 and GDF9 in each species. The mono- or polyovulatory status of the studied species was checked using the PanTHERIA database information on litter size [38]. We arbitrarily placed the litter size cutoff value at 1.6, placing the giant panda (mean litter size 1.62) in the polyovulating species. This is consistent with triplet observation [52] and mean litter size greater than two (2.14) observed in some giant panda populations [53]. The mono- or polyovulatory status of each species is presented in Figure 1 (13 mono-ovulating vs. 11 polyovulating) along with the reconstruction of ancestral litter size character by parsimony method using the Mesquite software [37]. Even if the ancestral character was predicted to be polyovulatory, the mono- or polyovulating character appeared several times independently in the tree, giving evidence for parallel evolution of this character in mammals (Fig. 1).
Variation in Substitution Rates Between Mono- and Polyovulating Species
To check if specific selective pressure in genes controlling prolificacy, particularly BMP15 and GDF9, could be linked to ovulation rate, we looked at the substitution rates along the amino acid sequence of these proteins for the 24 mammalian species, using parsimony method based on Protpars software included in the Phylip package [47, 48]. We compared the BMP15 and GDF9 substitution rate results to various members of the TGFβ family implicated (BMP5, BMP2, or inhibin alpha [INHA]), or not (BMP3, GDF10, or myostatin [MSTN]) in the ovarian function. We also checked for KIT, which is involved in ovarian folliculogenesis but is not a member of the TGFβ superfamily, showing an elevated substitution rate among vertebrates [54]. We also selected proteins neither members of the TGFβ superfamily nor implicated in folliculogenesis, such as interleukin 6 and GIP (Fig. 2).
The results demonstrated a great conservation of the C-terminal cystine-knot mature domain of each TGFβ member, the prodomain being more variable. Among all proteins tested, we found a major difference between the mono- and polyovulating groups only for the conservation of the sequences of BMP15 and GDF9 proteins. Indeed, some segments in BMP15 and GDF9 propeptide (two segments in BMP15 and three in GDF9) were more conserved in the mono-ovulating group as compared to the polyovulating group (large differential segments). These segments were situated at positions 154–193 and 240–292 relative to the human BMP15 amino acid sequence, and 1–59, 163–223, and 283–343 relative to the human GDF9 amino acid sequence (Fig. 2 and Supplemental Figs. S3 and S4). The BMP15 240–292 segment and the GDF9 283–343 segment both included the consensus cleavage site situated between the pro and the mature domains.
We also observed that mono-ovulating species had lower global substitution rates than the polyovulating animals in all the proteins we tested. This was probably due to the generation time in our 13 mono-ovulating species (mean 2242 days, median 1460 days; PanTHERIA database and Zhu et al. [52]) that was four times higher than the value observed in the 11 polyovulating animals studied (mean 537 days, median 350 days), resulting in less meiotic events and thus fewer mutation fixations in these species. In order to account for this global difference in substitution rates, we looked for gene regions exhibiting significant differences in substitution rates between mono- and polyovulating species (see Materials and Methods). This analysis revealed that the five segments previously highlighted included short areas at positions 155–163 and 248–264 in BMP15 and 1–34, 193–211, 288–296, and 322–337 in GDF9 significantly more conserved in mono-ovulating species, with P value < 0.01. We also noticed that two short zones including three and nine amino acids in KIT were also significantly more conserved in the mono-ovulating species. Altogether, these data provided evidence of a selective pressure on BMP15 and GDF9 resulting in lower sequence variability specifically in mono-ovulating species.
Homology Modeling of Homodimers and Heterodimers
To determine if the differentially conserved segments are important for structure/function of BMP15 and GDF9, we built homology models of the homodimer and heterodimer 3-D structures of the human and mouse BMP15 and GDF9 proteins (Fig. 3 and Supplemental Fig. S5). Human and mouse were chosen as prototypic examples of mono- and polyovulatory species and also because of an abundant literature on BMP15 and GDF9 biological action in these species. Because neither the BMP15 nor the GDF9 3-D structures were known, they were modeled by homology from the TGFB1 homodimer structure [49] using the Modeller software [50]. Because sequence identities were very low (around 20%), the sequences of BMP15 and GDF9 of both species were manually aligned with that of the pig TGFB1 (PDB 3RJR). Models were built separately for BMP15 homodimer, GDF9 homodimer, and BMP15-GDF9 heterodimer in both species. The analysis of these models showed that the differentially conserved segments between mono-ovulating and polyovulating species are localized in regions at the interface between monomers (Fig. 3). To evaluate the effect of the observed sequence variations on the stability of the different dimers, we computed the interaction energies using the PISA software [51] (Supplemental Table S14). Results show that in human the heterodimer is predicted to be more stable than the homodimers, whereas in mouse the most stable species is the GDF9 homodimer.
FIG. 3
Three-dimensional structure of BMP15 and GDF9 human heterodimer. Homology model of the three-dimensional structure of human BMP15 and GDF9 monomers (A) and the BMP15/GDF9 heterodimer (B) computed by Modeller using TGFB1 structure 3RJR [49] as template. Large differential segments between mono- and polyovulating species are represented in blue for BMP15 and orange for GDF9. Segments significantly more conserved in mono-ovulating species (P < 0.01) are represented in dark violet for BMP15 and red for GDF9. In A, mature domain is represented in dark gray for BMP15 and dark green for GDF9 and brackets show the position of the β4 and β5 strands.
Evolutionary Origin of BMP15 and GDF9
A global search for related BMP15 and GDF9 protein sequences by BLAST in public databases (Ensembl, NCBI, JGI) revealed that these genes were present in all groups of vertebrates. To find the common ancestor, we used elephant shark sequences to search for BMP15/GDF9-related protein in the invertebrate genomes. We identified three amphioxus (B. floridae) proteins: BMP2/4, already known as the common ancestor of BMP2 and BMP4 [55], the antidorsalizing morpho protein (ADMP) and a protein on the scaffold 167 (Bf_V2_167; GenBank accession no. XP_002595999), previously annotated as a BMP3/3b (BMP3b also known as GDF10) protein [56].
To determine the evolutionary origin of BMP15 and GDF9 and their relationship with the proteins retrieved in our BLAST search, a phylogenetic tree was deduced from maximum likelihood method (Fig. 4). The B. floridae BMP2/4 protein was clearly situated in the BMP2/BMP4 clade (blue), along with the BMP2/4 of the sea urchin (Strongylocentrotus purpuratus) and acorn worm (Saccoglossus kovalevskii) in accordance with their status of common ancestors of the BMP2 and BMP4 proteins. Regarding the B. floridae ADMP protein, it was situated, along with the S. kovalevskii ADMP protein, outside of the BMP15/GDF9, BMP3/GDF10, and BMP2/BMP4 clades. The protein present in scaffold 167 of B. floridae was suggested to belong to the BMP15/GDF9 clade (red) but was branched with a poor bootstrap value (56%). In this tree, the BMP3/GDF10 (yellow) and the BMP15/GDF9 (red) clades diverged from a common node supported by a bootstrap value of 73%. In order to confirm the common origin of BMP15/GDF9 and BMP3/GDF10, we searched for blocks of paralogs, i.e., paralogons, in the human genome using the Synteny Database [45]. Ten paralogous genes supported the duplication of GDF9 and BMP3 from a common ancestor with a sliding window of 50 genes (Fig. 5A). In addition, genes surrounding human BMP3 (15 genes on HSA4), GDF10 (17 genes on HSA10) and GDF9 (five genes on HSA5) had orthologous genes in the Bf_V2_167 scaffold of B. floridae (Fig. 5, B and C). Concerning BMP15, no synteny blocks were retrieved from the Synteny Database. Only the SHROOM4 gene (Ensembl accession no. ENST00000483955) lying close to BMP15 on HSAX could form a very small paralogon with SHROOM1 (Ensembl accession no. ENST00000488072) and GDF9 on HSA5 and with SHROOM3 (Ensembl accession no. ENST00000484236) and BMP3 on HSA4.
FIG. 4
Evolution of BMP15- and GDF9-related proteins. The evolutionary history of BMP15- and GDF9-related proteins was inferred by using the maximum likelihood method based on the JTT matrix-based model. Branch width was adjusted according to the bootstrap value (for the branch upstream of the node) calculated out of 500 replicates. A discrete gamma distribution was used to model evolutionary rate differences among sites. The rate variation model allowed for some sites to be evolutionarily invariable. The tree is drawn to scale, with branch lengths measured in the number of substitutions. BMP3/GDF10 (yellow) and BMP15/GDF9 (red) have a common origin. The protein present in scaffold 167 of B. floridae was suggested to belong to the BMP15/GDF9 branch but was supported by a poor bootstrap value (56%).
FIG. 5
Synteny and paralogy of B. floridae and human BMP15-related genes. A) Paralogy between human chromosome Hsa4 and Hsa5. Ten paralogous genes shared the same duplication history as GDF9 and BMP3 with a sliding window of 50 genes. B) Orthology between chromosome Hsa5 and the scaffold Bf_V2_167 of B. floridae. C) Orthology between chromosome Hsa10 and Hsa4 and the scaffold Bf_V2_167 of B. floridae.
We were able to conclude that the BMP15/GDF9 and the BMP3/GDF10 clades diverged from a common ancestral gene encoding a protein probably closely related to the B. floridae ADMP and the scaffold 167 proteins.
Substitution Rate Between Phylogenetic Tree Nodes
To determine the evolution of the selection pressure after the BMP3/GDF10 and BMP15/GDF9 common ancestor divergence, we checked for the evolution of the substitution rate in the branch of a phylogenetic tree, still using the parsimony method from Protpars (Fig. 6). The comparison of the substitution rates between the nodes from S. kowalevskii and B. floridae ADMP to eutherian BMP15 and GDF9 showed a progressive decrease in the substitution rate, particularly in the TGFβ domain (red box in Fig. 6), from the BMP15/GDF9 ancestral protein (node 2 in Fig. 6). Two other areas corresponding to the amino acid segments 34–56 and 66–85 of human BMP15 and 49–72 and 85–104 of human GDF9 (segments A and B in Fig. 6) appeared to be also progressively more conserved during evolution. Specifically for GDF9 sequences, we observed a highly variable zone (segment C in Fig. 6) from the BMP15/GDF9 duplication onward (nodes 3–7′). This area situated between amino acid positions 280 and 308 in human GDF9 was part of a differential segment associated to mono- or polyovulation highlighted in Figure 2. The relative stability of the substitution rate observed after the amphibian divergence (nodes 5 and 5′) might indicate that the BMP15 and GDF9 protein functions described in eutherians were already acquired in amniotes.
FIG. 6
Substitution rates of BMP15/GDF9 proteins in the different branches of the evolutionary tree. For each site the substitution rate was calculated from a moving average using a window size of 20 amino acids. The site position from start methionine is on the ordinate and the mean substitution rate on the abscissa. Each curve represents the substitution rates in the phylogenetic tree branch of the same color. An artifact due to a gap in the alignment, caused by a BMP3/GDF10-specific segment, appears from node 3. The position of the furin cleavage site is represented by a red dashed line. TGFβ domain (red box) and segments A (dark blue box) and B (violet box) show progressive stabilization during evolution. Segment C (light blue box) is highly variable in the branch between nodes 3 and 7′.
DISCUSSION
It was previously shown that BMP15 and GDF9 evolved more quickly than the other members of the TGFβ family, with evidence for positive selection in BMP15, especially in Hominidae [36]. Ancestral character reconstruction based on actual mammal litter size information on a panel of 24 species including 13 mono-ovulating and 11 polyovulating species showed parallel evolution leading to mono- or polyovulation status. We used a parsimony algorithm to obtain the substitution rate of each amino acid to investigate if the ovulatory status was linked to variations in BMP15 and GDF9 sequences, each being known to regulate ovulation rate in mammals. Indeed, our study highlighted specific segments significantly more conserved in BMP15 and GDF9 from mono-ovulating as compared to polyovulating species. Except for a weak signal in the KIT protein, we did not find any difference linked to ovulation status in the other proteins we tested, indicating a specific evolutionary pressure for this phenotype on BMP15 and GDF9. Nevertheless, the differential signal obtained in KIT may also be relevant to selective pressure on the ovarian function because it participates with BMP15 in the intraovarian signaling cross-talk between oocyte and somatic cells [57]. The approach of using substitution rate to map between genotype and phenotype was also successfully applied for gene evolution linked to male reproduction, when Semenogelin 2 was found to be significantly associated to mating systems in primates [58].
In the BMP15 and GDF9 primary sequences, the fragments differentially affected by the ovulatory status are mainly located in the prodomains. Among the eight mutations in BMP15 affecting ovulation rate in sheep, three are situated in the prodomain and locate exactly in these differentially selected segments (W155, K241, and Q290, related to human BMP15 sequence). In human, among the 17 mutations affecting BMP15 function associated with POI, 16 are also located in the prodomain [59], but only 3 are situated in the differentially selected segments (A180, I243, and L262). For GDF9, all mutations found in sheep are in the mature domain and do not relate with the selected segments we found in this study. In human GDF9, 2 mutations (among 12) associated with POI locate in a selected segment (S186 and V216). All these natural mutations present in sheep and human impaired protein function and thus are located in amino acids of great importance for protein activity (folding, cleavage, etc.). It is not surprising that important amino acids for BMP15 and GDF9 protein activity are also located outside the differentially conserved segments we highlighted. Recent studies have identified species-specific amino acid positions implicated in the difference of BMP15 or GDF9 expression and activity between mouse and human [24, 60]. These residues are located essentially in the mature domain and not in the segments we highlighted. We did not identify these regions essentially for two reasons: 1) this type of lineage-specific variation (mouse vs. human) led to no major difference in our substitution rates between numerous poly- and mono-ovulating species and 2) a single amino acid variation situated in a well-conserved area could not be detected in our analysis using a sliding window of 20 amino acids.
Interestingly, for both BMP15 and GDF9, one of the differentially selected segments included the consensus furin cleavage site separating the prodomain and the mature domain. This cleavage site is of great importance for the processing of the mature proteins and thus their cellular signaling. Particularly for BMP15, defects in the in vitro production of the mouse protein were observed, whereas the human BMP15 was perfectly produced [29]. However, when combined with the cleavage site and the prodomain of human BMP15, mouse BMP15 mature domain could be processed and express its biological activity [23]. Thus, we hypothesize that a high selective pressure in the environment of the furin cleavage site could lead to keeping an efficient processing of the mature domain and then the mono-ovulation status. In contrast, a more variable segment around the furin cleavage site that is characteristic of polyovulating species may lead to a less efficient processing and a lower biological activity of BMP15 and GDF9, driving the polyovulation phenotype.
Regarding the position of the differentially selected segment in the 3-D structure prediction of BMP15 and GDF9 homo- and heterodimers, it is also suspected that variations within these segments could strongly affect the interaction between monomers. Indeed, one of the more conserved segments in BMP15 and GDF9 (positions 154–193 in human BMP15 and 163–223 in human GDF9) from mono-ovulating species is situated in the interface between the two propeptides corresponding to the β4 and β5 strands of TGFB1 (brackets in Fig. 3A). These two β-strands conserved in all other TGFβ members except in GDF1 and GDF15 seem important for the interaction between propeptides and thus for dimerization [49]. The absence of these β-strands in GDF1 and GDF15 could be due to specific adaptation of these proteins for some function for which propeptide interaction between monomers is not obvious [6162–63]. The implication of the segments we highlighted in the interaction between monomers is also assessed by stability estimations obtained through the PISA software applied to the 3-D models of the different dimers on different species. Because 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. However, the models and the stability estimations tend to show that the observed sequence variations could induce a difference in the equilibrium between the different types of dimers, favoring the heterodimer in mono-ovulating species and the homodimers in polyovulating species. This tendency is in accordance with functional data previously published indicating a greater biological activity of BMP15-GDF9 heterodimer as compared to GDF9 homodimer in human [35]. Also in mouse, the BMP15 homodimer is barely active [23] whereas GDF9 homodimer and BMP15-GDF9 heterodimer have high activities [35]. Altogether, these results show that mono- and polyovulating species present a difference of selective pressure specifically in some segments of BMP15 and GDF9, two proteins implicated in the control of ovulation rate in the ovine species. We propose that the variations observed in these segments in polyovulating species alter both the processing and the dimerization of BMP15/GDF9, resulting in a decrease of the proteins' biological activity. Based on the physiological and molecular hypotheses we developed for the ovulation rate control in sheep [64], we assume that the decreased biological activity of BMP15 and/or GDF9 in polyovulating species affected the proliferation and the gonadotropin-dependent action on granulosa cells, as observed in BMP15 and/or GDF9 mutated hyperprolific ewes. Differences in the proteins' biological activity between mono- and polyovulating species could also be due to variations in ligand-receptor interactions, but the areas known to be implicated in those interactions are located in the mature domain [49] and are not included in the segments we highlight. At the molecular level, the mono-ovulation quota would be tightly controlled by the integrative biological action of BMP15 homodimer, BMP15/GDF9 heterodimer, and GDF9 homodimer. Alteration in the equilibrium between the three pathways due to species-specific sequences, to punctual mutations, or to change in expression (or receptivity) could lead to polyovulation. Whether the species-specific molecular consequence on granulosa cell proliferation or steroidogenesis passes through SMAD-dependent or SMAD-independent pathways needs to be fully addressed in the near future.
The present analysis considers only the evolutionary impact on the protein sequences of BMP15 and GDF9. This type of “static” evolution is time dependent and leads to adaption of the function of the protein due to mutation accumulation. However, “dynamic” processes with evolutionary impact on the regulatory sequences of the two genes may also be relevant in the explanation of the major role of BMP15 and GDF9 in the control of ovulation rate. The BMP15/GDF9 homodimers and heterodimer equilibrium could be influenced by the species-specific ratio of BMP15/GDF9 mRNA expression. Indeed, the BMP15/GDF9 expression ratio was shown to be different regarding the ovulatory status of ruminants versus rodents [65]. Thus, it is reasonable to think that the global effect of BMP15 and GDF9 genes on ovulation rate is due to a species-specific combination of static and dynamic evolution processes.
Looking at the evolution of BMP15 and GDF9 in deuterostomes, the topology of our maximum likelihood tree is in accordance with the phylogeny of Adoutte et al. [66] and Cameron et al. [67]. However, some artifacts are still present, as shown by the position of BMP2/4 of S. purpuratus, S. kovalevskii, and B. floridae in a same clade. These artifacts are probably due to long-branch attraction. We determined a common origin of BMP3/GDF10 and BMP15/GDF9. The existence of a common ancestor of these BMPs is evidenced by the presence of synteny blocks containing these genes. We can easily hypothesize the duplication of BMP2/4 in BMP4 and BMP2 occurring on the second round of whole-genome duplication [68, 69] because of the presence of the ancestral BMP2/4 gene in Petromyson marinus. In contrast, because of the lack of information present in the P. marinus taxa, it is difficult to evaluate when the duplication leading to BMP15/GDF9 on the one hand and BMP3/GDF10 on the other hand has occurred. Nevertheless, the elephant shark genome contains BMP15 and GDF9, indicating the presence of the two genes in Chondrichthyes and duplication before the vertebrate's common ancestor. The B. floridae protein encoded in scaffold 167 was previously annotated BMP3/3b [56]. However, based on our phylogenetic analysis, the absence of the fourth cysteine in the TGFβ domain, and the gene exon-intron structure, we assume that this protein is more likely related to BMP15 and GDF9. The common ancestor of BMP3/GDF10 and BMP15/GDF9 is probably closely related to ADMP found in S. kovalevskii and B. floridae and to the BMP15-related protein found in the scaffold 167 of B. floridae genome.
The substitution rate evolution from the predicted ancestor of BMP3/GDF10 and BMP15/GDF9 to eutherian BMP15/GDF9 demonstrates a progressive stabilization of the variations in the TGFβ cystine-knot domain carrying the biological activity of the two proteins. We also observed two other zones that seem to be progressively more conserved during evolution and are situated in the α1 (segment A in Fig. 6) and α2 (segment B in Fig. 6) helix of TGFB1 structure (conserved in our BMP15 and GDF9 structure) described as being particularly important to template the folding of TGFβ [49]. Interestingly, the first of these segments (segment A in Fig. 6) located between the amino acids 34 and 56 on human BMP15 and 49 and 72 on human GDF9 includes the five amino acids (positions 44–50) shown to mediate the BMP15 mature protein expression in human, mouse, and sheep [24]. Based on these observations, one might suggest that the actual roles of BMP15 and GDF9 seem to be acquired in amniotes. Specifically for GDF9, we observed the appearance of a highly variable area in the branch leading to vertebrates (segment C in Fig. 6). This area is part of one segment highlighted as more conserved in mono-ovulating mammals as compared to polyovulating ones. We can hypothesize a relaxation of the selective pressure occurring for this segment during divergence of GDF9 after the ancestral BMP15/GDF9 gene duplication, but we do not have any functional interpretation of this high variability.
In conclusion, we have shown that BMP15 and GDF9 evolved from the same ancestral gene along with BMP3 and GDF10. The functions of the proteins are probably already acquired in amniotes. The role of BMP15 and GDF9 in controlling ovulation rate has been formally proven in the ovine species. Here, through phylogenetic analysis we reveal data indicating that these two genes are affected by differential selective pressure leading to modifications of their biological activity. We hypothesize that this study partly explains the variability of the ovulation status observed among mammals.
REFERENCES
Notes
[1] Financial disclosure This work was sponsored by grants from the French “Agence Nationale de la Recherche” (ANR 2010 BLAN 1608 01, MONOPOLY) and the French “Centre National de la Recherche Scientifique” (CNRS PICS program). Presented in part at the 17th Evolutionary Biology Meeting, 17–20 September 2013, Marseille, France.
