# Spatio-temporal gene expression profiling during in vivo early ovarian folliculogenesis: integrated transcriptomic study and molecular signature of early follicular growth

Agnès Bonnet, Bertrand Servin, Philippe Mulsant, Béatrice Mandon-Pepin

► To cite this version:

Agnès Bonnet, Bertrand Servin, Philippe Mulsant, Béatrice Mandon-Pepin. Spatio-temporal gene expression profiling during in vivo early ovarian folliculogenesis: integrated transcriptomic study and molecular signature of early follicular growth. PLoS ONE, 2015, 5 (November), pp.1-25. 10.1371/journal.pone.0141482. [https://doi.org/10.1371/journal.pone.0141482](https://doi.org/10.1371/journal.pone.0141482)

HAL Id: hal-02634306

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

Submitted on 27 May 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
RESEARCH ARTICLE

# Spatio-Temporal Gene Expression Profiling during *In Vivo* Early Ovarian Folliculogenesis: Integrated Transcriptomic Study and Molecular Signature of Early Follicular Growth

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

CrossMark
click for updates

Agnes Bonnet¹,²,³*, Bertrand Servin¹,²,³, Philippe Mulsant¹,²,³, Beatrice Mandon-Pepin⁴

1 INRA, UMR 1388 GenPhySE (Génétique, Physiologie et Systèmes d'Elevage), F-31326 Castanet-Tolosan, France, 2 Université de Toulouse, INP, ENSAT, GenPhySE (Génétique, Physiologie et Systèmes d'Elevage), F-31326 Castanet-Tolosan, France, 3 Université de Toulouse, INP, ENVT, GenPhySE (Génétique, Physiologie et Systèmes d'Elevage), F-31076 Toulouse, France, 4 INRA, UMR1198 Biologie du Développement et de la Reproduction, F-78350 Jouy-en-Josas, France

* agnes.bonnet@toulouse.inra.fr

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

OPEN ACCESS

**Citation:** Bonnet A, Servin B, Mulsant P, Mandon-Pepin B (2015) Spatio-Temporal Gene Expression Profiling during *In Vivo* Early Ovarian Folliculogenesis: Integrated Transcriptomic Study and Molecular Signature of Early Follicular Growth. PLoS ONE 10(11): e0141482. doi:10.1371/journal.pone.0141482

**Editor:** Wei Yan, University of Nevada School of Medicine, UNITED STATES

**Received:** March 25, 2015

**Accepted:** October 8, 2015

**Published:** November 5, 2015

**Copyright:** © 2015 Bonnet et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All the raw RNA-seq data have been deposited in EMBL-EBI ArrayExpress http://www.ebi.ac.uk/arrayexpress/ under accession number E-MTAB-1587. All the relevant analyses are within the paper and its Supporting Information files.

**Funding:** This work "BioFoll: transcrits caractéristiques ou « biomarqueurs » de la croissance folliculaire basale chez la brebis" was supported by "Bioressources 2010" grant under INRA programmes. The funding was received by AB. The

## Abstract

### Background

The successful achievement of early ovarian folliculogenesis is important for fertility and reproductive life span. This complex biological process requires the appropriate expression of numerous genes at each developmental stage, in each follicular compartment. Relatively little is known at present about the molecular mechanisms that drive this process, and most gene expression studies have been performed in rodents and without considering the different follicular compartments.

### Results

We used RNA-seq technology to explore the sheep transcriptome during early ovarian follicular development in the two main compartments: oocytes and granulosa cells. We documented the differential expression of 3,015 genes during this phase and described the gene expression dynamic specific to these compartments. We showed that important steps occurred during primary/secondary transition in sheep. We also described the *in vivo* molecular course of a number of pathways. In oocytes, these pathways documented the chronology of the acquisition of meiotic competence, migration and cellular organization, while in granulosa cells they concerned adhesion, the formation of cytoplasmic projections and steroid synthesis. This study proposes the involvement in this process of several members of the integrin and BMP families. The expression of genes such as *Kruppel-like factor 9* (KLF9) and *BMP binding endothelial regulator (BMPER)* was highlighted for the first time during early follicular development, and their proteins were also predicted to be involved in gene regulation. Finally, we selected a data set of 24 biomarkers that enabled the
funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing Interests:** The authors have declared that no competing interests exist.

discrimination of early follicular stages and thus offer a molecular signature of early follicular growth. This set of biomarkers includes known genes such as *SPO11* meiotic protein cova- lently bound to DSB (*SPO11*), bone morphogenetic protein 15 (*BMP15*) and *WEE1* homo- log 2 (*S. pombe*) (*WEE2*) which play critical roles in follicular development but other biomarkers are also likely to play significant roles in this process.

## Conclusions

To our knowledge, this is the first *in vivo* spatio-temporal exploration of transcriptomes derived from early follicles in sheep.

## Introduction

In large mammalian species at birth, the ovaries contain a large and fixed reserve of non-grow- ing primordial follicles (oocytes surrounded by flattened pre-granulosa cells). Most of these fol- licles remain in the resting state until either their degeneration or their activation and growth towards the primary, secondary and tertiary stages (with an antral cavity). The gradual exit of primordial follicles starts shortly after formation of the primordial follicle pool and continues throughout the reproductive years [1]. This early follicular development is therefore important as it regulates the size of the remaining stock of primordial follicles and their fate, which in turn affects fertility and the reproductive life span.

Early follicular development is accompanied by an increase in oocyte diameter, a progres- sive acquisition of competence [2] and proliferation of granulosa cells (GCs). Early follicular development requires the appropriate expression of numerous genes at different developmen- tal stages and orchestrated communication between the two main compartments (oocytes and granulosa cells) [3, 4]. These compartments regulate follicle growth in an autocrine and para- crine manner via secreted factors and direct gap junctional communication.

Therefore, any integrated study of folliculogenesis must include changes in gene expression in all cell types, but such studies are faced with two major problems: restricted access to isolated stages and the limited supply of RNA. Indeed, the presence of all follicular stages in the ovary and the small size of preantral follicles, render the isolation of each early follicular stage even more problematic. As a result, relatively little is known about the key molecular mechanisms that underlie the complex biological processes of *in vivo* early folliculogenesis.

Until now, some transcriptomic studies focused on primordial-primary transition from neonatal rodent ovaries on day 0-2 and day 4 that were respectively enriched in primordial and primary follicles [5, 6], or from pure populations of oocytes from human primordial, inter- mediate and primary follicles [7]. Only two mouse studies explored the breadth of early follicu- lar development from oocytes and whole follicles, respectively [8, 9]. They showed that the principal change in expression in mouse species occurs at oocyte primordial-primary transition [8]. They revealed an over-representation of differentially expressed genes involved in protein synthesis and the cell cycle, in particular M-phase, throughout the early growth phase [8]. The differential expression of oocyte-secreted ligands involved in NOTCH, SHH, EGF, TGFβ, PDGF signaling pathways [8, 9] also suggests that multiple signaling pathways operate during follicular development. Nevertheless, these results obtained in poly-ovulating species cannot be fully exploited in mono-ovulating species (sheep, human), so recent studies addressed this issue in the latter. A first study focused on secondary and small antrum follicles in goat species and highlighted three main metabolic pathways: lipid, cell death and hematological pathways
[10]. However, our previous studies had been the first to describe differential gene expression between two isolated compartments, granulosa cells and oocytes, in sheep ovarian follicles [11, 12]. We identified enriched functional categories that reflect two distinct cell fates. We revealed the involvement of granulosa cell pathways such as SHH, WNT and RHO GTPase. Other signaling pathways such as VEGF, NOTCH and IGF1 suggested the existence of complex cell-cell interactions. The expression and the specificity of the two cell types are now characterized in sheep species [11].

Nevertheless, until now no studies had been able to describe the gene expression changes that occur during early follicular growth and in each compartment of the follicle. Consequently, the aim of the present study was to: 1- focus and decipher gene expression and functions that might be associated with early folliculogenesis from a developmental point of view (from primordial, primary, secondary and small antrum follicles) in sheep species, and 2- produce a set of biomarkers that are able to predict early follicular growth and could be used as *in vivo* signatures of early folliculogenesis.

For this purpose, we used the previous dataset of 15,349 gene expressions obtained by RNA-seq technology [12] in order to analyze the spatio-temporal gene expression. We were able to show that major changes in gene expression occur during the primary-secondary transition. We compiled a detailed description of the molecular mechanisms underlying early folliculogenesis in sheep which are probably transposable to other large mammals including humans.

# Materials and Methods

## RNA-seq experiment

This study was carried out in strict accordance with the institutional guidelines for research studies. The protocol was approved by the prefecture de la Haute Garonne (animal experimentation authorization no 31-297).

We briefly described the dataset previously published [11]:

- The lambs were euthanized at birth (<1 day). The ovaries were removed and embedded in O.C.T. embedding matrix, frozen in liquid nitrogen and stored at -80°C until use. Frozen sections were fixed and stained before laser capture microdissection [11].

- The follicular compartments (i.e. granulosa cells (GCs)) and oocytes (O) were selected for each follicular stage as a function of follicular diameter (primordial (PD<35µM), primary (PM: 35-50µM), secondary (SC: 60-120µM) follicles and small antral follicles (SA: 250-500µM), captured by laser capture microdissection and extracted using the PicoPure RNA Isolation Kit. Three/four independent replicates were obtained per sample.

- Thirteen tissue samples (fetal ovary, pituitary gland, hypothalamus, muscle, skin, heart, lung, intestine, stomach, kidney, spleen, liver, and theca) were collected in triplicate (from three separate animals) at the local slaughterhouse for RNA extraction. The 12 total RNA tissue samples (not including fetal ovary) were pooled in three separate groups.

- A total of 31 RNA LCM-derived (3 oocyte samples from primordial follicles (PDO), 4 oocyte samples from primary follicles (PMO), 4 oocyte samples from secondary follicles (SCO), 4 oocyte samples from small antral follicles (SAO), 4 granulosa samples from primordial follicles (PDG), 4 granulosa samples from primary follicles (PMG), 4 granulosa samples from secondary follicles (SCG) and 4 granulosa samples from small antral follicles (SAG)) and three multi-tissue RNA samples (3 independent pools of 12 total RNA tissue samples) were subjected to two rounds of T7 linear amplification.
- Sequencing was performed on an Illumina HiSeq2000 using the Illumina TruSeq SBS kit v2 (209 cycles including the index) to obtain paired-end reads (2x100 base pairs (bp)).

- The reads were mapped with bwa aln to the sheep genome sequence [13] (CSIRO Oarv2.0 released March 2011: http://www.livestockgenomics.csiro.au/sheep/). Overlapping reads were merged to produce genomic fragments.

The gene annotation of genomic fragments was performed by alignment with the Bos taurus UMD3.1 genome using Blat and filtered to conserve only the best genomic fragment for each gene.

## Statistical analysis

The significance of differential gene expression was determined using the DESeq package [14] of R software for multi-factorial design (R 2.14.0; DESeq release 1.6.1). Our experimental design included two factors: “stage” (with four levels: PD, PM, SC and SA) and “compartment” (two levels: O, GCs) [12]. To establish the differential expression profile throughout early follicular development we identified the genes whose intra-compartment variance was differentially expressed during at least one of the stages (nbinom GLM Test). For this, we specified two models: the stage model (count ~ stage) and the reduced model (count ~ 1) and applied them to each compartment dataset. For each compartment, the expression value of each gene was then compared between each follicular stage using pairwise comparison (nbinom Test). Last, differential gene expression during early follicular development was selected with a FDR <5% for the global intra-compartment effect (nbinom GLM Test) and pval <1% for the pairwise comparison (nbinom Test) with a fold change >2.

## Quantitative real-time PCR analysis of gene expression

Gene primer designs were based on the RNA-seq sequences as previously described [12] (S1 Table).

The LCM-derived aRNA samples were reverse transcribed as previously described [11]. The reactions were completed to 50µl and diluted at 1/20 prior to PCR. The assay for each gene consisted of four replicates per condition (except for PDO = 3) and negative controls. Gene expression was analyzed using 96.96 Dynamic Array™ IFCs and the BioMark™ HD System from Fluidigm, as previously described [11].

The PCR amplification efficiency of genes was determined from a foetal ovary cDNA pool (11 ng cDNA) that was serial diluted (1, 1:3; 1:3; 1:2; 1:2).

After determining the threshold cycle (Ct) for each LCM-derived aRNA sample, the PFAFFL method was applied to calculate the relative expression of each gene [15]. This relative expression was normalized using the corresponding geometric average of three reference genes with geNorm v3.4 [16]: *β-actin*, transmembrane p24 trafficking protein 4 (*TMED4*) and ribosomal protein L19 (*RPL19*) genes that were respectively slightly, moderately and highly expressed and not regulated during follicular development or in the compartment.

The significance of the relative expression data on genes involved in canonical pathways was tested using the one-way ANOVA model of the R statistical software system (the Comprehensive R Archive National, http://www.cran.r-project.org) after logarithm transformation of the data. For each gene, an ANOVA model was fitted using two factors: “stage” (4 levels) and “compartment” (2 levels) and their interactions. A backward variable selection procedure was applied, as previously described [11].
## Biological trends

Ingenuity® Pathway Analysis software (IPA; http://www.ingenuity.com) was used to examine the functional enrichment in differentially expressed genes. The significance of the association between the differential gene lists and functional categories was determined by a P value calculated using Fisher's exact test corrected for multiple testing (FDR <0.05 (Benjamini-Hochberg test)).

Pathway based analysis was performed using Webgestalt (http://bioinfo.vanderbilt.edu/webgestalt/ Kegg, Wikipathway and Pathwaycommons databases) and GeneCodis (http://genecodis.cnb.csic.es/ Kegg pathway and Panther databases) software programs. Significantly enriched pathways were identified using a hypergeometric test followed by a multiple testing correction of p-values using the Benjamini-Hochberg test (FDR <0.05). Then, to ensure that pathways obtained from computational biology were coherent and gave a good interpretation of experimental results [14, 17], only Webgestalt significant pathways identified using the two softwares (Webgestalt and GeneCodis) or from two databases, were conserved (Kegg/Wiki-pathway/ Pathwaycommons databases).

IPA downstream Effects Analysis was then used to identify the expected effects of our observed gene expression changes on the functions (expected to increase or decrease). This analysis examined genes in our dataset that were known to affect functions, compared the directions of change of the genes with those expected based on findings in the literature, and then issued a prediction for each function based on the direction of change (in the differential gene expression list).

In the same way, IPA Upstream Regulator analysis was used to understand the reasons for the gene expression changes observed. It identified the cascade of upstream transcriptional regulators that might explain the gene expression changes observed in the dataset and predicted the activity of their encoded protein (cf. Methods section). Taking account of the gene expression profile (expressed and/or differentially expressed) for each predicted upstream regulator provided more evidence for the biological mechanism. The expression and/or differential expression of the transcriptional regulator in the dataset produced more evidence for a particular biological mechanism. These analyses were described by Bonnet et al. [12].

## Biomarker selection

Biomarkers of the compartments and stages of early folliculogenesis were determined using the DESeq package [14] of R software for multi-factorial design (R 2.14.0; DESeq release 1.6.1). The genes mostly expressed in a specific cell type and stage (PDO, PMO, SCO, SAO, PDG, PMG, SCG, SAG) were selected using pairwise comparisons (nbinom Test) with other samples, including multi-tissue samples (FDR<5% and a fold change >3 for granulosa samples and a fold change >10 for oocytes samples (for example, PDO / PMO + SCO + SAO + PDG + PMG + SCG+SAG + MT)). Theca cells were included in multi-tissue samples to ensure the selection of specific expression in oocytes or granulosa cells and to prevent contamination.

The accuracy of biomarker selection was checked using qRT-PCR and the significance of differential gene expression was evaluated using Student's test after a fourth square transformation.

PLS regression was used to classify and discriminate the samples (PLS-DA). PLS-DA was implemented using the mixOmics package under R software.

## Linear mixed models for expression data

The linear mixed models were described in S1 Text. Briefly:
**Logistic regression model for the presence/absence of expression.** We first of all estimated a model based on the binary response and indicating the expression or lack of expression, with respect to the gene, developmental stage, cell type and replicate. This was modeled with a hierarchical logistic regression, including correlation parameters between random effects of different stages. All random effects were found to be highly significant ($p<10^{-7}$).

**Linear regression on the quantitative level of expression.** Secondly, we estimated a model on the level of expression for gene, developmental stage, cell type and replicate using a fourth-root transformation of the raw expression. In cases where no expression was found in any replicate, we included a single observation with null expression in the data.

**Equations for the prediction of developmental stage given a vector of expression data.** Given the above models estimated on the RT-PCR data, we derived equations that could predict the stage of the follicle, given the cell type and a vector of expression data on the genes. Specifically, we were dealing with a situation where we were given a vector of new observations (binary response and expression level), for a given cell type, and wanted to predict the stage reached by the corresponding follicle.

**Predictive ability of the models.** To test the predictive ability of the models, we implemented the predictive equations using a customized program and generated new vectors for observations through re-sampling: for a given cell type and for each stage, we generated a new data vector with one observation for each gene by sampling one value taken at random from the observed data. Thus each vector was a new set of observations drawn from the original data. We performed this re-sampling method 100 times for each stage. We then applied the predictive equations to each of the 100 vectors and recorded the posterior probabilities of each stage.

# Results and Discussion

## Global differential gene expression during early follicular development

This study was carried out using the gene profiles of 15,349 genes generated by RNA-seq experiment (http://www.ebi.ac.uk/arrayexpress/ under accession number E-MTAB-1587) [12]. This experiment used a laser capture microdissection method to split and produce granulosa cells or oocyte transcriptomes for each stage of early follicular development: primordial (PD), primary (PM), secondary (SC) and small antral (SA) follicle stages. The experiment included three/four biological replicates per sample. Previous study applied a generalized linear model (DESeq GLM) to this dataset to explore the differential expression and the specificity between each compartment (oocyte, granulosa cells). We also observed 19.6% of genes that were significantly differentially expressed (DEG) during early development (3,015 genes, FDR<5%) [12].

In this study, we specifically investigated the 3,015 DEG (S2 Table) in order to clarify the dynamics of transcription during early development and we observed more marked changes in gene expression in oocytes than in granulosa cells (2173/1192). A majority of the differentially expressed genes were down-regulated during early folliculogenesis in granulosa cells (66%), while the numbers of up- and down-regulated genes were similar in oocytes (Table 1).

An unsupervised hierarchical clustering was performed to explore differential expression and confirmed the relevance of these data (Fig 1). The branching distance of the first level of each dendrogram revealed a marked dissimilarity between primordial/primary follicles and secondary/small antrum follicles. The distance branching of the second and third levels of the dendrogram highlighted a clear separation between the other follicular stages. These clusterings depict the dynamics of expression during early follicular development and suggested important differential expression between primary and secondary stages.
<table><caption>Table 1. Summary of differential gene expressions.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Number of differential genes</th></tr><tr><th>Primary stage</th><th>Secondary stage</th><th>Small antrum stage</th></tr></thead><tbody><tr><td>oocyte</td><td>173</td><td>1005</td><td>2074</td></tr><tr><td>granulosa cells</td><td>171</td><td>732</td><td>1074</td></tr></tbody></table>

doi:10.1371/journal.pone.0141482.t001

In term of differential gene expression between two consecutive stages (S1A and S1B Fig), a weak differential expression could be seen during primordial/primary transition. The differences in gene expression between stages mainly started during primary-secondary transition, when 74% and 67.4% of differential genes displayed decreased expression in oocytes and granulosa cells, respectively. The number of differentially expressed genes was then maintained during oocyte secondary/small antrum transition (476 genes in the primary/secondary stages and 494 genes in the secondary/small antrum stages) whereas few differential gene expressions were found in granulosa cells during secondary/small antrum transition: 55 genes (S1A and S1B Fig).

These results therefore described a difference in expression dynamics between the two compartments during early development, with 1- important and constant changes in gene expression in oocytes that were in line with the extremely rapid growth of oocytes prior to the formation of antrum [18], and 2- fewer changes affecting granulosa gene expression, which was in line with the fact that most granulosa cell proliferation and differentiation occurs later during the antral phase when the oocytes have almost stopped growing [18]. Lastly, contrary to mouse species in which primordial/primary transition has been identified as the major transition [8], in sheep species we were able to highlight major changes of expression during primary/secondary transition (Fig 1 and S1 Fig).

## Validation of differential expression

Consistent with previous reports, our list of differentially expressed genes (3,015 genes) included genes known for their differential expression during early follicular development (for

<figure><img src="image_4.png" /><figcaption>**Fig 1. Heatmap display of unsupervised hierarchical clustering of all differential genes during early folliculogenesis.** This Figure shows two unsupervised hierarchical clusterings from oocyte (A) and granulosa cells (B) differential gene expressions. The first level of the dendrogram classified the primordial/ primary follicles and the secondary/small antrum follicles as two distinct clusters. The other levels revealed a lower dissimilarity between primordial/primary than small antrum/secondary stages. The genes are shown in lines and the relative expression means by follicle stages are shown in columns. Red, black and green represent the levels of expression: up, mean, and down, respectively. PDO: oocyte from primordial follicles, PMO: oocyte from primary follicles, SCO: oocyte from secondary follicles, SAO: oocyte from small antrum follicles, PDG: granulosa cells from primordial follicles, PMG: granulosa cells from primary follicles, SCG: granulosa cells from secondary follicles, SAG: granulosa cells from small antrum follicles.</figcaption></figure>

doi:10.1371/journal.pone.0141482.g001
example BMP15, *growth differentiation factor 9* (GDF9), spermatogenesis and *oogenesis specific basic helix-loop-helix 2* (SOHLH2), *NLR family pyrin domain containing 5* (NLRP5)). In addition, the expression profiles of 19 genes of interest involved in enriched canonical pathways (BMP, IGF1, Gap junction, WNT, FGF, PI3K, RAR activation, apoptosis.) were monitored using qRT-PCR, and statistical analysis confirmed the DE observed in the RNA-seq dataset for 14 of them (S3 Table).

## Global outlook relative to gene functions and pathways

Six global lists of differentially expressed genes (taking account of pairwise comparisons) were created to explore significant functions and pathways involved in early folliculogenesis (oocyte lists (PMO/PDO, SCO/PDO, SAO/PDO) and granulosa cell lists (PMG/PDG, SCG/PDG, SAG/ PDG)).

**Biological functions.** The biological functions associated with early follicular development were investigated using Ingenuity Pathway Analysis (IPA). Ten functional categories and nine main cellular function categories were identified as being significantly enriched in differentially expressed genes (FDR<0.05) and are presented respectively in S2 Fig and Fig 2.

At the primary stage, the DNA replication, recombination and repair category was specifically enriched in differentially expressed genes. These genes are essential to protecting genetic information during meiosis. Their corresponding proteins repair copy errors and block the resumption of meiosis. In the granulosa cell compartment, we observed the expression of several genes which were also involved in this functional category (DNA replication, recombination and repair category) but they did not change during follicular development. This granulosa cell gene expression stability was in agreement with a slow follicular grow process (130 days in sheep species) and a weak rise in the granulosa cell mitotic index up to early antral follicles in either sheep (from 60 µm (secondary follicles: 831 h doubling time) to 280 µm of follicular diameter (early small antrum follicles: 244 h doubling time)) [19] or human species [20]. Intense granulosa cell proliferation mainly occurred later during antral stage (from 280 µm follicular diameter up to 850 µm in sheep species) and was supported by a rapid rise in the mitotic index [19].

At the secondary stage, IPA analysis highlighted major functional change and two biological functions relative to protein metabolism (degradation, synthesis) as being mainly over-represented in the secondary stages of granulosa cells but not enriched in differentially expressed genes during oocyte development.

**Pathways.** Pathway based analyses were investigated using Webgestalt software and revealed significant pathway enrichment in the lists of differentially expressed genes (FDR<0.05; S4 Table). These pathways encompassed cell processing, signaling pathways and those involved in cell communication that were slightly enriched in primary follicles and mainly enriched from the secondary stage when compared to primordial stage (Fig 3). However, some pathways were only significantly enriched at the small antrum stage.

At the primary stage, we noted a slight enrichment in regulated genes involved in cell dialog such as cytokine-cytokine receptors and integrin cell surface interactions, chemokine signaling, focal adhesion and TGFB and WNT signaling. The PPAR signaling pathway, which is known to be involved in lipid metabolism, was specifically enriched in granulosa cells.

At the secondary stage, we underlined the enrichment of genes involved in oocyte cell cycle and meiosis (*cyclin-dependent kinase 1* (CDK1), *cyclin B1* (CCNB1), *Moloney sarcoma oncost gene* (MOS), *TTK protein kinase* (TTK), etc...) (Fig 3). We also observed pathways that were related to physical and molecular communication such as endocytosis, tight and gap junctions, ECM receptor interaction and neuronal signal transmission (axonal guidance). Fig 3 shows the
<figure><img src="image_5.png" /><figcaption>**Fig 2. Cellular function enrichment.** Functional enrichment analysis for sets of differentially expressed genes during early follicular development was performed *in silico* for each compartment using Ingenuity Pathway Analysis (IPA) software. Statistical significance was determined by a P value calculated using Fisher's exact test corrected for multiple testing correction (Benjamini-Hochberg test, FDR <0.05). This analysis identified nine main cellular function categories as being significantly enriched in differentially expressed genes. Bar colors correspond to enriched functions at primary stage (compared to primordial stage), secondary stage (compared to primary and primordial stages), and small antrum stage (compared to secondary, primary and primordial stages). The X axis corresponds to the level of significance of the function: -log(B-H p value). Granulosa cell functions are colored in red and oocyte functions are colored in blue. Numbers correspond to the number of focus genes contributing to the functions. Follicular stage abbreviations are described in Fig 1.</figcaption></figure>

doi:10.1371/journal.pone.0141482.g002

enrichment of additional signaling pathways such as the BMP receptor signaling network, mTOR, calcium and Jak-STAT signaling pathways.

At the small antrum stage, we identified an enrichment of genes involved in the oocyte DNA damage response and granulosa cell steroid hormone signaling processes (Fig 3). Finally, members of the PDGF receptor family (such as *platelet derived growth factor D* (PDGFD), *platelet-derived growth factor receptors* (PDGFRA and PDGFRB), etc...) and members of FGF

<figure><img src="image_6.png" /><figcaption>**Fig 3. Significantly enriched pathways highlighted in early folliculogenesis.** Pathway enrichment analysis for sets of differentially expressed genes during early follicular development was performed using Webgestalt software. Statistical significance was determined by multiple testing correction of p-values using Benjamini-Hochberg test (FDR <0.05). Thirty one enriched pathways are presented. Each graph corresponds to the percentage of gene enrichment in the pathway (%) according to different stages of follicles (PM (primary), SC (secondary) and SA (small antrum)). Granulosa canonical pathways are colored in red and oocyte canonical pathways are colored in blue.</figcaption></figure>

doi:10.1371/journal.pone.0141482.g003
signaling pathways (*fibroblast growth factor 1* (FGF1), *neural cell adhesion molecule 1* (NCAM1), *syndecan 2* (SDC2), etc...) were regulated in oocytes of small antrum follicles. Members of the G protein signaling pathway were regulated in granulosa cells while the insulin, NOTCH, and MAPK signaling pathways were enriched in differential genes in both the oocyte and granulosa cells compartments (S4 Table).

**Downstream effects.** In addition to the above studies, we used IPA to predict the downstream effect of differentially expressed genes on the functions (expected to be activated or inhibited). Apart from the categories linked to cell growth, this analysis showed that numerous categories were affected (S5 Table) which participated in cell communication and organization, cell growth and lipid metabolism.

Together the data on the principal processes in play, these analytical findings now offer an overview of early follicular development.

## Biological interpretation

At the end of the preantral stage, a large quantity of mRNA is stocked in oocytes which then become competent to resume meiosis. Moreover, granulosa cells may synthesize estrogens and respond to gonadotropin hormones.

The different IPA results (functions (Fig 2; S2 Fig), pathways (S4 Table), downstream effects (S5 Table)) were compiled to describe the molecular processes involved in achieving this status.

**Chronology of meiotic maturation.** Meiotic maturation is essential for the correct formation of metaphase II oocytes, and is mainly acquired during early folliculogenesis [21]. This complex process depends on the interplay in time and space between the cytoskeleton (spindle stability) and the cell cycle machinery, and on the coordinated exchange of signals with somatic cells. During this study, we identified a chronology in the expression of genes involved in the regulatory network governing the maintenance of meiotic arrest and the acquisition of meiotic competence (Fig 4).

During the primordial to primary transition, we observed an increase in the expression of genes including *SET domain and mariner transposase fusion gene* (*SETMAR*), *TTK*, *discs, large* (*Drosophila*) *homolog-associated protein 5* (*DLGAP5*) and *kinesin family member 18A* (*KIF18A*), etc., which are involved in the DNA replication, recombination and repair (Fig 4A). Thus, for the first time, we are able to identify an ovarian expression of the *SETMAR* gene (Fig 4A) that was shown to promote double strand break (DSB) [22] and non-homologous end joining (NHEJ) [23] repairs in human lineage cells. Its over expression in oocyte of primary follicles, by comparison with primordial follicles, is in line with the role of the protein to maintain genomic stability. TTK, DLGAP5 and KIF18A regulate the cell cycle by controlling chromosome alignment, spindle assembly and attachment. TTK (or MPS1) is a spindle assembly checkpoint protein that is necessary for chromosome segregation during oocyte meiosis I [24]. DLGAP5 stabilizes kinetochore fibers and promotes chromosome alignment on the spindle. Kinesins play essential roles in assembling spindles, separating centrosomes and attaching chromosomes to spindles [25] and could cause female sterility [26]. In Hela cells or regenerating mouse liver, the *DLGAP5* mRNA levels were found to be tightly regulated during cell cycle progression. In particular, the level of expression was elevated during the G2/M phase [27]. Consequently, the over-expression of *DLGAP5* mRNA in ovine oocytes from primary follicles accorded with their cell cycle arrest at the G2/M phase. Its protein also regulates KIF18A localization [28] whose gene was highlighted for the first time as differentially expressed during primordial/primary transition (3.62 fold change).

At the secondary stage, we identified an increased expression of genes such as *CDK1* and *CCNB1*, the constituents of oocyte maturation promoting factor (MPF), as well as the meiosis
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Meiotic maturation gene expression as a function of oocyte developmental stage. (A) Differential gene expression at primary stage, (B) differential gene expression at secondary stage, (C) differential gene expression at small antrum stage. The Y axis corresponds to normalized counts from RNA-seq data. \*: FDR <0.05, pairwise comparison pval <0.01.</figcaption></figure>

doi:10.1371/journal.pone.0141482.g004

inhibitor *WEE2* [29] (S4 Table, cell cycle; Fig 4B). MFP governs the G2/M transition of the cell cycle and must remain inactive to ensure the arrest of meiosis. The ability to resume meiosis is acquired when an oocyte obtains a sufficient amount of MPF. Our findings were in agreement with the initiation of acquisition of meiotic competence already described in mice at the secondary stage [8].

At the small antral stage, the acquisition of meiotic competence progressed (cell cycle, cell cycle checkpoints and oocyte meiosis pathways; S4 Table). The progression in each phase of the cell cycle is governed by cyclin-dependent kinases (Cdks) within Cyclin-Cdk complexes. These Cdks are negatively regulated by small polypeptides from the Cip/Kip and INK4 families such as *cyclin-dependent kinase inhibitor 1A* (CDKN1A/p21/cip1) and *cyclin-dependent kinase inhibitor 2C* (CDKN2C/p18INK4c), respectively [30]. Our data point out this cell cycle regulation with for example the over-expression of *S-phase kinase-associated protein 2*, *E3 ubiquitin protein ligase (SKP2)* that is involved in the degradation of CDKN1A [31], and then the down expression of *CDKN1A* (S4 Table, cell cycle; Fig 4C). CDKN2C gene also displayed an interesting temporal expression profile with decreased expression during primary/secondary transition and increased expression during secondary/small antrum transition. This gene profile illustrated the fine regulation of oocyte meiosis arrest during oocyte development. This profile accorded with the fact that the *CDKN2C* gene may function in a cell type-specific [30] and differentiation-specific manner [32]. Last we observed the over-expression of genes necessary to stabilize mouse oocyte meiotic arrest such as *cell division cycle 25A* (CDC25A) [33] or to maintain the fidelity of chromosome segregation (*shugoshin-like 1* (*S. pombe*) (*SGOL1*), *structural maintenance of chromosomes 3* (SMC3); S4 Table, oocyte meiosis).

**Cell-cell communication.** Follicular growth is dependent not only on close interactions between oocytes and granulosa cells but also between granulosa cells themselves. These interactions include molecular dialog and physical communications [12].

During oocyte development, the migration and cellular organization categories were affected by changes in gene expression.
At the primary stage, IPA highlighted genes involved in cell migration that might be associ-
ated with follicle activation in the growing phase (S5 Table; migration of endothelial cells:
increased expression of *adrenomedullin* (ADM), *calcitonin-related polypeptide alpha* (CALCA),
*thyroglobulin* (TG), *gastrin-releasing peptide* (GRP) genes and decreased expression of inhibitor
of DNA binding 3, dominant negative helix-loop-helix protein (ID3), *protein kinase D1*
(PRKD1), *secreted frizzled-related protein 4* (SFRP4) genes). For example, ADM [34] and
CALCA [35] genes are known to regulate cell proliferation and migration in vascular endothe-
lial and vascular smooth muscle cells.

At secondary and small antrum stages, the number of differentially expressed genes involved in cell migration increased (S5 Table; migration of cells: 182 differential expression genes (DEG)). For example, genes related to cell migration signaling pathways such as integrins (integrin alpha M (ITGAM), integrin beta 8 (ITGB8), caveolin 1 (CAV1), etc.), WNT (trans-forming growth factor, beta receptor II (TGFBR2), SFRP4, wingless-type MMTV integration site family member 5A (WNT5A), etc.) or Rho (Fig 3; S4 Table) were mainly down-regulated.

During follicular development and at the secondary stage, granulosa cells exhibited an inter-
esting subset of differentially expressed genes that was concerned with the adhesion (S5 Table;
adhesion of blood cells: 25 DEG) and the formation of cytoplasmic projections (S5 Table; for-
mation of neurites: 12 DEG). Lastly, a number of regulated genes encode for proteins display-
ing cell contact (neural precursor cell expressed, developmentally down-regulated 9 (NEDD9),
secreted phosphoprotein 1 (SPP1), integrin alpha 5–6 (ITGA5-6); S5 Table; adhesion of blood
cells) and chemokine properties as kit ligand (KITL) (S5 Table; cytokine-cytokine receptor
interaction).

Last, these results offer further evidence to reinforce the role of WNT and integrin signaling pathways in cell-cell communication during early follicular development (Fig 3).

The expression of antagonists that regulate the WNT pathway, such as dickkopf-like 1 (DKKL1), dickkopf WNT signaling pathway inhibitor 1 (DKK1) and secreted frizzled-related protein2 (SFRP2), increased during preantral folliculogenesis. We also identified the granulosa cell over expression of cadherin, EGF LAG seven-pass G-type receptor 2 (CELSR2) from primary follicles and two additional regulating genes from secondary follicles as catenin alpha 2 (CTNNA2), low density lipoprotein receptor-related protein 8 (LRP8) (S5 Table). CELSR2 is an atypical cadherin, that is able to transduce signals from the non canonical WNT/PCP signaling pathway [37] and plays an important role the control of tissue polarity in planar cells [38]. CTNNA2 is a catenin isoform that links cadherin adhesion receptors. LRP8 is a positive regulator of Wnt/β-catenin signaling that increases Wnt-induced transcriptional responses [36]. Gene expression profiles are presented Fig 5.

In addition to the regulated integrin family genes identified during oocyte development (ITGAM, ITGB8), we noted the regulation of other genes of this family during follicular progression such as *ITGA5-6* and *integrin beta-like 1 (ITGBL1)* (Fig 6). *ITGA5* expression has never been described in the ovary. *Itga6* is expressed in the granulosa cells of early follicles in the mouse [39]. In marmoset species, atretic tertiary follicular status is associated with a lack of *ITGA6* and *ITGB1* expression [40]. Our study showed that the expression of *ITGA6* increased at the secondary stage whereas that of *ITGA5* decreased. In the same way, *ITGBL1* increased. In accordance with marmoset data, *ITGA6* expression profile suggests that *ITGA6* may promote cell contact in order to preserve follicle health during early folliculogenesis in sheep.

**Granulosa cell steroid synthesis.** At the end of early follicular development, granulosa cells may be capable of producing steroids. Pig preantral follicles produced hormones such as estradiol and progesterone [41] that originate from cholesterol and are required for the development and maturation of follicles [42]. We demonstrated a gradual increase during follicular development in the expression of genes already known to participate in lipid and steroid
<figure><img src="image_8.png" /><figcaption>**Fig 5. Gene expression profiles of members of WNT signaling pathway during early folliculogenesis.** Significant differences in relative gene expression during early folliculogenesis of members of Wnt signaling pathway (S2 Table). The Y axis corresponds to normalized counts from RNA-seq data. Granulosa cells data are colored in red and oocyte data are colored in blue. \*: FDR <0.05, pairwise comparison pval <0.01.</figcaption></figure>

doi:10.1371/journal.pone.0141482.g005

biosynthesis (Fig 3). At the small antral stage granulosa cells over-expressed the *stearoyl-CoA desaturase* (SCD) gene that is involved in lipid metabolism which had previously been observed as being weakly expressed in rat primordial follicles and mainly expressed in granulosa cells of antral follicles [43]. Other differential genes concerned steroid synthesis such as steroid-5-alpha-reductase, alpha polypeptide 1 (SRD5A1), hydroxy-delta-5-steroid dehydrogenase, 3 beta- and steroid delta-isomerase 1 (HSD3β) and hydroxysteroid (17-beta) dehydrogenase (HSD17β) and steroid regulation such as the nuclear receptor subfamily 5, group A, member 2 (NR5A2) [44] and luteinizing hormone/choriogonadotropin receptor (LHCGR). In line with these changes, IPA predicted an increase in the steroid biosynthesis process (S5 Table). Gene expression profiles are summarized in Fig 7.

## Putative regulators involved during early folliculogenesis

IPA Upstream Regulator analysis was used in order to understand the cause of the gene expression changes observed. The relevance and usefulness of this analysis (S6 Table) were confirmed by the identification of important molecules and pathways known to affect or be involved in follicular growth. Indeed, our analysis confirmed the important roles of NOBOX oogenesis homeobox (NOBOX), KITL and PDGF at key stages in early development.

Genes that were differentially expressed at small antral stage enabled prediction of an increased activity of the NOBOX protein in ovine oocytes and of KITL protein in ovine
<figure><img src="image_9.png" /><figcaption>**Fig 6.** Dynamic representation of integrin family gene expression. Significant differences in relative gene expression during early folliculogenesis of members of integrin family (S2 Table). Granulosa cells data are colored in red and oocyte data are colored in blue. The Y axis corresponds to normalized counts from RNA-seq data. *: FDR <0.05, pairwise comparison pval <0.01.</figcaption></figure>

doi:10.1371/journal.pone.0141482.g006

granulosa cells (Fig 8A and 8B), and a decreased activity of PDGFC in oocytes. In mice, Nobox activates expression of the *BMP15*, *DNA (cytosine-5)-methyltransferase 1 (DNMT1)*, *GDF9*, *MOS*, *zygote arrest 1 (ZAR1)* genes in oocytes [45]. In our experiment, sheep oocytes exhibited an over-expression of the *BMP15*, *DNMT1*, *GDF9*, *MOS*, *ZAR1* genes at small antral stage when compared to the primordial stage and a stable expression of NOBOX during early folliculogenesis. This analysis therefore suggests an increase in NOBOX protein activity (Fig 8A). KITL protein acts at different levels of follicular development [46]. During early folliculogenesis, this protein promotes GC mitogenesis [47] via indirect mechanisms that may involve GDF-9 [3] and BMP-15 [48]. In agreement with these studies, we identified in ovine granulosa cells four down-expressed genes and 11 over-expressed genes known to be regulated by KITL (Fig 8B). In line with these changes, our RNAseq data also showed a significant over expression of *KITL* during follicular development. Last, we showed that *PDGFC* and *PDGFRA-B* are mainly expressed in granulosa cells (S2 Table) as has already been demonstrated [9], and are less expressed in oocytes. The predicted reduction in the action of PDGFC protein on oocytes at small antrum stage was in line with a decrease in *PDGFC* granulosa cell expression and *PDGFRA* and *PDGFRB* oocyte expression during early folliculogenesis in sheep (S6 and S4 Tables).

Two molecules that have been the subject of little study in the context of follicular development were predicted to significantly affect transcription during early folliculogenesis: KLF9
<figure><img src="image_10.png" /><figcaption>**Fig 7. Dynamic representation of gene expression involved in steroid biosynthesis process.**</figcaption></figure>

Significant differences in relative gene expression during early folliculogenesis of genes involved in steroid synthesis (S2 Table). \*: FDR <0.05, pairwise comparison pval <0.01.

doi:10.1371/journal.pone.0141482.g007

and BMPER (Fig 8C and 8D; S6 Table). Studies have shown that KLF9 stimulates *cyclin D1* (Ccnd1) [49] and *collagen, type I, alpha 1* (Col1a1) [50] expression and inhibits that of *calbindin 2* (Calb2) and *growth hormone receptor* (Ghr) [51]. During our experiment, we observed the opposite regulation of these genes. As a result, IPA predicted an inhibited status for the KLF9 protein (Fig 8C). In line with this, KLF9 gene expression decreased during early follicular development. KLF9 protein is involved in a large number of functions, including reproduction. It modulates the sensitivity of endometrial cells to oestradiol [52]. In pig granulosa cells, KLFs (KLF4-9-13) may recruit inhibitory complexes containing HDAC co-repressors in order to repress the transcription of *low density lipoprotein receptor* (LDLR) and *cytochrome P450, family 11* (CYP11) [53]. The role of KLF9 in oocytes still needs to be elucidated. Lastly, our findings suggest the involvement of BMPER as an oocyte regulator during primary transition to secondary follicles (Fig 8D). BMPER has been reported to antagonize the activities of BMP2-4-3-6 [54] and displays significant affinity with another BMP antagonist, chordin (CHRD) [55]. This gene is poorly documented in the ovary. In mice, the *bmper* gene is expressed from the primary follicles but is not detected in oocytes [56]. This gene is described here for the first time as being mainly expressed in sheep oocyte and up-regulated from the secondary stage (S2 Table). The downstream oocyte targets of BMPER are *actin alpha 2* (ACTA2), *activin A receptor type II-like 1* (ACVRL1), *matrix Gla protein* (MGP) and *platelet/endothelial cell adhesion molecule 1* (PECAM1), the expression of which decreased at the secondary stage in our study (Fig 8D). Further experiments are now necessary to understand the role of BMPER in folliculogenesis. In
<figure><img src="image_11.png" /><figcaption>**Fig 8.** Upstream regulators involved during early folliculogenesis. This Figure focuses on four upstream regulators of interest predicted to be activated or inhibited at key transitions during follicular development. Genes in our dataset are highlighted based on up- (green) or down- (red) expression and increasing intensity in line with the degree of fold change. The predicted action of the central gene is indicated as activated (orange) or inhibited (blue) with the degree of confidence increasing in line with color intensity. Arrowheads at the ends of the interactions indicate activation, whereas bars indicate inhibitory effects. Unbroken arrows and dashed arrows represent direct and indirect interactions between genes and upstream regulators, respectively. Predicted activation of (A) NOBOX in oocytes at small antrum stage, (B) KITL in granulosa cells at small antrum stage, (C) KLF9 in oocytes at small antrum stage, (D) BMPER in oocytes at secondary stage.</figcaption></figure>

doi:10.1371/journal.pone.0141482.g008

addition, we identified the differential expression of other signaling BMP antagonists in granulosa cells or oocytes such as *follistatin* (*FST*), *WAP*, *follistatin/kazal*, *immunoglobulin, kunitz and netrin domain containing 2* (*WFIKKN2*) and *gremlin 2* (*GREM2*) (*S3 Fig*), thus highlighting the complex regulation of the BMP signaling pathway during early folliculogenesis.

## Biomarkers of early folliculogenesis

**Biomarker selection.** A number of protocols have been developed for domestic species such as cattle, buffalo and sheep regarding the performance of *in vitro* functional analysis, manipulation and embryo production, or the evaluation of cryopreservation prior to cancer treatments. The identification of biomarkers is therefore an important step because they can be used to monitor follicular development and survival. For this purpose, we performed a pair-wise comparison between each stage dataset and the other datasets including multi-tissue samples (FDR<5%, FC>3–10) and selected the 29 most differentially expressed genes (see *Materials and Methods*). The expression of this gene selection was then checked by qRT-PCR. Statistical analysis confirmed the differential expression observed in the RNA-seq data for 24 out of the 29 genes tested (*S7 Table*). Lastly, the two sets of 24 expression profiles arising from RNAseq and RT-PCR were subjected to a PLS-DA analysis to evaluate its capacity to discriminate follicular stages (*S4 Fig*). SPLS-DA analysis is a partial least squares regression technique applied to categorical variables (here the stages of follicles). PLS-DA has often been used for classification and discrimination or variable selection. As expected, these two analyses (RNA-seq/qRT-PCR) gave similar discrimination of the follicle stages, making this biomarker set easily usable by qRT-PCR
technology. The two analyses showed some structure of variables without clearly separating the stages. Consequently we turned to using parametric models for the difference in expression levels across stages, and specifically the ability of our biomarkers to predict the follicular stage. A summary of the 24 significantly differential profiles is presented in Fig 9.

**Predictive ability of biomarkers.** Based on these RT-PCR expression data, we estimated linear mixed models of expression that could be explained by developmental stage and cell type in the 24 differentially expressed genes. We then used these estimated models to derive predictive equations of stage conditional on the expression data for each cell type (oocyte and granulosa cells). Details of these procedures are given in the Methods section.

The expressions of these genes were modeled jointly using linear mixed models. These data displayed both a binary nature (differences in the presence/absence of expression across stages and cell types) and a quantitative nature (different levels of expression). This behavior was clarified using two different models: 1- a model for the presence/absence of expression using a hierarchical logistic regression model, and 2- hierarchical linear regression to determine the quantitative level of expression.

Both models included parameters to fit gene and cell type effects within stages relative to either the probability (model 1) or quantitative level (model 2) of expression. With model 1, we found that across the genes, the probability of expression at the PD, PM and SC but not SA stages were highly correlated (>0.99). Within cell types, the correlations between the PD and PM stages on the one hand, and the SC and SA stages on the other, were also very high (>0.9). Using model 2, we found that the levels of expression across the genes between the PD and PM stages were highly correlated (>0.99) but that the expression levels between the PD or PM and SA stage were highly negatively correlated (<-0.99).

Using the estimated models, we then derived predictive equations for the follicular stage from a vector of observed expression levels on all genes. We then simulated 100 new observation vectors from the original data, for each stage (see Methods for details; S1 Text). Fig 10 presents the posterior probability of an expression vector arising from each of the possible stages, when the simulated vector was made up of observations from the PD, PM, SC or SA stages (from left to right), for granulosa (top) and oocyte (bottom) cells. In Fig 10, both the presence and absence of expression, and quantitative expression levels, are taken into account. The predictive ability of the model is very good, with the posterior probability of the correct stage being almost always equal to one for all stages. We noted however that this predictive ability was higher for oocytes than the granulosa cells. We also found that the binary response alone did not ensure good predictive ability for most stages (see S5 Fig). For example, when sampling a vector from the PD stage in oocytes, the posterior probabilities were spread almost evenly between the PD and PM stages.

Our results therefore suggest that our predictive equations could be used to obtain reliable data on the follicular stage of a cell from the expression data for this set of genes. However, we noted that this predictive ability might be over-estimated as the random vectors generated were obtained by random sampling of the data used to estimate the parameters for the predictive equations. Determining this predictive ability on independent data was not possible here, but this issue will be addressed in the context of a future study.

**Biological relevance of biomarkers.** The list of biomarkers shown in Fig 9 includes known genes such as *SPO11*, *BMP15*, *FST* and *WEE2* which play critical roles in meiosis, primary-secondary transition or early folliculogenesis. For example, *BMP15* and *FST* expression profiles indicated cell type and stage-specific gene expressions as previously described [11]. This confirms the level of purity of the LCM-derived aRNA samples, and the quality of our RNAseq analysis. Other biomarkers are also likely to play significant role in follicular development.
<figure><img src="image_12.png" /><figcaption>**Fig 9.** Distribution of the expression of the 24 biomarkers during early follicular development. These biomarkers correspond to the most strongly expressed genes in a specific cell type and stage (PDO, PMO, SCO, SAO, PDG, PMG, SCG, SAG). They were selected using pairwise comparisons (nbinom Test) comparing each stage to the others including multi-tissue samples (FDR<5% and a fold change >3 for granulosa samples and a fold change >-10 for oocyte samples). Their differential expression was confirmed by qRT-PCR. Gene names marked in red are genes expressed in granulosa cells and gene names marked in blue are genes expressed in oocytes. Gene names marked in black are genes expressed both in oocytes and granulosa cells.</figcaption></figure>

doi:10.1371/journal.pone.0141482.g009

<figure><img src="image_13.png" /><figcaption>**Fig 10.** Predictive power of the combination of logistic regression and the linear regression model. Predictive power of biomarkers using linear mixed model equations incorporating both the presence/absence of expression and its quantitative level. The scatter-plot shows the posterior probability that an expression vector arises from each of the possible stages, when the simulated vector is made up of observations from the PD, PM, SC or SA stages (from left to right), for granulosa cells (top) and oocytes (bottom).</figcaption></figure>

doi:10.1371/journal.pone.0141482.g010
At the primordial stage, we highlighted the oocyte expression of *potassium channel, two pore domain subfamily K, member 9 (KCNK9)*, *sparc/osteonectin, cwcv and kazal-like domains proteoglycan (testican) 1 (SPOCK1)*, *insulin-like growth factor 2 mRNA binding protein 1 (IGF2BP1)* and the granulosa cell expression of *microtubule-associated protein 1 light chain 3 gamma (MAP1LC3C)*. KCNK9 encodes a potassium channel family protein that has already been identified in the membrane of cow oocytes [57]. SPOCK1 encodes for a proteoglycan protein. It was initially identified in human testis and tissue distribution studies have shown that it is most prominently expressed in nervous tissue. It can block the substrate attachment and neurite outgrowth of neuronal cells [58] and may therefore regulate brain development. Until now, no studies have explored SPOCK1 in the ovary. IGF2BP1 belongs to a family of RNA-binding proteins implicated in mRNA localization, turnover and translational control. The spatio-temporal expression of IGF2BP1 protein has been described in the female gonads of mouse and humans [59]. Its expression is ubiquitous up to 16.5 days in the mouse embryonic stage and then becomes restricted to germ cells. There is usually no expression in adult tissues, with the exception of gonads. In adult mouse and human ovaries, IGF2BP1 protein has been detected in both resting and growing oocytes and less clearly detected in granulosa cells. IGF2BP1 mRNA is nevertheless a good marker of primordial oocytes in sheep. The MAP1LC3C gene is mainly expressed in granulosa cells of primordial follicles. This gene is involved in the formation of autophagosomal vacuoles, an autophagy activity that is fundamental to many processes across the reproduction spectrum, from development of the primordial follicle and spermatozoa to embryogenesis [60]. For example, autophagy is involved in maintaining the primordial follicles in a dormant state via non apoptotic mechanism [61, 62]. Our results reinforce these findings. Later, at antral stages, it was also shown that autophagy process is also implicated in follicular atresia in different species [63, 64] and may be regulated through cross talk with apoptosis proteins [65].

At the primary stage we identified expression of the *tectorin alpha (TECTA)* gene in sheep oocytes. TECTA protein is one of the major noncollagenous components of the tectorial membrane and plays a role in intracochlear sound transmission. This gene was not documented in other tissues.

At the secondary stage we highlighted the expression of the *cannabinoid receptor 1 (CNR1)* gene in granulosa cells. This gene encodes one of cannabinoid receptors known to regulate spermatic differentiation. Its invalidation influences sperm chromatin remodelling by reducing the displacement of histones [66]. Its role in granulosa cells still needs to be defined.

Finally at the small antral stage we highlighted the granulosa cell expression of the *synapse differentiation inducing 1 (SYNDIG1)* gene and the oocyte expression of the *SPO11*, *WEE2* and *B-cell translocation gene 4 (BTG4)* genes. SYNDIG1 encodes the SYNG1 protein that belongs to the interferon-induced transmembrane family of proteins. This protein is a regulator of excitatory synapse maturation and induces synapse maturation [67]. Its role during follicular development requires further investigation. Surprisingly, during this study, we identified a preferential expression of SPO11 in small antrum oocytes. It is assumed that SPO11 is necessary to generate chromatin breaks during leptotene of the first meiotic prophase in yeast and mammals [58, 68]. In early ovarian follicles, this process is completed and germ cell meiosis is arrested at the diplotene stage. The function of SPO11 after DNA double-strand breaks (DSBs) has been little studied and it would be interesting to investigate its role as from the end of prophase I. Thus Houmard et al showed in the human fetal ovary that the number of SPO11 transcripts was elevated after 12 weeks of gestation (corresponding to the end of meiosis) and remained elevated up to 18 weeks (initiation of follicle formation at 12–14 weeks) [69]. The detection of SPO11 expression in sheep oocytes of small antrum follicles was in line with this previous human study so we can hypothesize new roles for SPO11 during oocyte development. Wee2 is responsible for phosphorylating the CDK1 inhibitory site and maintaining meiotic arrest in oocytes. In macaque monkeys species, and as seen during our study, WEE2 mRNA appears to accumulate during folliculogenesis,
reaching the highest level in preovulatory follicles [29]. BTG4 belongs to the PC3/BTG/TOB family of cell-cycle inhibitors. It is able to induce G₁ arrest, and is highly expressed in the testis, oocytes, and preimplantation embryos [70]. BTG4 is also highly expressed in bovine testis and ovary/oocyte, which suggests a role for BTG4 in meiosis [71, 72].

This biomarker set thus offers an “*in vivo*” reference to predict follicular development by a
useable qRT-PCR technology.

**Conclusions**

For the first time, a detailed and comprehensive spatio-temporal transcriptome for the develop-
mental stages of early folliculogenesis has been presented in this study. Substantial changes in
gene expression occur in both oocytes and granulosa cells as the follicles develop. We have
described the specific gene expression dynamic during follicular development for each compart-
ment. We identified more gene expression change in oocytes than in granulosa cells, and revealed
important differential expressions between primary and secondary stages in sheep. Gene expres-
sion was little affected in granulosa cells at the last studied stage (small antrum stage).

Our analysis highlighted the *in vivo* molecular course of a number of processes. As early as the primary stage, a significant number of genes involved in the cell cycle regulation were found to be differentially expressed in oocytes, illustrating initiation of the acquisition of meiotic competence. Other changes in gene expression affected oocyte migration and cellular organization. At the secondary stage, genes related to cell migration signaling pathways such as integrins, WNT or Rho were mainly down-regulated in oocytes. In neighboring cells (granulosa cells), an interesting subset of genes was found to be involved in the adhesion and formation of cytoplasmic projections. The expression of integrin family members, involved in cell morphology, was notably affected in both compartments. At the small antral stage, the acquisition of meiotic competence progressed with an over-expression of genes required to stabilize mouse oocyte meiotic arrest. We demonstrated an increase in the expression of genes already known to participate in lipid and steroid biosynthesis in granulosa cells. Lastly, some important pathways were found to be associated with early follicular development, such as FGF signaling in oocytes and WNT and BMP signaling in both compartments.

Numerous questions remain regarding the mechanisms governing the regulation of early follicular growth. Our findings proposed the new involvement of a number of molecules such as KLF9 and BMPER in gene regulation during this growth. We also generated a data set of 24 biomarkers that could be useful to predict follicular growth. They also are likely to play important roles in the process of follicular growth. They enabled discrimination between early follicular stages and thus offer a molecular signature for early follicular growth. They are of importance to determining the quality and correct development of *in vitro* or *in vivo* protocols for the study of early folliculogenesis, the production of embryos or estimation of the cryogenic preservation of the ovarian cortex.

This study has therefore highlighted new arguments that enable a better understanding of
the molecular mechanisms involved in early ovarian follicular development and new insights
into both reproductive biology and regenerative medicine.

**Supporting Information**

S1 Fig. Analysis of changes in expression during stage transition. This Figure shows the number of differentially expressed genes (under- and over-expressed) in oocyte (A) and granulosa cells (B) at each follicular transition: primordial/primary, primary/secondary, secondary/small antrum.

(TIF)
**S2 Fig.** Functional enrichment during early follicular development. Genes differentially expressed during early development were evaluated *in silico* using Ingenuity Pathway Analysis (IPA) for each compartment (FDR <5%). The bar color corresponds to enriched functions at primary stage (compared to primordial stage), secondary stage (compared to primary and primordial stages), and small antrum stage (compared to secondary, primary and primordial stages). The X axis corresponds to the level of significance of the function: -log(B-H p value). Granulosa cell functions are colored in red and oocyte functions are colored in blue. Numbers correspond to the numbers of focus genes that contributed to the functions.

(TIF)

**S3 Fig.** Gene expression profiles of the genes involved in BMP signaling. Y axis corresponds to normalized counts from RNA-seq data. GC data are colored in red and oocyte data are colored in blue.

(TIF)

**S4 Fig.** (s)PLS-DA analysis of the biomarker set transcriptome. (s) PLS-DA was performed on the biomarker dataset after normalization using the DEseq R package to classify follicular stages according to gene expression. The Figure visualizes the first three components of the analysis from: (A) RNA-seq dataset, (B) qRT-PCR dataset.

(TIF)

**S5 Fig.** Predictive power of the logistic regression model. Predictive power of biomarkers using linear mixed model equations incorporating only the presence/absence of expression. The scatter-plot shows the posterior probability that an expression vector arises from each of the possible stages, when the simulated vector is made up of observations from the PD, PM, SC or SA stages (from left to right), for granulosa cells (top) and oocytes (bottom).

(TIF)

**S1 Table.** Primer sequences for real-time PCR.

(XLSX)

**S2 Table.** Genes significantly differentially expressed during early follicular development.

(XLSX)

**S3 Table.** QRT-PCR validation. The expression profiles of 19 genes of interest involved in enriched canonical pathways were monitored using qRT-PCR, and statistical analysis confirmed the DE observed in the RNA-seq dataset for 14 of them.

(DOCX)

**S4 Table.** Significantly enriched oocyte and granulosa cell canonical pathways during early follicular development. Significant pathway enrichment in differentially expressed genes was investigated using webgestalt software (FRD<0.05).

(XLSX)

**S5 Table.** Analysis of downstream effects. This analysis predicted the effect of change in gene expression on the functions using IPA software.

(XLSX)

**S6 Table.** Analysis of upstream effects. IPA Upstream Regulator analysis was used to obtain clues regarding the cause of the change in gene expression and to provide more evidence regarding the biological mechanism.

(XLSX)
S7 Table. qRT-PCR validation of gene expression in the biomarker set.

(XLSX)

S1 Text. Description of Linear Mixed models.

(DOCX)

Acknowledgments

We would like to thank Katia Feve for her help with the Fluidigm experiment. Special thanks
to V. Hawken for English corrections.

Author Contributions

Conceived and designed the experiments: AB PM BMP. Performed the experiments: AB. Analyzed the data: AB BS. Contributed reagents/materials/analysis tools: AB BS. Wrote the paper: AB BS PM BMP.

References

1. Choi Y, Rajkovic A. Genetics of early mammalian folliculogenesis. Cell Mol Life Sci. 2006; 63(5):579-90. PMID: 16416028.

2. Fair T, Hyttel P, Greve T. Bovine oocyte diameter in relation to maturational competence and transcriptional activity. Mol Reprod Dev. 1995; 42(4):437-42. doi: 10.1002/mrd.1080420410 PMID: 8607973.

3. Gilchrist RB, Ritter LJ, Armstrong DT. Oocyte-somatic cell interactions during follicle development in mammals. Anim Reprod Sci. 2004;82-83:431-46. PMID: 15271471.

4. Matzuk MM, Burns KH, Viveiros MM, Eppig JJ. Intercellular communication in the mammalian ovary: oocytes carry the conversation. Science. 2002; 296(5576):2178-80. doi: 10.1126/science.1071965 PMID: 12077402.

5. Kezele PR, Ague JM, Nilsson E, Skinner MK. Alterations in the ovarian transcriptome during primordial follicle assembly and development. Biol Reprod. 2005; 72(1):241-55. PMID: 15371273.

6. Dharma SJ, Modi DN, Nandedkar TD. Gene expression profiling during early folliculogenesis in the mouse ovary. Fertil Steril. 2009; 91(5 Suppl):2025-36. PMID: 18504043. doi: 10.1016/j.fertnstert.2008.02.088

7. Markholt S, Grondahl ML, Ernst EH, Andersen CY, Ernst E, Lykke-Hartmann K. Global gene analysis of oocytes from early stages in human folliculogenesis shows high expression of novel genes in reproduction. Molecular Human Reproduction. 2012; 18(2):96-110. doi: 10.1093/molehr/gar083 WOS:000299786200006. PMID: 22238370

8. Pan H, O'Brien MJ, Wigglesworth K, Eppig JJ, Schultz RM. Transcript profiling during mouse oocyte development and the effect of gonadotropin priming and development in vitro. Dev Biol. 2005; 286 (2):493-506. PMID: 16168984.

9. Yoon SJ, Kim KH, Chung HM, Choi DH, Lee WS, Cha KY, et al. Gene expression profiling of early follic-ular development in primordial, primary, and secondary follicles. Fertil Steril. 2006; 85(1):193-203. PMID: 16412753.

10. Magalhães-Padilha DM, Geisler-Lee J, Wischral A, Gastal MO, Fonseca GR, Eloy YR, et al. Gene Expression During Early Folliculogenesis in Goats Using Microarray Analysis. Biol Reprod. 2013. doi: 10.1095/biolreprod.112.106096 PMID: 23759311.

11. Bonnet A, Bevilacqua C, Benne F, Bodin L, Cotinot C, Liaubet L, et al. Transcriptome profiling of sheep granulosa cells and oocytes during early follicular development obtained by laser capture microdissection. BMC Genomics. 2011; 12:417. doi: 10.1186/1471-2164-12-417 PMID: 21851638; PubMed Central PMCID: PMCPMC3166951.

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

13. Archibald AL, Cockett NE, Dalrymple BP, Faraut T, Kijas JW, Maddox JF, et al. The sheep genome reference sequence: a work in progress. Anim Genet. 2010; 41(5):449-53. doi: 10.1111/j.1365-2052.2010.02100.x PMID: 20809919.

14. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biology. 2010; 11(10). doi: 10.1186/gb-2010-11-10-r106 WOS:000287378900008.
15. Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001; 29(9):e45. PMID: 11328886.

16. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002; 3(7):RESEARCH0034.

17. Tieri P, Nardini C. Signalling pathway database usability: lessons learned. Mol Biosyst. 2013; 9 (10):2401-7. doi: 10.1039/c3mb70242a PMID: 23942525.

18. Liu K, Rajareddy S, Liu L, Jagarlamudi K, Boman K, Selstam G, et al. Control of mammalian oocyte growth and early follicular development by the oocyte PI3 kinase pathway: new roles for an old timer. Dev Biol. 2006; 299(1):1-11. doi: 10.1016/j.ydbio.2006.07.038 PMID: 16970938.

19. Cahill LP, Mauléon P. Influences of season, cycle and breed on follicular growth rates in sheep. J Reprod Fertil. 1980; 58(2):321-8. PMID: 7191893.

20. Gougeon A. Regulation of ovarian follicular development in primates: facts and hypotheses. Endocr Rev. 1996; 17(2):121-55. PMID: 8706629.

21. Hamatani T, Yamada M, Akutsu H, Kuji N, Mochimaru Y, Takano M, et al. What can we learn from gene expression profiling of mouse oocytes? Reproduction. 2008; 135(5):581-92. doi: 10.1530/REP-07-0430 PMID: 18411407.

22. Fnu S, Williamson EA, De Haro LP, Brenneman M. Wray J, Shaheen M, et al. Methylation of histone H3 lysine 36 enhances DNA repair by nonhomologous end-joining. Proc Natl Acad Sci U S A. 2011; 108 (2):540-5. doi: 10.1073/pnas.1013571108 PMID: 21187428; PubMed Central PMCID: PMCPMC3021059.

23. Beck BD, Lee SS, Williamson E, Hromas RA, Lee SH. Biochemical characterization of metnase's endonuclease activity and its role in NHEJ repair. Biochemistry. 2011; 50(20):4360-70. doi: 10.1021/bi200333k PMID: 21491884; PubMed Central PMCID: PMCPMC3388547.

24. Hached K, Xie SZ, Buffin E, Cladière D, Rachez C, Sacras M, et al. Mps1 at kinetochores is essential for female mouse meiosis I. Development. 2011; 138(11):2261-71. doi: 10.1242/dev.061317 PMID: 21558374.

25. Stumpff J, Wagenbach M, Franck A, Asbury CL, Wordeman L. Kif18A and chromokinesins confine centromere movements via microtubule growth suppression and spatial control of kinetochore tension. Dev Cell. 2012; 22(5):1017-29. doi: 10.1016/j.devcel.2012.02.013 PMID: 22595673; PubMed Central PMCID: PMCPMC3356572.

26. Endow SA, Kull FJ, Liu H. Kinesins at a glance. J Cell Sci. 2010; 123(Pt 20):3420-4. doi: 10.1242/jcs.064113 PMID: 20930137; PubMed Central PMCID: PMCPMC2951464.

27. Tsou AP, Yang CW, Huang CY, Yu RC, Lee YC, Chang CW, et al. Identification of a novel cell cycle regulated gene, HURP, overexpressed in human hepatocellular carcinoma. Oncogene. 2003; 22 (2):298-307. doi: 10.1038/sj.onc.1206129 PMID: 12527899.

28. Ye F, Tan L, Yang Q, Xia Y, Deng LW. Murata-Hori M, et al. HURP regulates chromosome congression by modulating kinesin Kif18A function. Curr Biol. 2011; 21(18):1584-91. doi: 10.1016/j.cub.2011.08.024 PMID: 21924616.

29. Hanna CB, Yao S, Patta MC, Jensen JT, Wu X. WEE2 is an oocyte-specific meiosis inhibitor in rhesus macaque monkeys. Biol Reprod. 2010; 82(6):1190-7. doi: 10.1095/biolreprod.109.081984 PMID: 20200212; PubMed Central PMCID: PMCPMC2874501.

30. Cánepa ET, Scassa ME, Ceruti JM, Marazita MC, Carcagno AL, Sirkin PF, et al. INK4 proteins, a family of mammalian CDK inhibitors with novel biological functions. IUBMB Life. 2007; 59(7):419-26. doi: 10.1080/15216540701488358 PMID: 17654117.

31. Wei Z, Jiang X, Liu F, Qiao H, Zhou B, Zhai B, et al. Downregulation of Skp2 inhibits the growth and metastasis of gastric cancer cells in vitro and in vivo. Tumour Biol. 2013; 34(1):181-92. doi: 10.1007/s13277-012-0527-8 PMID: 23229098.

32. Li Y, Pal R, Sung LY, Feng H, Miao W, Cheng SY, et al. An opposite effect of the CDK inhibitor, p18 (INK4c) on embryonic stem cells compared with tumor and adult stem cells. PLoS One. 2012; 7(9): e45212. doi: 10.1371/journal.pone.0045212 PMID: 23049777; PubMed Central PMCID: PMCPMC3458833.

33. Oh JS, Susor A, Schindler K, Schultz RM, Conti M. Cdc25A activity is required for the metaphase II arrest in mouse oocytes. J Cell Sci. 2013; 126(Pt 5):1081-5. doi: 10.1242/jcs.115592 PMID: 23345398; PubMed Central PMCID: PMCPMC3635459.

34. Beltowski J, Jamroz A. Adrenomedullin—what do we know 10 years since its discovery? Pol J Pharma-col. 2004; 56(1):5-27. PMID: 15047974.
35. Dong YL, Reddy DM, Green KE, Chauhan MS, Wang HQ, Nagamani M, et al. Calcitonin gene-related peptide (CALCA) is a proangiogenic growth factor in the human placental development. Biol Reprod. 2007; 76(5):892-9. doi: 10.1095/biolreprod.106.059089 PMID: 17267696.

36. Katoh-Fukui Y, Tsuchiya R, Shiroishi T, Nakahara Y, Hashimoto N, Noguchi K, et al. Male-to-female sex reversal in M33 mutant mice. Nature. 1998; 393(6686):688-92. doi: 10.1038/31482 PMID: 9641679.

37. Boutin C, Goffinet AM, Tissir F. Celsr1-3 cadherins in PCP and brain development. Curr Top Dev Biol. 2012; 101:161-83. doi: 10.1016/B978-0-12-394592-1.00010-7 PMID: 23140629.

38. Kozyulina PY, Loskutov YV, Kozyreva VK, Rajulapati A, Ice RJ, Jones BC, et al. Pro-metastatic NEDD9 regulates individual cell migration via caveolin-1-dependent trafficking of integrins. Mol Cancer Res. 2014. doi: 10.1158/1541-7786.MCR-14-0353 PMID: 25319010.

39. Fujiwara H, Maeda M, Honda T, Yamada S, Ueda M, Kanzaki H, et al. Granulosa cells express integrin alpha 6: possible involvement of integrin alpha 6 in folliculogenesis. Horm Res. 1996; 46 Suppl 1:24-30. PMID: 8864745.

40. Giebel J, Rune GM. Relationship between expression of integrins and granulosa cell apoptosis in ovarian follicles of the marmoset (Callithrix jacchus). Tissue Cell. 1997; 29(5):525-31. PMID: 9364802.

41. Shuttleworth G, Broughton Pipkin F, Hunter MG. In vitro development of pig preantral follicles cultured in a serum-free medium and the effect of angiotensin II. Reproduction. 2002; 123(6):807-18. PMID: 12052235.

42. Tasaki H, Iwata H, Sato D, Monji Y, Kuwayama T. Estradiol has a major role in antrum formation of porcine preantral follicles cultured in vitro. Theriogenology. 2013; 79(5):809-14. doi: 10.1016/j.theriogenology.2012.12.009 PMID: 23347746.

43. Moreau C, Froment P, Tosca L, Moreau V, Dupont J. Expression and regulation of the SCD2 desaturase in the rat ovary. Biol Reprod. 2006; 74(1):75-87. doi: 10.1095/biolreprod.105.044545 PMID: 16207839.

44. Parker KL, Rice DA, Lala DS, Ikeda Y, Luo XR, Wong M, et al. Steroidogenic factor 1: An essential mediator of endocrine development. Recent Progress in Hormone Research, Vol 57: Reproductive Hormones & Human Health. 2002; 57:19-36. WOS:000179991300002.

45. Rajkovic A, Pangas SA, Ballow D, Suzumori N, Matzuk MM. NOBOX deficiency disrupts early folliculogenesis and oocyte-specific gene expression. Science. 2004; 305(5687):1157-9. PMID: 15326356.

46. Klinger FG, De Felici M. In vitro development of growing oocytes from fetal mouse oocytes: stage-specific regulation by stem cell factor and granulosa cells. Dev Biol. 2002; 244(1):85-95. doi: 10.1006/dbio.2002.0592 PMID: 11900461.

47. Reynaud K, Cortvrindt R, Smitz J, Driancourt MA. Effects of Kit Ligand and anti-Kit antibody on growth of cultured mouse preantral follicles. Mol Reprod Dev. 2000; 56(4):483-94. doi: 10.1002/1098-2795 (200008)56:4<483::AID-MRD6>3.0.CO;2-O PMID: 10911398.

48. Otsuka F, Yao Z, Lee T, Yamamoto S, Erickson GF, Shimasaki S. Bone morphogenetic protein-15. Identification of target cells and biological functions. J Biol Chem. 2000; 275(50):39523-8. PMID: 10998422.

49. Zhang XL, Simmen FA, Michel FJ, Simmen RC. Increased expression of the Zn-finger transcription factor BTEB1 in human endometrial cells is correlated with distinct cell phenotype, gene expression patterns, and proliferative responsiveness to serum and TGF-beta1. Mol Cell Endocrinol. 2001; 181(1-2):81-96. PMID: 11476943.

50. Simmen FA, Su Y, Xiao R, Zeng Z, Simmen RC. The Krüppel-like factor 9 (KLF9) network in HEC-1-A endometrial carcinoma cells suggests the carcinogenic potential of dys-regulated KLF9 expression. Reprod Biol Endocrinol. 2008; 6:41. doi: 10.1186/1477-7827-6-41 PMID: 18783612; PubMed Central PMCID: PMCPMC2542371.

51. Gowri PM, Yu JH, Shaufl A, Sperling MA, Menon RK. Recruitment of a repressosome complex at the growth hormone receptor promoter and its potential role in diabetic nephropathy. Mol Cell Biol. 2003; 23(3):815-25. PMID: 12529387; PubMed Central PMCID: PMCPMC140700.

52. Simmen RC, Pabona JM, Velarde MC, Simmons C, Rahal O, Simmen FA. The emerging role of Krüppel-like factors in endocrine-responsive cancers of female reproductive tissues. J Endocrinol. 2010; 204(3):223-31. doi: 10.1677/JOE-09-0329 PMID: 19833720; PubMed Central PMCID: PMCPMC2971688.

53. Natesampillai S, Kerkvliet J, Leung PC, Veldhuis JD. Regulation of Kruppel-like factor 4, 9, and 13 genes and the steroidogenic genes LDLR, StAR, and CYP11A in ovarian granulosa cells. Am J Physiol Endocrinol Metab. 2008; 294(2):E385-91. doi: 10.1152/ajpendo.00480.2007 PubMed Central PMCID: PMCPMC2747322. PMID: 18056793
54. Moser M, Binder O, Wu Y, Aitsebaomo J, Ren R, Bode C, et al. BMPER, a novel endothelial cell precursor-derived protein, antagonizes bone morphogenetic protein signaling and endothelial cell differentiation. Mol Cell Biol. 2003; 23(16):5664-79. PMID: 12897139; PubMed Central PMCID: PMCPMC166349.

55. Zakin L, Chang EY, Plouhinec JL, De Robertis EM. Crossveinless-2 is required for the relocalization of Chordin protein within the vertebral field in mouse embryos. Dev Biol. 2010; 347(1):204-15. doi: 10.1016/j.ydbio.2010.08.025 PMID: 20807528; PubMed Central PMCID: PMCPMC2975673.

56. Fenwick MA, Mansour YT, Franks S, Hardy K. Identification and regulation of bone morphogenetic protein antagonists associated with preantral follicle development in the ovary. Endocrinology. 2011; 152 (9):3515-26. doi: 10.1210/en.2011-0229 PMID: 21791559.

57. Hur CG, Choe C, Kim GT, Cho SK, Park JY, Hong SG, et al. Expression and localization of two-pore domain K(+) channels in bovine germ cells. Reproduction. 2009; 137(2):237-44. doi: 10.1530/REP-08-0035 PMID: 18987255.

58. Baudat F, Manova K, Yuen JP, Jasin M, Keeney S. Chromosome synapsis defects and sexually dimorphic meiotic progression in mice lacking Spo11. Mol Cell. 2000; 6(5):989-98. PMID: 11106739.

59. Hammer NA, Hansen T, Byskov AG, Rajpert-De Meyts E, Grøndahl ML, Bredkjaer HE, et al. Expression of IGF-II mRNA-binding proteins (IMPs) in gonads and testicular cancer. Reproduction. 2005; 130 (2):203-12. doi: 10.1530/rep.1.00664 PMID: 16049158.

60. Kanninen TT, de Andrade Ramos BR, Witkin SS. The role of autophagy in reproduction from gametogenesis to parturition. Eur J Obstet Gynecol Reprod Biol. 2013; 171(1):3-8. doi: 10.1016/j.ejogrb.2013.07.020 PMID: 23932305.

61. Escobar ML, Echeverría OM, Ortíz R, Vázquez-Nin GH. Combined apoptosis and autophagy, the process that eliminates the oocytes of atretic follicles in immature rats. Apoptosis. 2008; 13(10):1253-66. doi: 10.1007/s10495-008-0248-z PMID: 18690537.

62. Reddy P, Liu L, Adhikari D, Jagarlamudi K, Rajareddy S, Shen Y, et al. Oocyte-specific deletion of Pten causes premature activation of the primordial follicle pool. Science. 2008; 319(5863):611-3. PMID: 18239123. doi: 10.1126/science.1152257

63. Rosales-Torres AM, Avalos-Rodriguez A, Vergara-Onofre M, Hernandez-Perez O, Ballesteros LM, Garcia-Macedo R, et al. Multiparametric study of atresia in ewe antral follicles: histology, flow cytometry, internucleosomal DNA fragmentation, and lysosomal enzyme activities in granulosa cells and follicular fluid. Mol Reprod Dev. 2000; 55(3):270-81. PMID: 10657046.

64. Velentzas AD, Nezis IP, Stravopodis DJ, Papassideri IS, Margaritis LH. Apoptosis and autophagy function cooperatively for the efficacious execution of programmed nurse cell death during Drosophila virilis oogenesis. Autophagy. 2007; 3(2):130-2. PMID: 17183224.

65. Shimizu S, Kanaseki T, Mizushima N, Mizuta T, Arakawa-Kobayashi S, Thompson CB, et al. Role of Bcl-2 family proteins in a non-apoptotic programmed cell death dependent on autophagy genes. Nat Cell Biol. 2004; 6(12):1221-8. doi: 10.1038/ncb1192 PMID: 15558033.

66. Cacciola G, Chioccarelli T, Altucci L, Ledent C, Mason JI, Fasano S, et al. Low 17beta-estradiol levels in CNR1 knock-out mice affect spermatid chromatin remodeling by interfering with chromatin reorganization. Biol Reprod. 2013; 88(6):152. doi: 10.1095/biolreprod.112.105726 PMID: 23677985.

67. Díaz E. SynDIG1 regulation of excitatory synapse maturation. J Physiol. 2012; 590(Pt 1):33-8. doi: 10.1113/jphysiol.2011.213884 PMID: 21878521; PubMed Central PMCID: PMCPMC3300043.

68. Celerin M, Merino ST, Stone JE, Menzie AM, Zolan ME. Multiple roles of Spo11 in meiotic chromosome behavior. EMBO J. 2000; 19(11):2739-50. doi: 10.1093/emboj/19.11.2739 PMID: 10835371; PubMed Central PMCID: PMCPMC212740.

69. Houmard B, Small C, Yang L, Naluai-Cecchini T, Cheng E, Hassold T, et al. Global gene expression in the human fetal testis and ovary. Biol Reprod. 2009; 81(2):438-43. doi: 10.1095/biolreprod.108.075747 PMID: 19369649; PubMed Central PMCID: PMCPMC2849816.

70. Buanne P, Corrente G, Micheli L, Palena A, Lavia P, Spadafora C, et al. Cloning of PC3B, a novel member of the PC3/BTG/TOB family of growth inhibitory genes, highly expressed in the olfactory epithelium. Genomics. 2000; 68(3):253-63. doi: 10.1006/geno.2000.6288 PMID: 10995567.

71. Vallée M, Gravel C, Palin MF, Reghenas H, Stothard P, Wishart DS, et al. Identification of novel and known oocyte-specific genes using complementary DNA subtraction and microarray analysis in three different species. Biol Reprod. 2005; 73(1):63-71. doi: 10.1095/biolreprod.104.037069 PMID: 15744023.

72. Pennetier S, Uzbekova S, Guyader-Joly C, Humblot P, Mermillod P, Dalbiès-Tran R. Genes preferentially expressed in bovine oocytes revealed by subtractive and suppressive hybridization. Biol Reprod. 2005; 73(4):713-20. doi: 10.1095/biolreprod.105.041574 PMID: 15930321.