Characterization of the miRNA regulators of the human ovulatory cascade

Ovarian follicular development and ovulation are complex and tightly regulated processes that involve regulation by microRNAs (miRNAs). We previously identified differentially expressed mRNAs between human cumulus granulosa cells (CGCs) from immature early antral follicles (germinal vesicle - GV) and mature preovulatory follicles (metaphase II - M2). In this study, we performed an integrated analysis of the transcriptome and miRNome in CGCs obtained from the GV cumulus-oocyte complex (COC) obtained from IVM and M2 COC obtained from IVF. A total of 43 differentially expressed miRNAs were identified. Using Ingenuity IPA analysis, we identified 7288 potential miRNA-regulated target genes. Two hundred thirty-four of these target genes were also found in our previously generated ovulatory gene library while exhibiting anti-correlated expression to the identified miRNAs. IPA pathway analysis suggested that miR-21 and FOXM1 cooperatively inhibit CDC25A, TOP2A and PRC1. We identified a mechanism for the temporary inhibition of VEGF during ovulation by TGFB1, miR-16-5p and miR-34a-5p. The linkage bioinformatics analysis between the libraries of the coding genes from our preliminary study with the newly generated library of regulatory miRNAs provides us a comprehensive, integrated overview of the miRNA-mRNA co-regulatory networks that may play a key role in controlling post-transcriptomic regulation of the ovulatory process.

In humans, miRNAs have already been profiled in cumulus versus mural granulosa cells 8 , Polycystic Ovary Syndrome (PCOS) patients 16 and cumulus cells according to oocyte nuclear maturity 17 . miRNAs have been suggested as potential biomarkers in IVF 18 .
In the past several years, miRNA and mRNA integrated analysis (MMIA) has become a tool for examining the biological functions of miRNA expression. As the biological effects of miRNAs are due to their modulation of target RNA expression, accurate miRNA target prediction is essential to any study of miRNA function.
To understand the molecular and regulatory mechanisms necessary for follicular maturation and ovulation, we aimed to uncover the miRNA transcriptome involved in these processes. Recent advances in genomics allow a systematic approach to identify critical genes and regulated pathways involved in oocyte maturation and ovulation. Using whole transcriptome sequencing, we recently identified mRNAs that are differentially expressed between immature mid-antral follicles and mature pre-ovulatory follicles 3 . The resulting database provides unprecedented insight into the processes and pathways involved in follicular maturation and ovulation. To complete the identification of factors involved in the ovulatory process, the aim of this work was to generate a library of global miRNAs involved in this process and to link the new ovulatory miRNA library with the previously described mRNA library. This enabled us to identify new regulatory mechanisms responsible for the final follicular maturation and ovulatory processes.

miRNA profile differences in compact and expanded cumulus cells.
To elucidate the roles of miR-NAs in human folliculogenesis and ovulation, we generated a library of regulated miRNAs during the final stage of follicular maturation and ovulation. We used NanoString on RNA extracted from two groups of cumulus granulosa cells: (1) Compact cumulus cells surrounding GV oocytes that were acquired during IVM treatment (CCGV n = 4) and (2) Expanded cumulus cells enclosing M2 oocytes that were acquired during IVF treatment (CCM2, n = 3). One CCGV sample was an outlier and was removed from further downstream analysis. Thus, the final number of samples in each group was three.
A Venn diagram of these groups is shown in Fig. 1B, illustrating the numbers of selective and co-expressed miRNAs in the two conditions, CCGV and CCM2. A total of 45 miRNAs were detected in CCM2 and 11 miR-NAs were detected in CCGV. Most notably, 35 miRNAs were exclusively detected in the CCM2 group; only 1 was exclusively detected in the CCGV group and 10 miRNAs were shared between the two conditions (Table 1). hsa-miR-378e was exclusively expressed in the CCGV group, although the differential expression level was not statistically significant.
With NanoString data inspected, all quality control parameters were within the acceptable ranges. Unsupervised hierarchical cluster analysis based on the expression of all detected miRNAs shows that the biological replicates cluster to their corresponding groups (Fig. 1A).
Almost all of the expressed miRNAs were significantly differentially expressed (43 of the 46 identified miR-NAs) between CCGV and CCM2, with a fold change >2 and a false discovery rate (FDR) < 5%. Of these, 25 represent unique seed sequences. Our results indicate that all differentially expressed miRNAs were up-regulated in CCM2 cells.
Validation of selected miRNA expression levels as estimated by qPCR. We used qPCR analysis to validate the results obtained by NanoString. We selected three miRNAs that were identified as differentially expressed by NanoString. qPCR was performed in three different experiments, as shown as fold induction in Fig. 1C, and demonstrated agreement with NanoString relative expression data.

