# Genome-Wide Identification of a Regulatory Mutation in BMP15 Controlling Prolificacy in Sheep

Louise Chantepie, Loys L. Bodin, Julien Sarry, Florent Woloszyn, Florence

Plisson-Petit, Julien Ruesche, Laurence Drouilhet, Stéphane Fabre

► To cite this version:

Louise Chantepie, Loys L. Bodin, Julien Sarry, Florent Woloszyn, Florence Plisson-Petit, et al. Genome-Wide Identification of a Regulatory Mutation in BMP15 Controlling Prolificacy in Sheep. Frontiers in Genetics, 2020, 11, 10.3389/fgene.2020.00585. [hal-02909721](https://hal.inrae.fr/hal-02909721)

HAL Id: hal-02909721

https://hal.inrae.fr/hal-02909721

Submitted on 31 Jul 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

# Genome-Wide Identification of a Regulatory Mutation in BMP15 Controlling Prolificacy in Sheep

Louise Chantepie†, Loys Bodin, Julien Sarry, Florent Woloszyn, Florence Plisson-Petit, Julien Ruesche, Laurence Drouilhet and Stéphane Fabre*

GenPhySE, Université de Toulouse, INRAE, ENVT, Castanet-Tolosan, France

The search for the genetic determinism of prolificacy variability in sheep has evidenced several major mutations in genes playing a crucial role in the control of ovulation rate. In the Noire du Velay (NV) sheep population, a recent genetic study has evidenced the segregation of such a mutation named *FecL*^L^. However, based on litter size (LS) records of *FecL*^L^ non-carrier ewes, the segregation of a second prolificacy major mutation was suspected in this population. In order to identify this mutation, we have combined a case/control genome-wide association study with ovine 50k SNP chip genotyping, whole genome sequencing, and functional analyses. A new single nucleotide polymorphism (OARX:50977717T > A, NC_019484) located on the X chromosome upstream of the *BMP15* gene was evidenced to be highly associated with the prolificacy variability (*P* = 1.93E-11). The variant allele was called *FecX*^N^ and shown to segregate also in the Blanche du Massif Central (BMC) sheep population. In both NV and BMC, the *FecX*^N^ allele frequency was estimated close to 0.10, and its effect on LS was estimated at +0.20 lamb per lambing at the heterozygous state. Homozygous *FecX*^N^ carrier ewes were fertile with increased prolificacy in contrast to numerous mutations affecting *BMP15*. At the molecular level, *FecX*^N^ was shown to decrease *BMP15* promoter activity and supposed to impact *BMP15* expression in the oocyte. This regulatory action was proposed as the causal mechanism for the *FecX*^N^ mutation to control ovulation rate and prolificacy in sheep.

**Keywords:** GWAS, major gene, prolificacy, sheep, BMP15

## OPEN ACCESS

**Edited by:**
Luis Varona,
University of Zaragoza, Spain

**Reviewed by:**
Shahin Eghbalsaed,
Islamic Azad University, Iran
Brian Kirkpatrick,
University of Wisconsin–Madison,
United States

*Correspondence:
Stéphane Fabre
stephane.fabre@inrae.fr

†Present address:
Louise Chantepie,
Institut de l'Elevage,
Castanet-Tolosan, France

**Specialty section:**
This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics

**Received:** 16 January 2020
**Accepted:** 14 May 2020
**Published:** 19 June 2020

**Citation:**
Chantepie L, Bodin L, Sarry J, Woloszyn F, Plisson-Petit F, Ruesche J, Drouilhet L and Fabre S (2020) Genome-Wide Identification of a Regulatory Mutation in BMP15 Controlling Prolificacy in Sheep. *Front. Genet.* 11:585. doi: 10.3389/fgene.2020.00585

## INTRODUCTION

There is now an accumulation of evidence that the oocyte plays a central role in controlling the ovarian folliculogenesis from the early stages up to ovulation. Among the local factors produced by the oocyte itself, members of the bone morphogenetic protein/growth and differentiation factor (BMP/GDF) family play an integral role in this control (Persani et al., 2014). Among them, the most important are surely *BMP15* and *GDF9*. Knockout mice models gave the first evidence of the importance of these two oocyte-derived factors acting individually as homodimers and/or through a synergistic cooperation to control the ovarian function (Elvin et al., 1999; Yan et al., 2001).

In humans also, a focus was done on *BMP15* and *GDF9* about their implication in various ovarian dysfunctions. Indeed, numerous heterozygous missense mutations have been identified in both genes associated with primary or secondary amenorrhea in different cohorts of women
affected by primary ovarian insufficiency (POI) all over the world. Particularly, the 10-fold higher prevalence of *BMP15* variants among patients with POI compared with the control population supports the causative role of these mutations (Persani et al., 2014). Alterations of *BMP15* and *GDF9* were also searched in association with polycystic ovary syndrome (PCOS). Here, again, several missense variants were discovered in both genes, but the pathogenic role of these mutations remains controversial in the etiology of this syndrome. However, several studies have reported an aberrant expression of *BMP15* and *GDF9* in the ovary of PCOS patients (Teixeira Filho et al., 2002; Wei et al., 2014). Interestingly, some *BMP15* polymorphisms situated in the 5'UTR are significantly associated with the over-response to recombinant FSH applied during assisted reproductive treatment and with the risk to develop ovarian hyper-stimulation syndrome (OHSS, Morón et al., 2006; Hanevik et al., 2011). Finally, polymorphisms in *BMP15* and *GDF9* genes were also searched in association with dizygotic twinning in humans. If no convincing results were obtained for *BMP15*, some loss-of-function variants of *GDF9* were observed significantly more frequently in mothers of twins compared to the control population (Palmer et al., 2006).

In parallel, the search for the genetic determinism of ovulation rate and prolificacy variability in sheep has also highlighted the crucial role of *BMP15* and *GDF9* by evidencing numerous independent loss-of-function mutations, all altering the coding sequence of these two genes (Persani et al., 2014; Abdoli et al., 2016). Depending on the mutation and its hetero- or homozygous state, the phenotype controlled by these mutations in *BMP15* and *GDF9* goes from the early blockade of folliculogenesis and subsequent sterility to an extraordinary increase of ovulation rate (OR) and, thus, litter size (LS) of carrier ewes (Galloway et al., 2000; Hanrahan et al., 2004; Silva et al., 2011; Demars et al., 2013). Thus, sheep exhibiting an extremely high prolificacy are of great interest for identifying genes and mutations involved in molecular pathways controlling ovarian function. These animal models have a double interest in agriculture for genetic improvement of the prolificacy and in human clinics for providing valuable candidate genes in the genetic determinism of female infertility or subfertility as described above.

The Noire du Velay (NV) population is a French local sheep breed mainly reared in the Haute-Loire and Loire departments. Ewes present naturally out-of-season breeding ability, very good maternal characteristics, and a quite high prolificacy (mean LS = 1.62 lamb per lambing). Large variation in LS has been observed in this breed, and a recent genetic study has evidenced the segregation of an autosomal mutation named *FecgL* controlling this trait (Chantepie et al., 2018). This variant, located in intron 7 of the *B4GALNT2* gene and associated with its ectopic ovarian expression, was originally discovered in the Lacaune meat sheep breed, increasing OR and prolificacy (Drouilhet et al., 2013). For the segregation study, more than 2700 NV ewes with LS records were genotyped at the *FecL* locus (Chantepie et al., 2018). Surprisingly, the distribution of LS and the existence of highly prolific ewes among the *FecgL* non-carriers have suggested the possible segregation of a second prolificacy major mutation in this population as previously observed in the Lacaune breed carrying both *FecgL* and *FecXL* (Bodin et al.,

2007; Drouilhet et al., 2013). In order to validate this hypothesis, after specific genotyping excluding all other known mutations affecting OR and LS and segregating in French sheep populations, we have performed a genome-wide association study (GWAS) based on a case/control design. Completed by the whole-genome sequencing of two specifically chosen animals based on their genotypes, we have identified a new regulatory variant called *FecXN* proposed to affect the oocyte-dependent expression of *BMP15* in association with increased prolificacy in sheep.

# MATERIALS AND METHODS

## Animals

Ewes (*Ovis aries*) from the NV breed (n = 2266) were genotyped on blood DNA at the *FecL* locus as previously described (Chantepie et al., 2018). In order to test the hypothesis of the segregation of a second major mutation controlling LS in this breed, a first set of 80 ewes with at least 5 LS records (mean LS = 1.84, ranging from 1.00 to 3.50) were selected among the *FecL*⁺ homozygous genotype (n = 2151, mean LS = 1.58). Subsequently, for the NV breed, the effect of the *FecXN* mutation on LS was estimated on 2252 ewes, considering the genotype at the *FecL* locus. The presence of the *FecXN* mutation in other breeds was checked on a diversity panel of 725 animals from 26 French sheep breeds (Rochus et al., 2018; **Table 3**). For the BMC population, the effect of the *FecXN* mutation on LS was estimated on 2456 ewes. For gene expression analysis, 10 homozygous ewes at the *FecXN* locus (five carriers and five non-carriers of the *N* allele) were bought from private breeders (six NV and four BMC) and reared at the INRA experimental facility (agreement number: D3142901). All experimental procedures were approved (approval number 01171.02) by the French Ministry of Teaching and Scientific Research and local ethical committee C2EA-115 (Science and Animal Health) in accordance with the European Union Directive 2010/63/EU on the protection of animals used for scientific purposes.

## Biological Samples

All blood sampling from the numerous sheep breeds studied were collected from the jugular vein (5 ml per animal) by the Venoject system with EDTA and directly stored at -20°C for further use. Part of these blood samples (GWAS and diversity panel) was used for extraction of genomic DNA as described (Bodin et al., 2007). All other samples were used for direct genotyping on whole blood without DNA purification (Chantepie et al., 2018).

For ovary collection and oocyte isolation, the estrus cycles of all adult NV and BMC ewes were synchronized with intravaginal sponges impregnated with flugestone acetate (FGA, 30 mg, CEVA) for 14 days. Ovaries were collected at slaughtering during the follicular phase 36 h after FGA sponge removal. Cumulus-oocyte complexes (COC) were immediately recovered from all visible 1-3 mm follicles by aspiration using a 1 ml syringe with a 26 G needle and placed in McCoy's 5A culture medium (Sigma-Aldrich). COC were mechanically dissociated by several pipetting and washing cycles in 150 µl drops of McCoy's 5A medium, and finally, denuded oocytes devoid of granulosa cells were recovered
in 1X PBS. Only intact oocytes with a good homogeneity of the cytoplasm were grouped to obtain two to three pools of five oocytes per animal and stored at -80°C before RNA extraction.

## Genotyping Analyses

The *FecL*^L^ mutation (OAR11:36938224T > A, NC_019468) was genotyped directly on whole blood samples by the KAPA-KASP assay as previously described (Chantepie et al., 2018). As a prerequisite before GWAS, a set of 30 highly prolific *FecL*^+^/*FecL*^+^ ewes was controlled for the absence of other evidenced major mutations affecting sheep prolificacy in French populations. Using the same KAPA-KASP assay, *FecX*^L^ and *FecX*^Gr alleles in BMP15 were genotyped as described (Chantepie et al., 2018). *FecB*^B^ in the exon 7 of the *BMPR1B* gene (OAR6:29382188A > G, NC_019463.1) was genotyped using forced restriction fragment length polymorphism (RFLP) as described by Wilson et al. (2001).

Whole-genome genotyping was performed on 80 ovine genomic DNA samples using the OvineSNP50 Genotyping Beadchip from Illumina according to the manufacturer's protocol at the Laboratoire d'Analyses Geïneitiques pour les espèces Animales (LABOGENA, Jouy en Josas, France)<fnref n="1" />. From the data set, individuals with a call rate < 0.98 were excluded. SNP exclusion thresholds were call frequency < 0.95 and minor allele frequency (MAF) < 0.01 or a significant deviation from Hardy-Weinberg equilibrium (HWE) in the controls (p < 1.10⁻⁶). Non-polymorphic SNP positions and markers with no position on the OARv3.1 reference genome map were also discarded. Finally, from the available design of 54,241 SNPs available on the Illumina OvineSNP50 Beadchip and 80 selected NV ewes, the final data set was reduced to 47,446 SNPs analyzed in 79 individuals.

The *FecX*^N^ mutation (OARX: 50977717T > A, NC_019484) was genotyped by a RFLP analysis using the Mse1 restriction enzyme (New England Biolabs) after a first step of Terra PCR Direct Polymerase Mix amplification (Takara) using 1 µl sample of total blood. The accuracy of the *FecX*^N^ RFLP genotyping was validated by Sanger sequencing on few samples using the same amplification primers.

PCR amplifications were conducted independently for each locus studied on an ABI2400 thermocycler (Applied Biosystems) with the following conditions: 5 min at 94°C, 32 cycles of 30 s at the specific melting temperature, 30 s at 72°C and 30 s at 94°C, followed by 5 min at 72°C. The primers used in this study are listed in **Supplementary Table S2**.

## Whole-Genome Sequencing (WGS) Analysis

DNA sequencing libraries were constructed from 1 µg of genomic DNA using TruSeq DNA PCR-free Library Prep kit (Illumina). Sequencing was run on an Illumina HiSeq 2500 apparatus using a paired-end read length of 2 × 150 pb with the Illumina Reagent Kits as previously described (Demars et al., 2017). WGS was performed at the Genotoul-GeT core facility (INRA Toulouse)<fnref n="2" />. The raw reads of Illumina DNA

sequencing were preprocessed by removing adapter sequences. After quality control, the FastQ files and metadata were submitted to the European Nucleotide Archive (ENA) at EMBL-EBI (accession number PRJEB35553). Reads mapping and variants calling were performed using the local instance of Galaxy<fnref n="3" /> at the Toulouse Midi-Pyrénées bioinformatics platform<fnref n="4" />. The cleaned paired reads were combined and mapped against the ovine genome assembly (Oar_v3.1.86) using BWA-MEM (Galaxy version 0.7.17.1). The resulting BAM files were sorted using Samtools_sort (Galaxy version 1.0.0). Sorted and indexed BAM files were visualized through the Integrative Genome Viewer, IGV software version 2.4.10 (Robinson et al., 2011). A GFF3 annotation file was obtained from Ensembl (Ovis_aries.Oar_V3.1.78). We applied GATK version 3.5-0 to perform SNP and InDel discovery and genotyping across the two samples simultaneously using standard filtering parameters according to GATK Best Practices recommendations (DePristo et al., 2011; Van der Auwera et al., 2013). Variants effect and annotation were realized by SNPEff version 4.1, and filtering of interesting variants was performed using the SNPSift tool.

## In vitro Transcription and Translation of BMP15

The full-length cDNA of ovine BMP15 [1480 bp (-297, + 1183) referring to ATG start codon] with or without the *FecX*^N^ mutation in position -290 was generated from oocyte-derived RNA after a reverse transcription (RT) step (described in RNA extraction and RT paragraph, primers are listed in **Supplementary Table S2**). The resulting PCR products were inserted by TA cloning into pGEM-T Easy plasmid (Promega) possessing T7 and SP6 promoters. The orientation of insertion and exclusion of unexpected PCR-induced mutations were controlled by Sanger sequencing.

In vitro transcription and translation were realized from 500 ng of cDNA pGEM-T construct using the TnT T7 Quick Coupled Transcription/Translation kit (Promega) and Transcend Biotin-Lysine-tRNA following the manufacturer's protocol. Reactions for each construct were run in duplicate in six independent TnT experiments. The resulting BMP15 protein was detected using a Transcend non-radioactive translation detection system with chemiluminescent method (Promega) after reducing SDS-PAGE on a gradient (4-15%) polyacrylamide gel (Promega) and transfer onto nitrocellulose membrane. Chemiluminescent signal was captured by a ChemiDoc MP imaging system and images were analyzed with the Image Lab Software (Bio-Rad).

## BMP15 Promoter Activity

The promoter sequence of the ovine *BMP15* gene was amplified by PCR on genomic DNA from both homozygous *FecX*^+^ and homozygous *FecX*^N^ ewes. Two sizes of fragments were generated for cloning the *BMP15* promoter in front of the luciferase (Luc) reporter gene, a long (lg) form of 732 bp [(-743, -11) bp referring to ATG start codon] and a short (sh) form of 341bp [(-443, -102) bp]. The PCR products were engineered for digestion

<footnote n="1">¹www.labogena.fr</footnote>

<footnote n="2">²get.genotoul.fr</footnote>

<footnote n="3">³https://galaxyproject.org</footnote>

<footnote n="4">⁴http://sigenae-workbench.toulouse.inra.fr</footnote>
using *KpnI* and *HindIII* restriction enzymes (New England Biolabs) and inserted into the pGL4.23 vector (Promega). The four resulting constructs (IgBMP15⁺-Luc, IgBMP15N-Luc, shBMP15⁺-Luc, shBMP15N-Luc) were controlled by Sanger sequencing. Primers used to generate these constructs are listed in **Supplementary Table S2**. Twenty-four hours after seeding (3.10⁴ cells/well, 24 wells plate), CHO (Chinese hamster ovary) cells were transfected using Lipofectamine 3000 (Invitrogen) with 500 ng/well of pGL4.23 constructs either empty or containing BMP15 promoter fragment. Forty-eight hours after transfection, cells were lysed and assayed for luciferase activity (Luciferase reporter assay kit, Promega). Luminescence in relative light units (RLU) was measured by a Glomax microplate reader (Promega). Each construct was assayed in triplicate in six independent transfection experiments.

## RNA Extraction, Reverse Transcription, and Quantitative PCR

Total RNA from pools of five oocytes were extracted using the Nucleospin RNA XS kit according to the manufacturer's protocol (Macherey-Nagel) and including a DNasel treatment. The low quantity of RNA recovered did not allow quantification. So the equivalent of 1.25 oocyte was reverse-transcribed using SuperScript II reverse transcriptase (Invitrogen) and anchored oligo(dT)22 primer (1 µl at 10 µM). Primer design using Beacon designer 8.20 (Premier Biosoft), SYBR green real-time PCR cycling conditions using a QuantStudio 6 Flex Real-Time PCR system (Thermo Fisher Scientific) and amplification efficiency calculation [$E = e^{(-1/slope^c)}$] were as previously described in Talebi et al. (2018a). Primer sequences, amplicon length, and amplification efficiency are listed in **Supplementary Table S2**. RNA transcript abundance was quantified using the ΔCt method with the mean expression of *GAPDH* and *SDHA* as internal references and following the formula $R = [E_{ref}^{Ct_{ref}}/E_{target}^{Ct_{target}}]$. The two reference genes were validated by the Bestkeeper algorithm (Pfaffl et al., 2004).

## Data Analysis

Single-marker association analyses were conducted using a Fisher's exact test, and a Bonferroni correction has been applied to check for significance levels. The chromosome- and genome-wide values have been established as mentioned by Balding (2006). Statistical analyses were done using PLINK1.9 software under a case/control design (Chang et al., 2015). Among the 79 data sets of 47,446 SNPs analyzed, the LS trait was considered as the case when mean LS ≥ 2.18 (n = 39) and control when LS ≤ 1.45 (n = 40). Haplotypic association analysis on the X chromosome were performed using FastPhase software (Scheet and Stephens, 2006). Empirical significance levels were calculated using maximum statistic permutation approach [max (T), n = 1000].

Allele effect on LS was estimated in NV and BMC breeds on data extracted from the French national database for genetic evaluation and research managed by the Institut de l'Elevage (French Livestock Production Institute) and the CTIG (Centre de Traitement de l'Information Génétique, Jouy-en-Josas, France).

Only females born after 2000 were retained (27,754 NV ewes with 122,110 LS records and 110,848 BMC ewes with 461,405 LS records) with their pedigree over five generations. The *FecX*N genotype effect on the subset of 79 case/control animals was assessed by one-way ANOVA, followed by the Newman-Keuls *post hoc* test. For the large animal cohort analyses, the linear mixed models used were as similar as possible to those of the national genetic evaluation system (Poivey et al., 1995). In the present study, the following fixed effects were considered: (i) the genotype at the *FecX* locus; (ii) the month of birth (12 levels); (iii) a physiological status effect combining parity, age at first lambing, rearing mode, and postpartum interval (44 levels); and (iv) a combination of the flock year and season effect. Two random effects were added to the model: a permanent environmental effect and an animal additive genetic effect. Moreover, an additional fixed effect of the reproduction type was considered for the BMC breed for which some hormonal treatments are used each year (87 and 13% for natural and induced estrus, respectively). For the NV breed, because the *FecL*L allele is also segregating in the population (Chantepie et al., 2018), the effect of the genotype at this locus (2252 known and 25,502 unknown genotypes) as well as its interaction with the genotypes at the *FecX* locus were considered. All these models were fitted using the ASReml software (Gilmour et al., 2009).

The comparison between *FecX* alleles for BMP15 protein quantification was analyzed using Student's *t*-test, using Welch's correction for unequal variance. For reporter luciferase assays, differences between constructs were analyzed by one-way ANOVA followed by Newman-Keuls *post hoc* test. QPCR data for *BMP15* and *GDF9* expression in oocytes were analyzed by two-way ANOVA considering genotype and breed effects. *P* < 0.05 was the threshold for statistical significance. All these experimental data are presented as means ± SEM and were analyzed using Prism 6 (GraphPad Software Inc.).

## RESULTS

### Genetic Association Analyses

A first set of genomic DNA from 30 NV ewes without the *FecL*L prolific allele at the *B4GALNT2* locus (LS records ranging from 2.00 to 3.00) was genotyped for already known mutations affecting sheep prolificacy at the three other loci, *BMPR1B*, *GDF9*, and *BMP15*. Using a specific RFLP assay (*BMPR1B*, Wilson et al., 2001) or Sanger sequencing of coding parts (*GDF9* and *BMP15*, Talebi et al., 2018b), none of the known mutations was evidenced (data not shown). Thus, to establish the genetic determinism of the remaining LS variation in this population, 80 ewes were genotyped by an Illumina Ovine SNP50 Genotyping Beadchip. The allele frequencies of the most highly prolific ewes (cases, *n* = 40, mean LS = 2.47) and lowly prolific ewes (controls, *n* = 40, mean LS = 1.23) were compared to identify loci associated with LS using GWAS according to the procedures described in section "Materials and Methods." Finally, genotype data were obtained from 79 animals (39 cases, 40 controls). Six markers located on OARX were significantly associated with LS variation at the genome-wide level after Bonferroni correction
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

**FIGURE 1** | Genome- and chromosome-wide association results. (A) Genome-wide association results for litter size in the NV sheep population. Manhattan plot shows the combined association signals [-log₁₀(*p*-value)] on the y-axis vs. SNPs position in the sheep genome on the x-axis and ordered by chromosome number (assembly OARv3.1). Red line represents the 5% genome-wide threshold. (B) OARX chromosome-wide association results. The curve shows the combined association signals [-log₁₀(*p*-value)] on the y-axis vs. SNPs position on the X chromosome on the x-axis (assembly OARv3.1). Red line represents the 5% chromosome-wide threshold. The *BMP15* gene location is indicated by a red arrowhead.

(Figure 1A and Table 1). Importantly, at the chromosome-wide level, a cluster of 26 significant markers encompassed the location of the *BMP15* candidate gene (Figure 1B). In order to

better characterize this locus on the X chromosome, we have determined, for each individual, the most likely linkage phase across 80 markers (10 Mb) including the significant region.
<table><caption>**TABLE 1** | Markers significantly associated with litter size.</caption><thead><tr><th>SNP</th><th>Chromosome</th><th>Position<sup>a</sup></th><th>MAF<sup>b</sup></th><th>P<sub>Unadj</sub><sup>c</sup></th><th>P<sub>Chrom</sub><sup>d</sup></th><th>P<sub>Genome</sub><sup>e</sup></th></tr></thead><tbody><tr><td>OARX_51294776.1</td><td>OARX</td><td>53,756,339</td><td>0.28</td><td>3.52E-11</td><td>4.18E-08</td><td>1.67E-06</td></tr><tr><td>s27837.1</td><td>OARX</td><td>53,825,247</td><td>0.30</td><td>5.83E-10</td><td>6.92E-07</td><td>2.77E-05</td></tr><tr><td>s73460.1</td><td>OARX</td><td>53,905,939</td><td>0.44</td><td>9.59E-10</td><td>1.14E-06</td><td>4.55E-05</td></tr><tr><td>s39212.1</td><td>OARX</td><td>53,852,735</td><td>0.44</td><td>3.46E-09</td><td>4.11E-06</td><td>1.64E-04</td></tr><tr><td>OARX_52608221.1</td><td>OARX</td><td>52,367,253</td><td>0.41</td><td>1.68E-08</td><td>2.00E-05</td><td>7.99E-04</td></tr><tr><td>s46003.1</td><td>OARX</td><td>80,222,479</td><td>0.44</td><td>2.56E-07</td><td>3.03E-04</td><td>1.21E-02</td></tr><tr><td>OARX_55032299.1</td><td>OARX</td><td>48,942,926</td><td>0.30</td><td>1.57E-06</td><td>1.86E-03</td><td>NS</td></tr><tr><td>OARX_72164491.1</td><td>OARX</td><td>74,388,397</td><td>0.23</td><td>1.71E-06</td><td>2.03E-03</td><td>NS</td></tr><tr><td>OARX_49135019.1</td><td>OARX</td><td>42,475,099</td><td>0.35</td><td>1.80E-06</td><td>2.13E-03</td><td>NS</td></tr><tr><td>OARX_111306030.1</td><td>OARX</td><td>920,41,520</td><td>0.27</td><td>2.25E-06</td><td>2.67E-03</td><td>NS</td></tr><tr><td>OARX_102620828.1</td><td>OARX</td><td>82,796,975</td><td>0.41</td><td>3.02E-06</td><td>3.59E-03</td><td>NS</td></tr><tr><td>s31917.1</td><td>OARX</td><td>58,202,482</td><td>0.44</td><td>3.12E-06</td><td>3.70E-03</td><td>NS</td></tr><tr><td>OARX_72351736.1</td><td>OARX</td><td>74,590,448</td><td>0.19</td><td>5.63E-06</td><td>6.68E-03</td><td>NS</td></tr><tr><td>OARX_72263548.1</td><td>OARX</td><td>74,498,463</td><td>0.17</td><td>6.45E-06</td><td>7.66E-03</td><td>NS</td></tr><tr><td>OARX_49564109.1</td><td>OARX</td><td>42,876,169</td><td>0.49</td><td>8.16E-06</td><td>9.69E-03</td><td>NS</td></tr><tr><td>DU400878_520.1</td><td>OARX</td><td>73,847,207</td><td>0.27</td><td>9.96E-06</td><td>1.18E-02</td><td>NS</td></tr><tr><td>s54281.1</td><td>OARX</td><td>58,993,959</td><td>0.24</td><td>1.20E-05</td><td>1.42E-02</td><td>NS</td></tr><tr><td>s05229.1</td><td>OARX</td><td>58,346,644</td><td>0.46</td><td>1.71E-05</td><td>2.03E-02</td><td>NS</td></tr><tr><td>s27938.1</td><td>OARX</td><td>53,275,559</td><td>0.35</td><td>1.77E-05</td><td>2.10E-02</td><td>NS</td></tr><tr><td>OARX_54104393.1</td><td>OARX</td><td>49,870,983</td><td>0.36</td><td>2.03E-05</td><td>2.41E-02</td><td>NS</td></tr><tr><td>OARX_72236232.1</td><td>OARX</td><td>74,464,263</td><td>0.18</td><td>2.23E-05</td><td>2.64E-02</td><td>NS</td></tr><tr><td>OARX_111349974.1</td><td>OARX</td><td>92,085,118</td><td>0.18</td><td>2.23E-05</td><td>2.65E-02</td><td>NS</td></tr><tr><td>OARX_53703822.1</td><td>OARX</td><td>51,193,144</td><td>0.35</td><td>2.59E-05</td><td>3.08E-02</td><td>NS</td></tr><tr><td>OARX_43227227.1</td><td>OARX</td><td>36,235,514</td><td>0.32</td><td>3.41E-05</td><td>4.04E-02</td><td>NS</td></tr><tr><td>OARX_51842287.1</td><td>OARX</td><td>53,162,079</td><td>0.49</td><td>3.43E-05</td><td>4.08E-02</td><td>NS</td></tr><tr><td>OARX_102654502.1</td><td>OARX</td><td>82,837,158</td><td>0.20</td><td>3.65E-05</td><td>4.33E-02</td><td>NS</td></tr></tbody></table>

After haplotype clusterization, a specific segment of 3.5 Mb (50,639,087–54,114,793 bp, OARv3.1 genome assembly) was identified to be more frequent in highly prolific cases than in controls ($f_{cases}$ = 0.51 vs. $f_{controls}$ = 0.037, $P$ = 1.92E-11), Chi-square test (Figure 2). This identified segment contained the BMP15 gene (50,970,938–50,977,454 bp, OARv3.1) well known to play a crucial role in ovarian function and to be a target of numerous mutations in its coding region controlling prolificacy (Persani et al., 2014).

## Characterization of the Mutation

Although the BMP15 gene could be considered as a positional and functional candidate gene, no mutation was evidenced by Sanger sequencing of the BMP15 coding regions of the most prolific ewes studied. In order to find the potential causal mutation, we sequenced the whole genome of two specifically chosen ewes based on the shortest haplotype within the region (homozygous reference vs. homozygous variant) and their opposite extreme phenotypes (LS 1.1 vs. 2.8).

Variant search analysis and annotations through the GATK toolkit was limited to the OARX: 50,639,087–54,114,793 bp region. We detected 60 SNPs and 90 small insertions and deletions (INDELs) with quality score > 30 (Supplementary

**Table S1).** Among them, we particularly focused on the 85 variants located within annotated genes (upstream, exon, intron, splice acceptor or donor, and downstream localization). After filtering these 85 variants for allele sharing with other breeds based on SheepGenome DB<fnref n="5" /> and 68 publicly available domestic sheep genomes (International Sheep Genomics Consortium)<fnref n="6" />, none of them were removed, all being NV breed-specific. Finally, and based on prolificacy gene knowledge, we were particularly interested in one SNP (T > A) identified in the upstream region of the BMP15 gene at position 50,977,717 bp on OARX v3.1. We then developed an RFLP assay to specifically genotype for this polymorphism. Among the 79 animals of the GWAS, 31 ewes were heterozygous and six homozygous for the A variant allele. As shown in **Table 2**, most of the A carrier ewes were in the highly prolific Case group (34 among 39), and only three set in the Control group. When associating the LS performance of the 79 ewes to their genotype at the OARX: 50977717T > A SNP, the A non-carriers exhibited a mean LS of 1.36, heterozygous T/A a mean LS of 2.32, and homozygous A/A a mean LS of 2.73, indicating that the A allele of this polymorphism was strongly associated with increased LS in NV (T/A or A/A vs. T/T,

<footnote><sup>a</sup>Position of markers are based on the OARv3.1 assembly in bp. <sup>b</sup>MAF, minor allele frequency. <sup>c</sup>P<sub>Unadj</sub> corresponds to exact unadjusted p-value for the Fisher's test.</footnote>

<footnote><sup>d</sup>P<sub>Chrom</sub> corresponds to p-value after chromosome-wide Bonferroni correction. <sup>e</sup>P<sub>Genome</sub> corresponds to p-value after genome-wide Bonferroni correction (NS, non-significant).</footnote>

<footnote n="5">⁵http://sheepgenomesdb.org/</footnote>

<footnote n="6">⁶http://www.sheephapmap.org/</footnote>
<figure><img src="image_5.png" /></figure>

$P < 1E^{-3}$, one-way ANOVA). Furthermore, this polymorphism appears to be in total linkage disequilibrium with the six more significant markers from the GWAS analysis (Figure 3). Genotype information at the OARX: 50977717T > A locus was introduced in the GWAS analysis. This SNP appeared as the most significant marker associated with the prolificacy phenotype ($P_{unadjusted} = 1.93E^{-11}$, $P_{Chromosome-wide corrected} = 1.62E^{-14}$ and $P_{Genome-wide corrected} = 9.13E^{-07}$), suggesting that it could be the

causal mutation (Supplementary Figure S1). In accordance with the *Fec* gene nomenclature, the mutant allele identified upstream of the *BMP15* gene in NV sheep was named *FecX*<sup>N</sup>.

As described for other prolific alleles, such as *FecB*<sup>B</sup>, *FecX*<sup>G</sup>, *FecG*<sup>H</sup>, *FecX*<sup>Gr</sup>, and *FecL*<sup>L</sup>, a given mutation can segregate in several sheep populations (Davis et al., 2002; Mullen et al., 2013; Chantepie et al., 2018; Ben Jemaa et al., 2019). We have tested the *FecX*<sup>N</sup> allele presence in a diversity of 26 sheep
<table><caption>**TABLE 2** | Distribution of OARX:50977717T > A SNP genotypes and associated LS in case and control groups.</caption><thead><tr><th>Group</th><th>Genotype</th><th>TT</th><th>TA</th><th>AA</th></tr></thead><tbody><tr><td rowspan="2">Low prolific control</td><td>n</td><td>37</td><td>3</td><td></td></tr><tr><td>Raw mean LS</td><td>1.22</td><td>1.39</td><td></td></tr><tr><td rowspan="2">Highly prolific case</td><td>n</td><td>5</td><td>28</td><td>6</td></tr><tr><td>Raw mean LS</td><td>2.41</td><td>2.43</td><td>2.73</td></tr><tr><td>Total</td><td></td><td>42</td><td>31</td><td>6</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**FIGURE 3** | Linkage disequilibrium (LD) plot. Generated by Haploview, pairwise LD between SNP markers (OARX: 50977717T > A and the five genome-wide significant LS-associated markers) is represented by the D' value. Strong LD is represented in dark red boxes (D' = 100) and weaker LD (D' < 100) in lighter red boxes.</figcaption></figure>

breeds representing 725 animals (Rochus et al., 2018). Among the breeds tested, the *FecX*<sup>N</sup> genotyping has confirmed the segregation of this mutation in NV breed and revealed its presence in the Blanche du Massif Central (BMC) and Lacaune breeds (**Table 3**). Additionally, the *FecX*<sup>N</sup> variant was absent from the Ensembl variant database<sup>7</sup> compiling information from (i) dbSNP, (ii) whole-genome sequencing information from the NextGen project (180 animals from various Iranian and Moroccan breeds), and (iii) the International Sheep Genome Consortium (551 animals from 39 breeds all over the world).

## *FecX*<sup>N</sup> Genotype Frequency and Effect on Prolificacy

Large cohorts of ewes, chosen at random, were genotyped in order to accurately estimate the allele frequencies in the NV and the BMC populations (**Table 4**). The frequency of the N prolific

allele at the *FecX* locus was similar in both populations, 0.11 and 0.10, with a distribution of 19.4 and 17.6% heterozygous, 1.5 and 1% homozygous carriers in NV and BMC, respectively. The genotype frequencies were consistent with the Hardy-Weinberg equilibrium (HWE) in both breeds (NV *P* = 0.28 and BMC *P* = 0.76).

Based on the raw mean LS observations, the *FecX*<sup>N</sup> carrier ewes clearly exhibited increased LS compared to non-carriers in both populations (**Table 4**). The L prolific allele at the *FecL* locus is also segregating in NV (Chantepie et al., 2018). Results of the linear mixed model showed that, for the NV breed, one copy of the *FecX*<sup>N</sup> allele significantly increased LS by + 0.22, and two copies increased LS by +0.65, and a single copy of the *FecL*<sup>L</sup> allele increased LS by + 0.41, and two copies by + 0.72. Based on the 80 ewes genotyped heterozygous at both loci, it appeared that the effects of *FecX*<sup>N</sup> and *FecL*<sup>L</sup> on LS were not fully additive, the expected LS being significantly slightly reduced by -0.05 (0.58 instead of 0.63) due to a significant interaction between the two loci (*P* < 0.001) (**Figure 4A**). For the BMC population, compared to *FecX*<sup>+</sup>/*FecX*<sup>+</sup> ewes, *FecX*<sup>N</sup>/*FecX*<sup>+</sup> exhibited increased LS by + 0.18 and *FecX*<sup>N</sup>/*FecX*<sup>N</sup> by + 0.30 under natural estrus (**Figure 4B**). The use of PMSG for estrus synchronization increased LS significantly among *FecX*<sup>+</sup>/*FecX*<sup>+</sup> ewes (+ 0.23) and *FecX*<sup>N</sup>/*FecX*<sup>+</sup> ewes (+ 0.18) while the effect on *FecX*<sup>N</sup>/*FecX*<sup>N</sup> ewes was negative (-0.13). The combined effect of the first copy of the *FecX*<sup>N</sup> allele and the use of PMSG treatment was not fully additive, the interaction being significant (*P* < 0.001) although low (0.35 instead of 0.41) (**Figure 4B**).

## Functional Effects of the *FecX*<sup>N</sup> Mutation

As described above, *FecX*<sup>N</sup> is located upstream of the coding region of the *BMP15* gene when referencing to the ovine genome v3.1 (ensembL.org) or v4.0 (ncbi.nlm.nih.gov). In both versions of the ovine genome, the *BMP15* gene annotation begins at the ATG start site and *FecX*<sup>N</sup> is located 290 pb upstream, possibly in the 5'UTR and/or the proximal promoter region. As a first approach, we took advantage of RNA sequencing data from ovine oocytes publicly available at EMBL-EBI (Bonnet et al., 2013). After reads mapping against the ovine genome (v3.1) using STAR2 aligner within the Galaxy pipeline and visualization with Integrative Genome Viewer (IGV), **Figure 5** shows the location of *FecX*<sup>N</sup> within the possible 5'UTR of the *BMP15* gene when expressed in the oocyte. Consequently, we have first tested the potential functional impact of *FecX*<sup>N</sup> on the *in vitro* stability and translatability of the *BMP15* mRNA. Thus, the reference (T, *FecX*<sup>+</sup>) and variant (A, *FecX*<sup>N</sup>) forms of the ovine BMP15 cDNA (-297, +1183 referring to ATG start codon) were cloned in a pGEM-T vector for subsequent *in vitro* T7 promoter-dependent transcription/translation experiment using reticulocyte lysate solution. As shown in **Figure 6**, the Western blotting of the BMP15 proteins produced from both forms, and their chemiluminescent quantification revealed that the *FecX*<sup>N</sup> mutation had no significant impact on the overall stability and translatability of the *BMP15* mRNA in this condition.

As a second hypothesis, we have tested the *FecX*<sup>N</sup> impact on the *BMP15* promoter activity. Two promoter regions were

<footnote>7http://www.ensembL.org</footnote>
<table><caption>**TABLE 3** | *FecX*<sup>N</sup> genotype distribution from a diversity panel of French ovine breeds.</caption><thead><tr><th rowspan="2">Breed</th><th rowspan="2">Total</th><th colspan="2">Genotype<sup>a</sup></th><th rowspan="2">Breed</th><th rowspan="2">Total</th><th colspan="2">Genotype<sup>a</sup></th></tr><tr><th>+/+ (+/Y)</th><th>N/+ (N/Y)</th><th>+/+ (+/Y)</th><th>N/+ (N/Y)</th></tr></thead><tbody><tr><td>Berrichon du Cher</td><td>29</td><td>29</td><td></td><td>Mourerous</td><td>26</td><td>26</td><td></td></tr><tr><td>Blanche Massif Central</td><td>31</td><td>27</td><td>4</td><td>Mouton Vendéen</td><td>30</td><td>30</td><td></td></tr><tr><td>Causse du Lot</td><td>32</td><td>32</td><td></td><td>Noire du Velay</td><td>28</td><td>26</td><td>2</td></tr><tr><td>Charmoise</td><td>31</td><td>31</td><td></td><td>Préalpes du sud</td><td>27</td><td>27</td><td></td></tr><tr><td>Charollais</td><td>29</td><td>29</td><td></td><td>Rava</td><td>29</td><td>29</td><td></td></tr><tr><td>Corse</td><td>30</td><td>30</td><td></td><td>Romane</td><td>29</td><td>29</td><td></td></tr><tr><td>Ile de France</td><td>28</td><td>28</td><td></td><td>Romanov</td><td>26</td><td>26</td><td></td></tr><tr><td>Lacaune (meat)</td><td>42</td><td>40</td><td>2</td><td>Rouge de l'Ouest</td><td>28</td><td>28</td><td></td></tr><tr><td>Lacaune (dairy)</td><td>40</td><td>40</td><td></td><td>Roussin</td><td>30</td><td>30</td><td></td></tr><tr><td>Limousine</td><td>30</td><td>30</td><td></td><td>Suffolk</td><td>20</td><td>20</td><td></td></tr><tr><td>Manech tête rousse</td><td>29</td><td>29</td><td></td><td>Tarasconnaise</td><td>32</td><td>32</td><td></td></tr><tr><td>Martinik</td><td>22</td><td>22</td><td></td><td>Texel</td><td>21</td><td>21</td><td></td></tr><tr><td>Merinos d'Arles</td><td>26</td><td>26</td><td></td><td></td><td></td><td></td><td></td></tr></tbody><tfoot><tr><td></td><td></td><td></td><td></td><td>Total</td><td>725</td><td>717</td><td>8</td></tr></tfoot></table>

<table><caption>**TABLE 4** | *FecX*<sup>N</sup> frequencies and effects on LS in NV and BMC breeds.</caption><thead><tr><th rowspan="3">FecX genotype</th><th colspan="6">Breed</th></tr><tr><th colspan="3">NV (n = 2323)</th><th colspan="3">BMC (n = 2456)</th></tr><tr><th>+/+</th><th>N/+</th><th>N/N</th><th>+/+</th><th>N/+</th><th>N/N</th></tr></thead><tbody><tr><td>Number of animals</td><td>1839</td><td>450</td><td>34</td><td>1999</td><td>432</td><td>25</td></tr><tr><td>Frequency (%)</td><td>79.2</td><td>19.4</td><td>1.5</td><td>81.3</td><td>17.6</td><td>1</td></tr><tr><td>Raw mean LS</td><td>1.66</td><td>1.93</td><td>2.45</td><td>1.56</td><td>1.79</td><td>1.87</td></tr><tr><td colspan="7"><strong>ANOVA solutions<sup>a</sup></strong></td></tr><tr><td>FecL genotype +/+</td><td>0.00</td><td>0.22 (0.02)</td><td>0.65 (0.08)</td><td>NE</td><td>0.00</td><td>0.18 (0.02)</td><td>0.30 (0.09)</td></tr><tr><td>L/+</td><td>0.41 (0.02)</td><td>0.58 (0.04)</td><td>0.56 (0.16)</td><td>IE</td><td>0.23 (0.01)</td><td>0.35 (0.03)</td><td>0.17 (0.11)</td></tr><tr><td>L/L</td><td>0.72 (0.08)</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

tested [(-743, -11) bp and (-443, -102) bp referring to ATG start codon] cloned in front of the luciferase reporter gene and transiently expressed in CHO cells cultured *in vitro*. As shown by the luciferase assays (Figure 7), the *FecX*<sup>N</sup> variant was able to significantly reduce the luciferase activity in the context of both the long or the short BMP15 promoters, indicating the possible inhibitory impact of *FecX*<sup>N</sup> on BMP15 gene expression. To go further with this hypothesis, *in vivo* BMP15 gene expression was measured directly on isolated oocytes pools from NV and BMC homozygous ewe carriers and non-carriers of the *FecX*<sup>N</sup> allele. Even if the mean expression of BMP15 checked by real-time qPCR was decreased by twofold in the oocytes of *FecX*<sup>N</sup> carriers, this difference appeared as not significant (P = 0.166, genotype effect, two-way ANOVA) due to a large variability between animals. In parallel, the mean expression of the second oocyte-specific prolificacy major gene GDF9 seemed unaffected by the presence of *FecX*<sup>N</sup> (Figure 8).

DISCUSSION

The present study identified the g.50977717T > A variant on the ovine chromosome X upstream of the BMP15 gene as the most likely causative mutation for the increased prolificacy of the NV ewes. The highly significant genetic association with the extreme LS phenotype, the significant effect of the A variant on increasing prolificacy by + 0.2 lamb per lambing in a large set of NV ewes, also found in the BMC genetic background, and the demonstrated action on BMP15 transcriptional activity all support the causality of this mutation named *FecX*<sup>N</sup>.

The BMP15 gene is at the top of the list of candidate genes controlling ovarian function, ovulation rate, and thus prolificacy in the ovine species with 10 independent mutations identified out of the 17 already known. Indeed, eight SNPs and two small INDELS, all within the open reading frame, were evidenced affecting prolificacy and BMP15 function. Among these mutations, two SNPs and the two INDELS impaired the

<footnote><sup>a</sup>+/+ or N/+ females, N/Y or N/+ hemizygous males.</footnote>

<footnote><sup>a</sup>Estimated increased LS (± standard error) compared to L and N non-carriers from the linear mixed models with interaction between genotypes at FecX and FecL loci for Noire du Velay (NV) and between FecX genotypes and estrus type (NE, natural estrus; IE, Induced estrus) for Blanche du Massif Central (BMC). Global effects for genotypes, estrus type and their interactions were all significant (P < 0.001).</footnote>
<figure><img src="image_7.png" /></figure>

**FIGURE 4** | Estimated effect of *FecX*<sup>N</sup>, (A) Estimated effect of *FecX*<sup>N</sup> knowing *FecL*<sup>L</sup> genotype in Noire du Velay (NV). (B) Estimated effect of *FecX*<sup>N</sup> according to induced or natural estrus in Blanche du Massif Central (BMC). Estimated increased LS and the corresponding standard errors from the linear mixed models with interaction between genotypes at the *FecX* and *FecL* loci for NV, and between *FecX* genotype and estrus type for BMC are given relative to *L* and *N* alleles non-carrier ewes (+/+). Different letters indicate significant differences (*P* < 0.01).

<figure><img src="image_8.png" /></figure>

**FIGURE 5** | Localization of *FecX*<sup>N</sup> in the *BMP15* upstream region. Integrative Genomics Viewer (IGV v2.4.10) snapshot of the ovine *BMP15* gene upstream region alignments (coordinates from OAR v3.1 assembly) with reads from whole genome DNA sequencing (DNASeq) of a homozygous *FecX*<sup>N</sup> carrier ewe (green line) and total mRNA sequencing of ovine oocytes from small antral follicles (RNASeq, Bonnet et al., 2013), indicating a possible localization of *FecX*<sup>N</sup> in the 5'UTR region of *BMP15* (ATG start codon is indicated by an arrowhead and *BMP15* gene annotation is on the minus strand).

protein production either by generating premature stop codon (*FecX*<sup>H</sup>, Galloway et al., 2000; *FecX*<sup>G</sup>, Hanrahan et al., 2004) or by breaking the reading frame (*FecX*<sup>R</sup>, Martinez-Royo et al., 2008, *FecX*<sup>Bar</sup>, Lassoued et al., 2017). The six other SNPs generate non-conservative amino acid substitutions, one being not functionally tested (L252P, Amini et al., 2018), and all others leading to a loss of function of BMP15 ranging from inhibited protein production (*FecX*<sup>L</sup>, Bodin et al., 2007) and impaired interaction with GDF9 (*FecX*<sup>I</sup> and *FecX*<sup>B</sup>, Liao et al., 2004) to altered

cell signaling activity (*FecX*<sup>Gr</sup> and *FecX*<sup>O</sup>, Demars et al., 2013). In contrast with the 10 mutations described above, the *FecX*<sup>N</sup> variant evidenced in the present study is not located in the open reading frame of *BMP15* and does not alter the protein sequence. However, no other polymorphism genetically linked to *FecX*<sup>N</sup> was found in the *BMP15* coding sequence when checked by whole genome or local Sanger sequencing of the *BMP15* gene from *FecX*<sup>N</sup> carrier animals. Of course, this does not rule out the possibility of a polymorphism lying in another gene nearby
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

**FIGURE 6** | Effect of *FecX*^N^ mutation on the ovine BMP15 protein produced *in vitro*. *In vitro* transcription and translation were realized in duplicate in six independent experiments from antisense *BMP15* cDNA carrying the wild-type allele (antisense +) as negative control or sense *BMP15* cDNA with (sense N) or without *FecX*^N^ (sense +). (A) A representative western blot experiment revealing the BMP15 protein. (B) Chemiluminescent signal of BMP15 was captured by a ChemiDoc MP imaging system and images were quantified (relative to sense +) and analyzed with the Image Lab Software (Bio-Rad). Means (± SEM) were not significantly different (NS) after Student's T-test with Welch's correction.

<figure><img src="image_11.png" /></figure>

**FIGURE 7** | Functional effect of *FecX*^N^ mutation on the *BMP15* promoter activity. *In vitro* reporter luciferase assay from CHO cells transiently transfected with empty vector or wild-type *BMP15* promoter (*BMP15*+) or mutant *BMP15* promoter (*BMP15*^N^). Two fragments were generated, a long (lg) form of 732 bp (-743, -11 bp referring to ATG start codon) and a short (sh) form of 341 bp (-443, -102 bp). Results are expressed as means ± SEM of the relative light unit (RLU) from six independent transfection experiments in triplicate. Asterisk indicates significant difference after one-way ANOVA \*P < 0.05; \*\*P < 0.01; \*\*\*P < 0.0001. MP, minimal promoter of pGL4.23; Luc 2, luciferase reporter gene; hatched bar, *BMP15* promoter; red line, *FecX*^N^ mutation.

with a still unknown role in the ovarian function and prolificacy. Nevertheless, we did not find any polymorphism (SNP and INDEL) altering the coding sequence of genes annotated in the significantly LS-associated genetic region of 3.5 Mb on OARX (Supplementary Table S1), leaving *BMP15* as the most obvious candidate.

Whatever the version of the ovine reference genome (Oar_v3.1, Oar_v4.0, or even the last Oar_rambouillet_v1.0), the annotation of the *BMP15* gene always starts at the ATG-initiating codon. Using publicly available transcriptome data from ovine oocytes RNAseq analysis, we were able to show that *FecX*^N^ located 290 bp upstream of *BMP15* could stand
<figure><img src="image_12.png" /><figcaption>**FIGURE 8** | Effect of *FecX*<sup>N</sup> mutation on *BMP15* and *GDF9* expression in ovine oocyte. Quantitative real-time PCR results of *BMP15* (**A**) and *GDF9* (**B**) expression in oocyte pools from growing (1–3 mm) follicles (two pools of five oocytes per animal) from homozygous *FecX*<sup>N</sup> carrier (*N/N*, *n* = 5) and non-carrier ewes (+/+, *n* = 5) during the follicular phase of the estrus cycle. Results are expressed as means ± SEM of the mRNA relative level for each genotype, using *GAPDH* and *SDHA* as internal references. Raw data were analyzed by two-way ANOVA. Means (± SEM) were not significantly different (NS).</figcaption></figure>

in its 5'UTR region. From our *in vitro* functional analyses, *FecX*<sup>N</sup> was not demonstrated to influence the translatability of the *BMP15* mRNA, but on the contrary, it was shown to decrease the *BMP15* promoter activity. Little is known about transcription factors able to regulate *BMP15* expression. Several regulatory elements were evidenced in pig *BMP15* promoter hosting consensus-binding sites for LHX8, NOBOX, and PITX1 transcription factors. However, only LHX8 was demonstrated as functionally activating the porcine BMP15 promoter activity (Wan et al., 2015). In human, a regulatory mutation in the 5'UTR of *BMP15* (c.-9C > G) was associated to non-syndromic premature ovarian failure (Dixit et al., 2006), and also iatrogenic ovarian hyperstimulation syndrome (Morón et al., 2006). This mutation was shown to enhance the fixation of the PITX1 factor transactivating the BMP15 promoter (Fonseca et al., 2014). However, the *FecX*<sup>N</sup> position does not fit with the syntenic location of porcine LHX8 and human PITX1 binding sites on the ovine *BMP15* promoter. Using the MatInspector promoter analysis tool (Genomatix), we were only able to hypothesize an alteration by *FecX*<sup>N</sup> of a putative TATA-box-like sequence (TTAATA > TTATA). Unfortunately, our electromobility shift assay attempts using CHO nuclear extracts failed to demonstrate the binding of any factor at the *FecX*<sup>N</sup> position, preventing us from defining the precise molecular mechanism by which *FecX*<sup>N</sup> decreases *BMP15* promoter activity.

The inhibition of the promoter activity suggests a transcriptional regulatory role of *FecX*<sup>N</sup> on *BMP15* gene expression. That is the reason why *BMP15* mRNA accumulation in homozygous *FecX*<sup>N</sup>/*FecX*<sup>N</sup> oocytes was checked and compared to non-carrier oocytes. The twofold decrease of

*BMP15* mean expression in oocytes of *FecX*<sup>N</sup> carriers was not statistically significant. However, the moment we have chosen during the follicular phase of the late folliculogenesis for the comparative analysis between *FecX*<sup>+</sup> and *FecX*<sup>N</sup> oocytes from 1 to 3 mm antral follicles without controlling the atretic status might not be optimal to visualize a highly significant differential expression of *BMP15*. To this could be added an induced effect on oocyte by our experimental procedures, progestagen treatment for estrus cycle synchronization and/or oocyte denudation (Rose et al., 2013; Menchaca et al., 2018). The *BMP15* gene expression in ovine oocytes begins during the primary stage of follicular development, and its expression increases up to the antral stages where it stays stable (McNatty et al., 2005; Bonnet et al., 2011; Kona et al., 2016; Bonnet and Fabre, unpublished data). Moreover, the streak ovaries phenotype of infertile ewes carrying homozygous mutations in *BMP15* have evidenced its crucial role in controlling the primary to secondary follicle transition (Galloway et al., 2000; Bodin et al., 2007; Lassoued et al., 2017). Consequently, it would certainly be appropriate to follow the *BMP15* expression in *FecX*<sup>N</sup> carrier ewes from these early stages of folliculogenesis to better decipher the mutation impact on ovarian physiology. Nevertheless, the fact that *FecX*<sup>N</sup> could inhibit the *BMP15* gene expression fits well with the physiological and molecular models associating BMP system loss of function and increased sheep prolificacy (Fabre et al., 2006; Demars et al., 2013).

One copy of the *FecX*<sup>N</sup> allele significantly increased by +0.30 to +0.50 the raw mean LS of NV ewes. When corrected for different environmental effects and more particularly for the genotype at the *FecL* locus, the estimated effect of *FecX*<sup>N</sup> on LS was +0.22 lamb per lambing for the first copy and + 0.43 for the second copy. The *FecX*<sup>N</sup> effect is in the range of already known prolific alleles in various sheep breeds (Jansson, 2014) and explains about 20% of the genetic variance of prolificacy without considering information at the *FecL* locus and up to 50% when polygenic effect is corrected for the presence of *FecL*<sup>L</sup>, confirming its status of major gene. The effect of *FecX*<sup>N</sup> on LS seems relatively independent of the genetic background. Indeed, the estimated positive effect of the first copy of *FecX*<sup>N</sup> on prolificacy in NV was confirmed in the BMC breed with + 0.18 lamb per lambing based on natural estrus. Moreover, the same robust effect was observed even in the presence of PMSG for synchronizing the estrus cycles preceding the lambing. The same observation is made for other mutations controlling sheep prolificacy. For instance, the *FecL*<sup>L</sup> allele exhibited a similar effect on LS in NV (+0.41, present study; +0.42, Chantepie et al., 2018) Lacaune, +0.47 (Martin et al., 2014), and D'man (+0.30, Ben Jemaa et al., 2019), and this was also observed for the *FecB*<sup>B</sup> allele introgressed in several populations (Kumar et al., 2008).

By genotyping a diversity panel, we also evidenced the presence of 2 *FecX*<sup>N</sup> carrier animals in the Lacaune meat strain, which will require further genotyping of numerous animals. If this is confirmed, the Lacaune meat breed will be another population, as Belclare, in which three different natural prolific mutations are segregating (Hanrahan et al., 2004; Bodin et al.,
2007; Drouilhet et al., 2013). The presence of *FecgL*^L^ in both NV and Lacaune and the presence of *Fecx*^N^ in NV, BMC, and Lacaune also raise the question of the origin of these mutations. From population structure analysis, it was shown that NV, BMC, and Lacaune shared the same origin within the European southern sheep populations that may explain the segregation of the same mutations in these populations (Rochus et al., 2018).

## CONCLUSION

In conclusion, through a case/control GWAS strategy and genome sequencing, we have identified in the NV breed a second prolific mutation named *Fecx*^N^ affecting the expression of the *BMP15* gene, a well-known candidate gene controlling OR and LS in sheep. This work confirms the relevance of the whole genome approaches to decipher the genetic determinism of the prolificacy trait. Homozygous *Fecx*^N^/*Fecx*^N^ animals were still hyperprolific as previously observed for *Fecx*^Gr^ and *Fecx*^O^ but in contrast with sterile animals observed for the seven other *Fecx* homozygous variants in *BMP15*. As an upstream regulatory mutation, *Fecx*^N^ also contrasts with these nine other prolific causal mutations all evidenced in the coding part of *BMP15* and altering the protein function. Thanks to this new sheep model, the genetic etiology of ovarian pathologies in women could be improved by searching polymorphisms, not only in the coding region, but also in the regulatory parts driving the *BMP15* expression within the oocyte.

## DATA AVAILABILITY STATEMENT

The datasets generated for this study can be found in the European Nucleotide Archive (ENA) at EMBL-EBI (accession number PRJEB35553).

## ETHICS STATEMENT

All experimental procedures were approved (approval number 01171.02) by the French Ministry of Teaching and Scientific Research and local ethical committee C2EA-115 (Science and Animal Health) in accordance with the European Union Directive 2010/63/EU on the protection of animals used for scientific purposes. Written informed consent for participation was not obtained from the owners because breeders in selection schemes are affiliated to selection organisms and give their tacit agreement for animal participation to these organisms collaborating with research institutes.

## AUTHOR CONTRIBUTIONS

LC conducted the investigation, validation, formal analysis and visualization of the data, and wrote the original draft. LB participated in the conceptualization, formal analysis,

data curation, and final writing, review, and editing of the manuscript. JS, FW, FP-P, JR, and LD participated in the investigation and getting the biological resources. LD helped in writing, review, and editing of the final manuscript. SF was implicated in the conceptualization, supervision, funding acquisition, investigation, formal analysis and writing, review, and editing of the final manuscript. All authors contributed to the article and approved the submitted version.

## FUNDING

LC was supported by a Ph.D. grant co-funded by the APIS-GENE through the Proligen project and the European Funds for Regional Development (FEDER) through the Interreg POCTEFA programme in the framework of the PIRINNOVI project (EFA103/15). Part of the NV sampling was supported by the DEGERAM project co-funded by the FEDER Massif Central, the Règions: Aquitaine, Midi-Pyrènèes, Limousin and Auvergne, and the French government.

## ACKNOWLEDGMENTS

This manuscript has been released as a preprint at BioRxiv (Chantepie et al., 2019). We thank Claire Chantaduc, Didier Cathalan, and Kévin Chile from ROM Seilection managing the NV and BMC populations for their precious help in the planning of blood sampling. We grateful to the breeders who made their animals available for this study.

## SUPPLEMENTARY MATERIAL

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

**FIGURE S1** | Genome-wide and chromosome-wide association results integrating the SNP OARX: 50977717T > A. (A) The SNP OARX: 50977717T > A genotypes were included in the SNP50 Beadchip data for genome-wide association analysis for litter size in the NV sheep population. Manhattan plot shows the combined association signals [-log₁₀(*p*-value)] on the y-axis vs. SNPs position in the sheep genome on the x-axis and ordered by chromosome number (assembly OARv3.1). Red line represents the 5% genome-wide threshold. (B) OARX chromosome-wide association results. The curve shows the combined association signals [-log₁₀(*p*-value)] on the y-axis vs. SNPs position on the X chromosome on the x-axis (assembly OARv3.1). Red line represents the 5% chromosome-wide threshold. In both panels, the position of the SNP OARX:50977717T > A is indicated by a red dot. In (B), the *BMP15* gene location is indicated by a red arrowhead.

**TABLE S1** | List of variants found in the OARX: 50639087-54114793 region. Listing of 60 SNPs and 90 small INDELs with quality score > 30.

**TABLE S2** | List of primers used in the study. Locations of primers are based on the OARv3.1 ovine genome assembly available on ensembl.org.
REFERENCES

Abdoli, R., Zamani, P., Mirhoseini, S. Z., Ghavi Hossein-Zadeh, N., and Nadri, S. (2016). A review on prolificacy genes in sheep. *Reprod. Domest. Anim.* 51, 631-637. doi: 10.1111/rda.12733

Amini, H.-R., Ajaki, A., Farahi, M., Heidari, M., Pirali, A., Forouzanfar, M., et al. (2018). The novel T755C mutation in BMP15 is associated with the litter size of Iranian Afshari. Ghezel, and Shal breeds. *Arch. Anim. Breed.* 61, 153–160. doi: 10.5194/aab-61-153-2018

Balding, D. J. (2006). A tutorial on statistical methods for population association studies. *Nat. Rev. Genet.* 7, 781–791. doi: 10.1038/nrg1916

Ben Jemaa, S., Ruesche, J., Sarry, J., Woloszyn, F., Lassoued, N., and Fabre, S. (2019). The high prolificacy of D'man sheep is associated with the segregation of the FecLL mutation in the B4GALNT2 gene. *Reprod. Domest. Anim.* 54, 531–537. doi: 10.1111/rda.13391

Bodin, L., Di Pasquale, E., Fabre, S., Bontoux, M., Monget, P., Persani, L., et al. (2007). A novel mutation in the bone morphogenetic protein 15 gene causing defective protein secretion is associated with both increased ovulation rate and sterility in lacaune sheep. *Endocrinology* 148, 393–400. doi: 10.1210/en.2006-0764

Bonnet, A., Bevilacqua, C., Benne, F., Bodin, L., Cotinot, C., Liaubet, L., et al. (2011). Transcriptome profiling of sheep granulosa cells and oocytes during early follicular development obtained by Laser Capture Microdissection. *BMC Genomics* 12:417. doi: 10.1186/1471-2164-12-417

Bonnet, A., Cabau, C., Bouchez, O., Sarry, J., Marsaud, N., Foissac, S., et al. (2013). An overview of gene expression dynamics during early ovarian folliculogenesis: specificity of follicular compartments and bi-directional dialog. *BMC Genomics* 14:904. doi: 10.1186/1471-2164-14-904

Chang, C. C., Chow, C. C., Tellier, L. C., Vattikuti, S., Purcell, S. M., and Lee, J. J. (2015). Second-generation PLINK: rising to the challenge of larger and richer datasets. *Gigascience* 4:7. doi: 10.1186/s13742-015-0047-8

Chantepie, L., Bodin, L., Sarry, J., Woloszyn, F., Plisson-Petit, F., Ruesche, J., et al. (2018). Presence of causative mutations affecting prolificacy in the noire du velay and mouton vendeen sheep breeds. *Livestock Sci.* 216, 44–50. doi: 10.1016/j.livsci.2018.07.007

Chantepie, L., Bodin, L., Sarry, J., Woloszyn, F., Plisson-Petit, F., Ruesche, J., et al. (2019). Genome-wide identification of a regulatory mutation in BMP15 controlling prolificacy in sheep. *BioRxiv* [Preprint]. doi: 10.1101/858860

Davis, G. H., Galloway, S. M., Ross, I. K., Gregan, S. M., Ward, J., Nimbkar, B. V., et al. (2002). DNA tests in prolific sheep from eight countries provide new evidence on origin of the Booroola (FecB) mutation. *Biol. Reprod.* 66, 1869–1874. doi: 10.1095/biolreprod66.6.1869

Demars, J., Cano, M., Drouilhet, L., Plisson-Petit, F., Bardou, P., Fabre, S., et al. (2017). Genome-wide identification of the mutation underlying fleece variation and discriminating ancestral hairy species from modern woolly sheep. *Mol. Biol. Evol.* 34, 1722–1729. doi: 10.1093/molbev/msx114

Demars, J., Fabre, S., Sarry, J., Rossetti, R., Gilbert, H., Persani, L., et al. (2013). Genome-wide association studies identify two novel BMP15 mutations responsible for an atypical hyperprolificacy phenotype in sheep. *PLoS Genet.* 9:e1003482. doi: 10.1371/journal.pgen.1003482

DePristo, M. A., Banks, E., Poplin, R., Garimella, K. V., Maguire, J. R., Hartl, C., et al. (2011). A framework for variation discovery and genotyping using next-generation DNA sequencing data. *Nat. Genet.* 43, 491–498. doi: 10.1038/ng.806

Dixit, H., Rao, L. K., Padmalatha, V. V., Kanakavalli, M., Deenadayal, M., Gupta, N., et al. (2006). Missense mutations in the BMP15 gene are associated with ovarian failure. *Hum. Genet.* 119, 408–415. doi: 10.1007/s00439-006-0150-0

Drouilhet, L., Mansanet, C., Sarry, J., Tabet, K., Bardou, P., Woloszyn, F., et al. (2013). The highly prolific phenotype of Lacaune sheep is associated with an ectopic expression of the B4GALNT2 gene within the ovary. *PLoS Genet.* 9:e1003809. doi: 10.1371/journal.pgen.1003809

Elvin, J. A., Yan, C., Wang, P., Nishimori, K., and Matzuk, M. M. (1999). Molecular characterization of the follicle defects in the growth differentiation factor 9-deficient ovary. *Mol. Endocrinol.* 13, 1018–1034. doi: 10.1210/mend.13.6.0309

Fabre, S., Pierre, A., Mulsant, P., Bodin, L., Di Pasquale, E., Persani, L., et al. (2006). Regulation of ovulation rate in mammals: contribution of sheep genetic models. *Reprod. Biol. Endocrinol.* 4:20. doi: 10.1186/1477-7827-4-20

Fonseca, D. J., Ortega-Recalde, O., Esteban-Perez, C., Moreno-Ortiz, H., Patiño, L. C., Bermúdez, O. M., et al. (2014). BMP15 c.-9C>G promoter sequence variant may contribute to the cause of non-syndromic premature ovarian failure. *Reprod. BioMed. Online* 29, 627–633. doi: 10.1016/j.rbmo.2014.07.018

Galloway, S. M., McNatty, K. P., Cambridge, L. M., Laitinen, M. P., Juengel, J. L., Jokiranta, T. S., et al. (2000). Mutations in an oocyte-derived growth factor gene (BMP15) cause increased ovulation rate and infertility in a dosage-sensitive manner. *Nat. Genet.* 25, 279–283. doi: 10.1038/77033

Gilmour, A. R., Gogel, B. J., Cullis, B. R., and Thompson, R. (2009). *ASReml User-Guide Release 3.0*. Hemel Hempstead: VSN International Ltd.

Hanevik, H. I., Hilmarsen, H. T., Skjelbred, C. F., Tanbo, T., and Kahn, J. A. (2011). A single nucleotide polymorphism in BMP15 is associated with high response to ovarian stimulation. *Reprod. BioMed. Online* 23, 97–104. doi: 10.1016/j.rbmo.2011.02.015

Hanrahan, J. P., Gregan, S. M., Mulsant, P., Mullen, M., Davis, G. H., Powell, R., et al. (2004). Mutations in the genes for oocyte-derived growth factors GDF9 and BMP15 are associated with both increased ovulation rate and sterility in cambridge and belclare sheep (*Ovis aries*). *Biol. Reprod.* 70, 900–909. doi: 10.1095/biolreprod103.203093

Jansson, T. (2014). *Genes Involved in Ovulation Rate and Litter Size in Sheep*. Bachelor Thesis, Swedish University of Agricultural Sciences, Uppsala.

Kona, S. S. R., Praveen Chakravarthi, V., Siva Kumar, A. V. N., Srividya, D., Padmaja, K., and Rao, V. H. (2016). Quantitative expression patterns of GDF9 and BMP15 genes in sheep ovarian follicles grown in vivo or cultured in vitro. *Theriogenology* 85, 315–322. doi: 10.1016/j.theriogenology.2015.09.022

Kumar, S., Mishra, A. K., Kolte, A. P., Arora, A. L., Singh, D., and Singh, V. K. (2008). Effects of the Booroola (FecB) genotypes on growth performance, ewe's productivity efficiency and litter size in *Garole×Malpura* sheep. *Anim. Reprod. Sci.* 105, 319–331. doi: 10.1016/j.anireprosci.2007.03.012

Lassoued, N., Benkhil, Z., Woloszyn, F., Rejeb, A., Auina, M., Rekik, M., et al. (2017). FecX (Bar) a Novel BMP15 mutation responsible for prolificacy and female sterility in Tunisian Barbarine Sheep. *BMC Genet.* 18:43. doi: 10.1186/s12863-017-0510-x

Liao, W. X., Moore, R. K., and Shimasaki, S. (2004). Functional and molecular characterization of naturally occurring mutations in the oocyte-secreted factors bone morphogenetic protein-15 and growth and differentiation factor-9. *J. Biol. Chem.* 279, 17391–17396. doi: 10.1074/jbc.M401050200

Martin, P., Raoul, J., and Bodin, L. (2014). Effects of the FecL major gene in the lacaune meat sheep population. *Genet. Sel. Evol.* 46:48. doi: 10.1186/1297-9686-46-48

Martinez-Royo, A., Jurado, J. J., Smulders, J. P., Martí, J. I., Alabart, J. L., Roche, A., et al. (2008). A deletion in the bone morphogenetic protein 15 gene causes sterility and increased prolificacy in rasa aragonesa sheep. *Anim. Genet.* 39, 294–297. doi: 10.1111/j.1365-2052.2008.01707.x

McNatty, K. P., Juengel, J. L., Reader, K. L., Lun, S., Myllymaa, S., Lawrence, S. B., et al. (2005). Bone morphogenetic protein 15 and growth differentiation factor 9 co-operate to regulate granulosa cell function in ruminants. *Reproduction* 129, 481–487. doi: 10.1530/rep.1.00517

Menchaca, A., Cuadro, F., Dos Santos-Neto, P. C., Bosolasco, D., Barrera, N., de Brun, V., et al. (2018). Oocyte developmental competence is improved by relatively greater circulating progesterone concentrations during preovulatory follicular growth. *Anim. Reprod. Sci.* 195, 321–328. doi: 10.1016/j.anireprosci. 2018.06.010

Morón, F. J., de Castro, F., Royo, J. L., Montoro, L., Mira, E., Sáez, M. E., et al. (2006). Bone morphogenetic protein 15 (BMP15) alleles predict over-response to recombinant follicle stimulation hormone and iatrogenic ovarian hyperstimulation syndrome (OHSS). *Pharmacogenet. Genomics* 16, 485–495. doi: 10.1097/01.fpc.0000215073.44589.96

Mullen, M. P., Hanrahan, J. P., Howard, D. J., and Powell, R. (2013). Investigation of prolific sheep from UK and Ireland for evidence on origin of the mutations in BMP15 (FecX(G), FecX(B)) and GDF9 (FecG(H)) in Belclare and Cambridge sheep. *PLoS One* 8:e53172. doi: 10.1371/journal.pone.0053172

Palmer, J. S., Zhao, Z. Z., Hoekstra, C., Hayward, N. K., Webb, P. M., Whiteman, D. C., et al. (2006). Novel variants in growth differentiation factor 9 in mothers of dizygotic twins. *J. Clin. Endocrinol. Metab.* 91, 4713–4716. doi: 10.1210/jc.2006-0970

Persani, L., Rossetti, R., Di Pasquale, E., Cacciatore, C., and Fabre, S. (2014). The fundamental role of bone morphogenetic protein 15 in ovarian function and
its involvement in female fertility disorders. *Hum. Reprod. Update* 20, 869–883.
doi: 10.1093/humupd/dmu036

Pfaffl, M. W., Tichopad, A., Prgomet, C., and Neuvians, T. P. (2004). Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: bestkeeper–Excel-based tool using pair-wise correlations. *Biotechnol. Lett.* 26, 509–515. doi: 10.1023/b:bile.0000019559.84305.47

Poivey, J. P., Tiphine, L., Berny, F., and Jullien, E. (1995). Estimation of breeding values under BLUP animal model methodology for meat sheep. *Rencontres Recherches Ruminants* 2, 453–456.

Robinson, J. T., Thorvaldsdóttir, H., Winckler, W., Guttman, M., Lander, E. S., Getz, G., et al. (2011). Integrative genomics viewer. *Nat. Biotech.* 29, 24–26. doi: 10.1038/nbt.1754

Rochus, C. M., Tortereau, F., Plisson-Petit, F., Restoux, G., Moreno-Romieux, C., Tosser-Klopp, G., et al. (2018). Revealing the selection history of adaptive loci using genome-wide scans for selection: an example from domestic sheep. *BMC Genomics* 19:71. doi: 10.1186/s12864-018-4447-x

Rose, R. D., Gilchrist, R. B., Kelly, J. M., Thompson, J. G., and Sutton-McDowall, M. L. (2013). Regulation of sheep oocyte maturation using cAMP modulators. *Theriogenology* 79, 142–148. doi: 10.1016/j.theriogenology.2012.09.020

Scheet, P., and Stephens, M. (2006). A fast and flexible statistical model for large-scale population genotype data: applications to inferring missing genotypes and haplotypic phase. *Am. J. Hum. Genet.* 78, 629–644. doi: 10.1086/502802

Silva, B. D. M., Castro, E. A., Souza, C. J. H., Paiva, S. R., Sartori, R., Franco, M. M., et al. (2011). A new polymorphism in the growth and differentiation factor 9 (GDF9) gene is associated with increased ovulation rate and prolificacy in homozygous sheep. *Anim. Genet.* 42, 89–92. doi: 10.1111/j.1365-2052.2010.02078.x

Talebi, R., Ahmadi, A., Afraz, F., Sarry, J., Plisson-Petit, F., Genêt, C., et al. (2018a). Transcriptome analysis of ovine granulosa cells reveals differences between small antral follicles collected during the follicular and luteal phases. *Theriogenology* 108, 103–117. doi: 10.1016/j.theriogenology.2017.11.027

Talebi, R., Ahmadi, A., Afraz, F., Sarry, J., Woloszyn, F., and Fabre, S. (2018b). Detection of single nucleotide polymorphisms at major prolificacy genes in the Mehraban sheep and association with litter size. *Ann. Anim. Sci.* 18, 685–698. doi: 10.2478/aoas-2018-0014

Teixeira Filho, F. L., Baracat, E. C., Lee, T. H., Suh, C. S., Matsui, M., Chang, R. J., et al. (2002). Aberrant expression of growth differentiation factor-9 in oocytes of women with polycystic ovary syndrome. *J. Clin. Endocrinol. Metab.* 87, 1337–1344. doi: 10.1210/jcem.87.3.8316

Van der Auwera, G. A., Carneiro, M. O., Hartl, C., Poplin, R., del Angel, G., Levy-Moonshine, A., et al. (2013). From FastQ data to high confidence variant calls: the genome analysis toolkit best practices pipeline. *Curr. Protoc. Bioinformatics* 43, 11.10.1–11.10.33. doi: 10.1002/0471250953.bi1110s43

Wan, Q., Wang, Y., and Wang, H. (2015). Identification and Analysis of Regulatory Elements in Porcine Bone Morphogenetic Protein 15 Gene Promoter. *Int. J. Mol. Sci.* 16, 25759–25772. doi: 10.3390/ijms161025759

Wei, L.-N., Huang, R., Li, L.-L., Fang, C., Li, Y., and Liang, X.-Y. (2014). Reduced and delayed expression of GDF9 and BMP15 in ovarian tissues from women with polycystic ovary syndrome. *J. Assist. Reprod. Genet.* 31, 1483–1490. doi: 10.1007/s10815-014-0319-8

Wilson, T., Wu, X. Y., Juengel, J. L., Ross, I. K., Lumsden, J. M., Lord, E. A., et al. (2001). Highly prolific Booroola sheep have a mutation in the intracellular kinase domain of bone morphogenetic protein IB receptor (ALK-6) that is expressed in both oocytes and granulosa cells. *Biol. Reprod.* 64, 1225–1235. doi: 10.1095/biolreprod64.4.1225

Yan, C., Wang, P., DeMayo, J., DeMayo, F. J., Elvin, J. A., Carino, C., et al. (2001). Synergistic roles of bone morphogenetic protein 15 and growth differentiation factor 9 in ovarian function. *Mol. Endocrinol.* 15, 854–866. doi: 10.1210/mend.15.6.0662

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Chantepie, Bodin, Sarry, Woloszyn, Plisson-Petit, Ruesche, Drouilhet and Fabre. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.