Identification of differentially expressed miRNA target genes -in silico analysis. Because miR-
NAs act by regulating the expression of target genes, precise miRNA target prediction is important for elucidating miRNA function. The identification of miRNA target genes was performed using the microRNA target filter tool of QIAGEN's Ingenuity Pathway Analysis software, taking into account only experimentally validated and highly confident predictions of miRNA target interactions. We found 7288 putative miRNA target genes of which 727 (294 unique) were experimentally validated (Supp . Tables S3 and S4).
Identification of differentially expressed miRNAs that negatively correlated with ovulatory target genes -experimental analysis. To further dissect the putative ovulatory miRNA targets, we combined our previously generated mRNA library 3 and miRNA expression data using the microRNA target filter tool of QIAGEN's Ingenuity Pathway Analysis 19 . A Venn diagram depicting the relationships between the putative miRNA targets and experimentally differentially expressed mRNAs is depicted in Fig. 1D.
We chose to focus on negatively correlated targets. Since all of the differentially expressed CCM2 miRNAs were upregulated, all chosen paired mRNA targets were downregulated. The 696 downregulated mRNAs in the CCM2 group (compared to CCGV group) were paired with the 7288 putative miRNA targets. Of the 696 downregulated mRNAs, 234 were targets of 22 differentially expressed miRNAs (of the 25 having unique seed sequences). In total, these 234 genes formed 479 miRNA-target gene pairs with an inverse correlation of expression (Table 2).
To understand the putative roles of these miRNA target genes, two lists of putative regulated genes were generated. The first list included all of the putative miRNA targets (in silico analysis) and the second list included only the negatively correlated differentially expressed targets (experimental analysis). This contributes to both a broad and targeted view of the ontologies and pathways related to the roles of miRNAs in ovulation. Gene Ontology (GO) and pathway analyses were performed using the GeneAnalytics tool 20 .
Gene Ontology and pathway classification of the in silico analysis of miRNA target genes. For the unfiltered miRNA targets (7288 genes), we identified 78 high-score GO terms ( Fig. 2A   Each column represents a sample and each row represents a transcript. Expression level of each miRNA in a single sample is depicted according to the color scale. (B) Venn diagram of selective and co-expressed miRs in CCGV and CCM2 samples. A total of 45 miRs were expressed in CCM2 (pink) and 11 miRs were expressed in CCGV (green). Most notably, 35 miRs were exclusively expressed in the IVF group and 10 miRs were shared by both samples. (C) Total mRNA was purified from CCs denuded from GV COC aspirated during IVM procedures (CCGV) and CCs denuded from M2 COC (CCM2) aspirated during IVF procedures. The miRNAs were subjected to qPCR in duplicate with the examined genes and RNU6B primers. Gene expression was calculated relative to the RNU6B level in the same sample and expression levels were compared using Student's t-test. The difference reached a level of significance (P < 0.05) for all tested genes. NanoString (green) and qPCR (blue) results are presented as log2fold change between CCM2 and CCGV samples. (D) Venn diagram showing the relationship between the putative miRNA targets (brown) and experimentally differentially expressed mRNAs (from Yerushaslmi et al. 3 , upregulated mRNA (purple), downregulated mRNA (red)). Of all the putative miRNA targets, 234 were negatively correlated.  Table 1. miRNAs identified by the NanoString nCounter miRNA expression assay in human cumulus granulosa cells. Forty-three of the miRNAs were differentially expressed with a fold change >2 and an FDR < 5%. *miRNA detected as one of the 10 most abundant in human cumulus cells 8 . **Differential expression between CCGV and CCM2 samples was not significant.

Count of
For the subset of negatively correlated miRNA targets (234 genes), only four high scoring and four medium scoring GO terms were identified. The GO terms included DNA replication, ascending aorta morphogenesis, extracellular matrix organization and mitotic cell cycle. The top medium scoring GO terms included collagen fibril, basement membrane organization and signal transduction ( Fig. 2C and Supp. Table S7). The nine high scoring pathways included RB in cancer, integrated pancreatic cancer pathway, diseases associated with O-glycosylation, cell cycle, E2F-mediated DNA replication, degradation of extracellular matrix and ERK signaling ( Fig. 2D and Supp. Table S8).
Upstream regulator analysis of negatively correlated miRNA target genes. We used Ingenuity's IPA upstream regulator analytic feature to identify a cascade of upstream transcriptional regulators that can explain the observed gene expression changes of the negatively correlated genes beyond the identified miRNAs. These regulators can be transcription factors (TFs) and any gene or small molecule that has been observed experimentally to affect gene expression. TF and miRNA may mutually regulate each another or regulate a shared target gene and, as such, enhance the robustness of gene regulation. Hence, TF-miRNA regulatory network analysis will be helpful to decipher gene expression regulatory mechanisms of the ovulatory process. A total of 68 significant upstream regulators were identified (Z-score > ±2, Supp. Table S9). The top ten inhibited regulators and activated regulators are presented in Tables 3 and 4, respectively. As expected, among the activated regulators, we identified 5 miRNAs, 3 of which were also differentially expressed in our analysis (Supp . Table S9).
We plotted several of the upstream regulators and their targets, as shown in Fig. 3. We identified miR-21 to be involved in the regulation of CDC25A, DDAH1, ECT2, MSH2, NUSAP1, PBK, PRC1, RAD51AP1, STMN1 and TOP2A (Fig. 3A). Interestingly, most genes are involved in reproductive disease (highlighted in pink).

Discussion
A number of previous studies have examined mRNA and/or miRNA expression in cumulus granulosa cells. However, no study has examined global miRNA expression in human cumulus granulosa cells during final follicular maturation and ovulation. Moreover, this is the first study to link the miRNA-target gene pairs with an inverse correlation of expression. This approach allows the identification of the regulated mRNA-miRNA networks during final follicular maturation and ovulation. Differential expression analysis revealed dramatic changes in miRNA expression during the CC maturation process; all differentially expressed miRNAs were up-regulated in CCM2 cells.
Compared to published miRNA expression in CCs obtained from M2 COC during IVF treatment 8 , we observed that 8 of the 10 most abundant miRNAs were also present in the CCM2 sample and 43 of the 46 expressed miRNAs were identified by Velthut-Meikas et al.
Forty-two of the differentially expressed miRNAs were previously identified in granulosa cells. Among them, 14 are involved in steroidogenesis or PCOS, 14 are involved in apoptosis (eight of which belong to the let-7 family), 12 have been implicated in bovine follicular development and 20 were described as differentially expressed between human cumulus and mural granulosa cells. miR-4454 was not previously described in granulosa cells and was recently implicated in the pathogenesis of cartilage degeneration by promoting inflammatory, catabolic, and cell death activity in chondrocytes 21 . We believe that miR-4454 is implicated in the control of final follicular maturation perhaps by promoting similar proinflammatory-like activities in granulosa cells.
The differentially expressed miRNAs were predicted to regulate 696 target genes that were differentially expressed in our previously generated ovulatory cDNA library, of which, 234 displayed an inverse correlation.
During the final stages of ovulation, several upregulated processes were previously identified by the analysis of global mRNA expression, such as cellular movement, inflammatory response, immune cell trafficking, tissue development and lipid metabolism 3,22,23 . However, several processes, such as tissue morphology, DNA replication and cell cycle, are downregulated 3,22,23 . These downregulated processes were also represented in our analyses of miRNA targets. We observed that a large number of the putative miRNA target genes are related to cell cycle and DNA replication processes. These are probably downstream to the FSH blocking effect following the LH surge, thus leading to decreased proliferation 24 . Interestingly, all of the differentially expressed miRNAs were upregulated in our analysis. This result suggests that the main role of miRNA in the late ovulatory process is more toward control and inhibition of granulosa cell proliferation and less toward other processes found when analyzing gene, rather than miRNA, expression.
The top identified pathways among the miRNA negatively correlated targets include several pathways that pertain to cellular proliferation and survival, including "RB in cancer", "E2F regulation of DNA replication" and "cell cycle". Among them, the retinoblastoma protein (RB) was identified as an important factor in ovarian physiology. Conditional deletion of the RB gene in murine ovarian granulosa cells leads to increased follicular recruitment and premature ovarian failure 25 .
We also identified the O-glycosylation of proteins pathway, which might pertain to cumulus matrix maintenance since cumulus complexes from O-glycan-deficient oocytes were smaller and contained fewer CCs, although fertility was not impaired 26 .
As expected, more potential novel pathways, including NFAT (nuclear factor of activated T-cells) and PAK (p21-activated kinase) signaling pathways, were found in the broader analysis, which included all miRNA targets (and not only the anti-correlated ones).  NFAT has not been described in granulosa cells, although it has been implicated in GnRH signaling 37 . NFAT was also described as a regulator of COX2 in endometrial stromal cells 38 and other tissues and may play a role in prostaglandin regulation during ovulation.
PAK signaling was implicated in cell migration by altering actin cytoskeletal dynamics downstream to Rac/ Cdc42, although this activity was not described in granulosa cells and may be important during cumulus expansion and ovulation 39 . miR-21 is highly expressed and upregulated by LH in murine granulosa cells 40 and ovine follicles 12 . miR-21 suppression results in granulosa cell apoptosis 40 . Our results further corroborate miR-21 upregulation in humans during the ovulatory process and, by combining our negatively correlated ovulatory gene expression data, we suggest a mechanism of action for miR-21 function in granulosa cells. We found 10 different transcripts regulated by miR-21 that were downregulated in our library, all of which are involved in the cell cycle and apoptosis. These findings expand our knowledge of the processes that lead to the observed inhibition of proliferation on granulosa cells during the final stages of the ovulatory process. Some of the targets were already identified in granulosa cells, including CDC25a (apoptosis, cell cycle 41 ), MSH2 (apoptosis 42 ), PRC1 (apoptosis 43 ) and STMN1 (cell migration, apoptosis 44 ). Other targets that were not previously described in granulosa cells include TOP2A (apoptosis, cell cycle), DDAH1 (apoptosis), ECT2 (cell cycle), RAD51AP1 (cell cycle), PBK (apoptosis) and NUSAP1 (cell migration, apoptosis). The known and previously unknown targets provide insight into miR-21 activity in granulosa cells. The negative correlation between TOP2A and miR-21 was experimentally validated (Figs 1C and 4B).
Forkhead box M1 (FOXM1), a member of the large family of Forkhead box transcription factors, is highly expressed in proliferating cells and plays pivotal roles in embryonic and fetal development, DNA replication and mitosis 45 . The role of FOXM1 in ovulation and granulosa cells has not yet been established. FOXM1 is downregulated in granulosa cells obtained from obese patients undergoing IVF compared to lean patients 46 and was also found as an inhibited upstream regulator of negatively correlated miRNA target genes in our analysis. Our observation suggests that the reciprocal activation of miR-21 and suppression of FOXM1 synergize to cause the observed inhibition of granulosa cell proliferation during the final follicular maturation just prior to ovulation. Interestingly, high expression of miR-21 is implicated in increased proliferation in cancer and inflammation 47 . However, the inhibition of breast cancer cell growth by 3, 3′-diindolylmethane is mediated by the upregulation of miR-21 and CDC25A and FOXM1 suppression 48 , consistent with that observed in cumulus cells during ovulation 3 . These in silico findings between FOXM1, TOP2A, and miR-21 were experimentally validated (Figs 1C and 4B).
The let-7 family of miRNAs is believed to act as tumor suppressors 49 and has been implicated in granulosa cell atresia in porcine 50 . This family was among the most upregulated miRNAs during final follicular maturation. Several of its negatively correlated targets were previously identified in granulosa cells, including AURKB (cell cycle 51 ), CDC25A (apoptosis, cell cycle 41 ), EZH2 (histone methylation 52 ), and THBS1 (adhesion, angiogenesis 53 ). Other targets that were not described in granulosa cells include BRCA2 (DNA repair), CDC6 (DNA replication), FANCD2 (DNA repair), MCM3 (DNA replication), PPP1R12B (myosin phosphatase), RRM2 (DNA synthesis) and CHEK1 (cell cycle). Most negatively correlated let-7 targets are related to DNA replication/repair and cell cycle progression. Thus, the putative role of let-7 upregulation during final follicular maturation may be in mediating the observed inhibition of granulosa cell proliferation.  19 to examine the crosstalk between the TGFB1 signaling pathway and miR-34a and miR-16a-5p (seed of miR-424-5p). We revealed a number of genes and pathways, most notably VEGFA, that are reciprocally regulated by these negatively correlated transcription regulators. (B) Total mRNA was purified from CCs denuded from GV COC and M2 COC aspirated during IVF procedures. The mRNAs were subjected to qPCR in duplicate with the examined genes and ACTB primers. Gene expression was calculated relative to the ACTB level in the same sample and expression levels were compared using Student's t-test. The difference reached p = 0.006 for FOXM1, p = 0.01 for TOP2A and p = 0.08 for CD47. RNAseq (green) and qPCR (blue) results are presented as log2-fold change between COC M2 and COC GV samples. We determined that miR-34a is upregulated in CCM2 compared to CCGV. miR-34a is a tumor suppressor gene that suppresses ovarian cancer proliferation and motility by targeting AXL receptor tyrosine kinase 54 . miR-34a levels are increased in sheep granulosa cells during the follicular to luteal transition 12 . In addition, miR-34a may be involved in granulosa cell apoptosis in pig ovaries by targeting the inhibin B gene 55 . The integrated mRNA-miRNA analysis suggests that miR-34a acts through the regulation of E2F1 (transcription), CD47 (adhesion), MCM3 (proliferation), VEGFA (angiogenesis) and TGFB1 signaling during final follicular maturation and luteinization.
Furthermore, the roles of the TGFB1 signaling pathways in folliculogenesis have been extensively studied 56 . The Ingenuity's IPA upstream regulator analytic identified 38 inhibited TGFB1 targets of the 234 putative miRNA negatively correlated targets, making it the most significant miRNA-regulated pathway. When TGFB1 targets were plotted along with miR-16-5p (seed of miRNA-424-3p) and miR-34a targets, we determined that VEGFA expression was affected by all of these regulators. It was previously shown by us 3 and others 2 that VEGFA expression is reduced just prior to ovulation and increased in the luteal phase. This identified regulatory network suggests a novel cooperative mechanism for VEGF transient downregulation.
In conclusion, integrated analysis between the expression of coding genes from our preliminary study with the newly generated library of regulatory miRNAs enabled us to better understand the regulation of coding ovulatory genes by non-coding transcripts and their function from the single-molecule level to whole pathways.

Methods
Ethical approval. The current study was approved by the Sheba Medical Center Institutional Review Board (IRB) Committee (ethical approval number 8707-11-SMC), and written informed consent was obtained from each patient.
All experiments were performed in accordance with relevant guidelines and regulations.

IVF protocol.
Ovarian stimulation was carried out as previously described 57,58 according to the "short antag- For the miRNA NanoString validation experiments by qPCR, cumulus cells were obtained from 3-4 different women and pooled to generate a single replicate. Each woman donated CGCs from one M2 COC. The women's ages and infertility etiologies were similar to those described in Table S1.
For functional validation of upstream regulator analysis, cumulus cells were obtained from 3-4 different women and pooled to generate a single replicate. Each woman (total −21 women) donated CGCs of one M2 or one GV COC. Due to the limited availability of GV COC obtained from IVM, we used cumulus cells of GV COC obtained during IVF. The age of the women was 28-40 years (average 36.5 years). The etiology of infertility included male factor and PGD. IVM protocol. Four women (ages 28-40 years) undergoing IVM were selected for this study (see Supp. Table S1 for details). In the final analysis, only 3 women were included; one woman (Patient code 59 A in Table S1) was excluded because she was an outlier. Each woman donated one single COC. IVM cycles were carried out as previously described 59 . Briefly, sonographic assessment of the antral follicle count and of endometrial thickness was carried out on day 3 of a spontaneous menstrual cycle. The serum concentrations of estradiol and progesterone were also determined. Next, 150 IU/day rFSH were administered to the patients for 3 days. A second evaluation was performed on day 6. An injection of 10000 IU hCG (Pregnyl; Organon, Oss, Holland) was administered subcutaneously when the endometrial thickness was ≥5 mm and the leading follicle was at least 12 mm. Oocyte retrieval was carried out 36 hours later. Compact CCs were obtained from surplus germinal vesicle (GV) oocytes that were acquired during IVM treatment (CCGV group, provided by Dr. Ruben Fadini from the Biogenesi Reproductive Medicine Center, Istituti Clinici Zucchi, Monza, Italy).
For NanoString analysis, four women (ages 28-40 years) undergoing IVM were selected (see Supp. Table S1 for details). In the final analysis, only 3 women were included; one woman (Patient code 59 A in Table S1) was excluded because she was an outlier. Each woman donated one single COC. For the miRNA NanoString validation experiments, cumulus cells were obtained from 3-4 different women and pooled to generate a single replicate. Each woman donated CGCs from one GV COC. The women's ages and infertility etiologies were similar to those described in Table S1.
Cumulus granulosa cell collection. CGCs were obtained through oocyte denudation in the course of intracytoplasmic sperm injection (ICSI) procedures. After oocyte retrieval, CGCs of each oocyte were removed using hyaluronidase (SAGE) and a glass denudation pipette (Swemed). The CGCs were washed in PBS and centrifuged at 5000 × g for 5 minutes at room temperature. The resulting pellets were stored at −80 °C until RNA isolation. Quantitative real-time PCR. cDNA synthesis for the detection of mature miRNAs was performed with the miScript Reverse Transcription Kit (QIAGEN) using a blend of oligo(dT) and random primers (iScript ™ cDNA Synthesis Kit). cDNA synthesis for the quantification of mRNA expression was performed with the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems). miRNA expression was determined using the miScript SYBR Green Kit (QIAGEN) and mRNA expression was determined using Fast SYBR Green Master Mix (Applied Biosystems). The StepOnePlus Real-Time PCR System (Applied Biosystems) was used to detect amplification.
Expression levels were normalized to ACTB (for mRNA) and RNU6B (for miRNA). qPCR results were analyzed with StepOne software. Relative gene expression was calculated using the delta-delta Ct method 60 . Details of the primers are shown in Table S2.
NanoString and bioinformatics sample analysis. Evaluation of miRNA expression was carried out using NanoString technology, which is not based on sequencing but rather based on a digital molecular barcoding system. Each barcode is attached to a single target-specific probe corresponding to a specific miRNA. The output from this technology is the digital count of 800 tested miRNAs.
A total of 3 μl (20-120 ng) of each RNA sample was prepared according to the manufacturer's instructions. Mature miRNAs were ligated to a species-specific tag sequence (miRtag) via a thermally controlled splinted ligation. miRNA assays were performed according to the NanoString miRNA Assay Manual. Hybridizations were carried out by combining 5 µl of each miRNA assay with 20 µl of nCounter Reporter probes in hybridization buffer and 5 µl of nCounter Capture probes for a total reaction volume of 30 µl. miRNA expression was assessed using the NanoString nCounter system (NanoString Technologies, Seattle, WA, USA), which enables multiplexed direct digital counting of 800 human miRNA molecules. The raw count data from the 6 samples were analyzed using the edgeR 61 and limma R 62 packages. First, TMM normalization was applied to the data followed by voom transformation. To assess the number of miRNAs expressed under each condition, we required that probes be detectable in all replicates of the tested condition. As detectable, we defined probes having a normalized expression level that exceeded that of the "highest expressed" negative control probe. This resulted in 54 expressed miRNAs. Next, linear models were used to assess differential expression using the limma method. Pairwise comparisons were performed using moderated t statistics. To perform multiple group comparisons, one-way ANOVA was applied, except the residual mean squares were moderated between genes (see the limma package for details 62 ). miRNAs were considered differentially expressed if their FDR was <5%.
The datasets generated during and/or analyzed during the current study are available from the corresponding author upon reasonable request.
Bioinformatics analysis of the library. The resultant differentially expressed miRNAs were analyzed through the use of IPA (Ingenuity Pathway Analysis QIAGEN Inc. https://www.qiagenbioinformatics.com/products/ingenuitypathway-analysis) 19 . Putative miRNA targets found using Ingenuity miRNA target filter for experimentally validated and putative predicted targets (high confidence level) -in silico analysis.
Using whole transcriptome sequencing, we previously generated an mRNA ovulatory genes library 3 . The library was comprised of the same groups as in the current study (CCGV and CCM2), encompassing all differentially expressed mRNAs between immature mid-antral follicles and mature pre-ovulatory follicles.
Our previously generated mRNA library and the miRNA putative targets were combined using the microRNA target filter to generate a subset of the negatively correlated miRNA targets -experimental analysis.
Both mRNA target lists (in silico analysis target list and experimental analysis target list) were analyzed for GO terms and pathways using the GeneAnalytics tool 20 .
The DE miRNAs and negatively correlated targets were also analyzed using the Ingenuity Upstream Regulator Analysis tool to predict upstream molecules, including miRNA and transcription factors, that may have caused the observed gene expression changes.

Statistics.
Each experiment was performed at least three times. Data are expressed as the mean ± standard error of the mean (SEM) and were evaluated using Student's t-test with a two-tailed distribution, with two samples equaling variance, or with ANOVA (ANalysis Of VAriance) for more than two variances using the post hoc Tukey test assuming equal variance, or the Games-Howell test for unequal variance. When appropriate, the Kruskal-Wallis non-parametric comparison test was used. For all statistical analyses, SPSS 22 software (IBM, Armonk, NY, USA) was used. P-value < 0.05 was considered statistically significant.