Transcriptomics of cumulus cells – a window into oocyte maturation in humans

Background Cumulus cells (CC) encapsulate growing oocytes and support their growth and development. Transcriptomic signatures of CC have the potential to serve as valuable non-invasive biomarkers for oocyte competency and potential. The present sibling cumulus-oocyte-complex (COC) cohort study aimed at defining functional variations between oocytes of different maturity exposed to the same stimulation conditions, by assessing the transcriptomic signatures of their corresponding CC. CC were collected from 18 patients with both germinal vesicle and metaphase II oocytes from the same cycle to keep the biological variability between samples to a minimum. RNA sequencing, differential expression, pathway analysis, and leading-edge were performed to highlight functional differences between CC encapsulating oocytes of different maturity. Results Transcriptomic signatures representing CC encapsulating oocytes of different maturity clustered separately on principal component analysis with 1818 genes differentially expressed. CCs encapsulating mature oocytes were more transcriptionally synchronized when compared with CCs encapsulating immature oocytes. Moreover, the transcriptional activity was lower, albeit not absent, in CC encapsulating mature oocytes, with 2407 fewer transcripts detected than in CC encapsulating immature (germinal vesicle - GV) oocytes. Hallmark pathways and ovarian processes that were affected by oocyte maturity included cell cycle regulation, steroid metabolism, apoptosis, extracellular matrix remodeling, and inflammation. Conclusions Herein we review our findings and discuss how they align with previous literature addressing transcriptomic signatures of oocyte maturation. Our findings support the available literature and enhance it with several genes and pathways, which have not been previously implicated in promoting human oocyte maturation. This study lays the ground for future functional studies that can enhance our understanding of human oocyte maturation.


Background
Cumulus cells (CC) provide somatic support to the maturing oocyte, and together they comprise the functional unit known as Cumulus-Oocyte-Complex (COC) [1][2][3][4]. Understanding oocyte maturation and associated pathologies can help improve fertility treatments, as well as culture media conditions in the lab. Cumulus cells can be collected without compromising the oocyte, and their transcriptomic signatures are valuable non-invasive biomarkers for processes within the oocyte [5,6].
Oocyte maturation is contingent on rapid transcription and translation, governed by paracrine and autocrine signaling prior to ovulation [7][8][9]. Once matured, the MII oocyte is less transcriptionally active than its precursors, relying on stored mRNA transcripts that were acquired throughout its maturation, to undergo successful fertilization and support early embryo development until embryonic genome activation [10][11][12]. Moreover, in addition to the stored transcripts, there is active transportation of transcripts from cumulus cells to the growing oocyte through trans-zonal projections [13,14]. These projections are critical for both oocyte and cumulus cell differentiation [15].
Currently, oocyte assessment relies mainly on morphological criteria that provide little insight on oocyte quality and competence [16]. Furthermore, available techniques for maturing oocytes in-vitro are inefficient and do not provide good alternatives in cases were pathologies of oocyte maturation lead to retrieval of multiple immature eggs despite adequate controlled ovarian stimulation. This is why molecular investigation of processes in COCs responsible for nuclear and cytoplasmic maturation, as well as pathologies that could arise, are key in improving patient treatment and outcomes. The objectives of this study were: 1) to profile the transcriptome of CC from mature MII and immature GV oocytes, from the same treatment cycle, 2) to use the above profile to validate existing transcriptomic literature exploring human oocyte maturation, and 3) to provide a comprehensive list of genes impacted by in vitro maturation which can be used for future studies exploring human oocyte maturation.

Collected samples and patient demographics
A total of 22 cumulus-oocyte complexes (COCs) were collected from 11 patients, for RNA sequencing, with a mean age of 33.9 ± 1.6 years old, mean BMI of 25.4 ± 2.3 kg/m 2 , mean Anti-Mullerian Hormone (AMH) levels of 29.8 ± 12.5 pmol/L, mean Day 3 FSH levels of 7.9 ± 1.4 mIU/mL, and the mean number of collected oocytes was 10.6 ± 1.9 (Table 1). All samples had sufficient number of sequencing reads, high average quality scores, and high sequence alignment rates sufficient for differential gene expression analysis, as per guidelines previously published for quality control of RNAseq experiments [17] (Supplemental Table S2).
Samples clustered according to the degree of maturity of the encapsulated oocyte A total of 6220 genes were detected in the MII-CC cohort, 202 of which were unique, and 8627 genes were detected in the GV-CC cohort, 2609 of which were unique (Fig. 1a). Unsupervised hierarchical clustering demonstrated a separation of the GV-CC cohort from the MII-CC cohort. Notably, the MII-CC cohort clustered more tightly than the GV-CC cohort, indicating that with maturation there is decreased inter-sample variability (Fig. 1b). This was further demonstrated by principal component analysis (PCA) (Fig. 1c) with 16.3% of variability in the dataset corresponding to oocyte maturity (PC1).
Differential expression reveals marked differences in gene expression according to degree of maturity of the encapsulated oocyte In order to ensure the sample size and sequencing depth were sufficient to capture biologically significant differences between the MII and GV CCs, a post-hoc power analysis was conducted showing that the minimum required biological replicates is 6 and the minimum sequencing depth is 10 million reads per replicate, both of which were exceeded in this study [18]. When comparing MII-CC with GV-CC cohorts, 1818 genes were differentially expressed (Supplemental Table S3), which comprise 10.3% of annotated RefSeq genes (2 < FC < − 2 and FDR < 0.05). Of these, 40 genes changed by 10-fold or more, 207 genes changed by 5 to 10-fold, and the rest (1571 genes) changed by 2 to 5-fold. When accounting for the direction of differential expression, 1031 genes were significantly downregulated (5.9% of annotated genes) and 787 genes were significantly upregulated (4.5% of annotated genes) in the MII-CC cohort compared with GV-CC (Fig. 1d). The top 20 genes enriched in the different maturation stages are reported in Table 2, all differentially expressed genes are reported in Supplemental Table S3.
Novel findings from this study enhance available literature exploring processes that lead to synchronized oocyte maturity When comparing our differentially expressed genes to those previously reported in the literature, 42 genes were associated with oocyte maturation in both this study and in previous literature (Table 3). Forty-five genes previously correlated with oocyte maturation were not differentially expressed in the current study (Supplemental Table S4) [27]. Three thousand five hundred and fifty-four genes were differentially expressed in our study and have not been annotated previously in studies exploring oocyte maturation. Of these, 129 are known to be regulated by at least one gonadotropin, 82 have been previously reported to be regulated by LH alone, 16 by FSH alone and 31 by both LH and FSH (Supplemental Table S5).
Pathway analysis and leading-edge analysis revealed significant differences between maturational stages Gene set enrichment analysis (GSEA) was conducted to determine the pathways and cellular processes that are altered throughout oocyte maturity, thus allowing for the interpretation of the complex interactions between differentially expressed genes ( Fig. 2) [28]. Furthermore, leading edge analysis (LEA) helped determine which genes were driving pathway enrichment scores, and to better understand the major biological differences between the two cohorts.
Pathways and processes that were primarily enriched in downregulated genes in the MII-CC cohort included Nuclear maturation, chromatin remodeling and replication initiation, faithful chromosome segregation, apoptosis and inflammation (Fig. 2). Furthermore, specific genes identified as biologically significant genes for oocyte maturation including nuclear maturation (Fig. 3a), chromatin remodeling and DNA replication initiation (Fig. 3b), and apoptosis and inflammation (Fig. 3c) are further highlighted in Fig. 3.
Pathways and processes that were primarily enriched in upregulated genes in the MII-CC cohort included Extracellular matrix (ECM) components and its remodeling enzymes, and steroid metabolism and processing (Fig. 2). Genes identified as biologically significant genes for oocyte maturation including ECM remodeling (Fig. 3d), and steroid metabolism (Fig. 3e) are further highlighted in Fig. 3.
Leading edge analysis identified several genes involved in cell cycle control (CDK1, CCNB1, CCNB2, CCNA2, BUB1, and CDC20), DNA replication initiation (MCM2-7), and centromere assembly and organization (CENPF), among others, as the genes that were most significantly driving the gene set enrichment analysis (Supplemental Table S6). GSEA on DE genes known to be regulated by at least one gonadotropin revealed two major pathways; transcriptional regulation of tp53 (apoptosis), overall enriched in downregulated genes, and metabolic biosynthesis, overall enriched in upregulated genes. LEA identified several genes involved in cell cycle control and cell death indicating that the MII-CC cohort, in response to LH and FSH decreases cell death signaling and increases biosynthesis.

Validation of NGS results by qPCR
We selected 16 differentially expressed genes as determined by RNAseq and are known to be involved in various pathways of CC expansion and oocyte maturation. For all selected genes, similar fold changes were observed using qPCR as were observed using NGS (Fig. 4a). In addition, when assessing multiple CCs from 7 additional patients at both the MII and GV stage (from 2/3 GV-CC and 2/3 MII-CC per patient) (a total of 34 COCs), the expression of all tested genes was consistent within each patient. Interestingly, the expression of all genes used for validation across all MII-CC were consistent both within and between patients (ΔCt SEM (range) = 0.26(0.12-0.56)). GV-CC were also consistent within and between patients, however to a lesser extent (ΔCt SEM (range) = 0.40(0.21-0.57)) ( Fig. 4b).
This was further validated using Pearson's R correlation which also demonstrated that the cohort with the highest similarity was MII-CC within the same patient, with a R 2 of 0.954 ± 0.023 followed by GV-CC within the same patient, with a R 2 of 0.906 ± 0.036. The correlation coefficients and significance are outlined in Table 4. Taken together, this indicates that the heterogeneity of the population of CCs within the ovary is low, and sampling one mature MII-CC and one immature GV-CC for NGS is representative of the cohort of CCs.

Discussion
This study is novel in the choice of cohorts for comparison. We included CC encapsulating oocytes arrested at the GV stage, despite being exposed to adequate stimulation, and compared them with CC encapsulating MII oocytes that matured in vivo, from the same patients during an IVF cycle. This comparison allowed us to fine-tune our understanding of oocyte maturation invivo. Previous human oocyte maturation studies analyzed COCs from in-vitro maturation cycles [16,[29][30][31][32][33]. This is why their findings are more relevant for  processes in in-vitro maturation per se, but they do not tease out factors specifically associated with failed maturation despite adequate in-vivo exposure to controlled ovarian hyper stimulation (COH). Forty-two differentially expressed genes in our study, have been previously associated with oocyte maturation and cumulus cell expansion in IVF treatments (Table 3). Following our extensive literature search, we identified 45 genes, which were previously correlated with human oocyte maturation, but not captured by our study design, possibly because they may impact IVM alone and may not reflect in vivo maturation (Supplemental Table S4). A third group of 3554 genes captured in the current study but not in previous studies which represents a novel group of genes that should be further explored as they have not been previously implicated in human oocyte maturation (Supplemental Table S3). Of these, 129 genes have been previously shown to be regulated by either LH, FSH, or both (Supplemental Table S5).
In this study, several factors and their regulators involved in nuclear maturation and cell cycle control were differentially expressed between cumulus cells encapsulating oocytes of different maturity, reiterating findings from previous studies [19,20,29,34]. These include cell cycle regulators (BIRC5, BUB1, BUB1B, CCNA2, CCNB, CDK1, FBXO5 MAD2L1, and PTTG1) and components of the centromere (CENPA, CENPE, and CENPH) [29]. In our MII-CC cohort we observed downregulation of MCM2-7, which form the hexameric pre-replication protein complex. This complex is involved in initiating replication forks and recruiting other DNA replication related proteins. We also observed downregulation of TOP2A, which relaxes supercoiled and circular DNA molecules. Reinforcing available literature that states that while crucial at the MI stage for chromatin remodeling [21,22], its activity decreases in mature oocytes [23].
Apoptosis was also attenuated in the MII-CC cohort, further supporting decreased cell turnover with advanced maturity. Related pathways including Wnt pathway and Akt-pathway were affected, as demonstrated by downregulation of SFRP4, a potent inhibitor of Wnt signaling [24], and upregulation of OSMR, an activator of Akt-mediated proliferation [25]. These findings corroborate previous literature reporting downregulation of SFRP4 during oocyte maturation [26,35,36], and upregulation of OSMR in bovine preovulatory follicles posttriggering by gonadotropins [37].
Extracellular matrix remodeling was also altered between the two maturity cohorts, as evident by members of the matrix metalloproteinases (MMP) family and their inducers (MMP11 and SPARC1L). Again, this supports previous literature showing significant decrease of MMP11 in granulosa cells following hCG administration [38]. This effect is further demonstrated by increased expression of TNC, NID2, and SPOCK2 -all ECM proteins and MMP substrates [26,[39][40][41]. Notably, well characterized ECM remodeling enzymes, ADAMTS1 and SERPINE2, were also differentially expressed, aligning with previous studies [42,43]. Both play critical roles in follicular remodeling during follicular growth and rupture [44], by metabolizing Versican and Hyaluronan which lead to cumulus cell matrix expansion and attenuation [45].
Another key process enhanced in follicular niche maturation is inflammation, which is crucial for ovulation. Upon gonadotropin stimulation, the follicle wall is weakened, thereby facilitating its eventual rupture [46]. In our MII-CC cohort, we observed marked upregulation of genes associated with inflammation, including members of the Interleukin and TGF-beta families. Among the genes upregulated in our MII-CC cohort were IL18R1 which promotes cumulus cell expansion [47], and TGFBR3 which promotes cellular differentiation, migration, adhesion and extracellular matrix production [48,49]. IL6ST which is part of the cytokine receptor complex (gp130) was also upregulated in the MII-CC cohort, consistent with previous studies in non-human primates and equine models [50,51].
Key players that emerged in our cohort as being significant for cumulus cells to facilitate oocyte maturation are AREG, EREG, PTGS2, and STAR. Two factors at the heart of this complex process are AREG and EREG, which have been shown to mediate the LH signal driving cumulus expansion and oocyte maturation [19,32,52]. They also activate the EGF receptor (EGFR) which in turn releases matrix metalloproteinases (MMPs) and promotes cumulus expansion [52,53]. Furthermore, in conjunction with progesterone, AREG and EREG Fig. 2 Pathway analysis of significantly differentially expressed genes. GSEA (Gene Set Enrichment Analysis) reveals that the MII-CC cohort significantly downregulate pathways involved in chromatin assembly, apoptosis, mitotic cell cycle control, and DNA repair processing (in blue) and significantly upregulate pathways involved in lipid biosynthesis, steroid metabolism, inflammation, and leukocyte activation (in orange). A total of 60 gene sets were enriched in upregulated genes and 223 gene sets were enriched in downregulated genes at FDR q-value < 0.05. The size of the node corresponds to the number of genes in each gene set enhance PTGS2 (also upregulated in our MII-CC cohort) via EGF to increase prostaglandin production and maintenance of chromosomal spindles [33,[54][55][56]. In addition, AREG mediates hCG-induced STAR expression (also upregulated in our MII-CC cohort), which plays a key role in steroid and progesterone production in human granulosa cells [57], and is a potential predictive biomarker for nuclear maturation [58] and oocyte quality [33]. It is important to note, that despite being well defined as key in ovarian maturation [32,52,59], EREG has not been found to be differentially expressed in previous genomic signature studies addressing this question. This further highlights the importance of our study design in better refining the pathophysiology of oocyte maturation.
Other critical factors confirmed by this study are IL1, FSHR, BDNF, HSD11B1, and HSD17B1, all of which are implicated in the control of steroid synthesis and paracrine response to steroids. IL1 (both alpha and beta subunits), which stimulates steroidogenesis, was upregulated in the MII-CC cohort with a concurrent decreased expression of FSHR in the same cohort, substantiating what was previously observed in rodents and humans [60,61]. BDNF, which modulates granulosa cell function via FSHR-coupled signaling pathway, to affect aromatase-mediated steroidogenesis, was also downregulated in our MII-CC cohort [62].
HSD11B1, the enzyme responsible for cortisone production, an essential substrate for steroid hormone synthesis, was upregulated in our MII-CC cohort. A companion enzyme, HSD17B1, catalyzes the last step in estrogen metabolism converting E1 of low estrogenic  activity to E2 of high activity using cortisone as a substrate [63]. HSD17B1 has not been captured in previous human studies, but was downregulated in our MII-CC cohort, consistent with the results seen in a previous bovine study [64], and further highlighting the advantage of our study design. Both LHCGR and FSHR were differentially expressed in the MII-CC cohort when compared to the GV-CC cohort. To further explore how this may impact the transcriptional profiles, we identified all genes differentially expressed between these two cohorts that are known to be regulated by LH and/or FSH. GSEA based on these 129 genes, identified two major pathways: regulation of apoptosis (via p53) and biosynthesis of various biomolecules.
Overall, apoptosis was enriched in downregulated genes. Interestingly, several major players in the regulation of apoptosis, including BIRC5, TP53, HMGB1, HMGB2, and SFRP4 are also known to be regulated by LH and/or FSH [24,35,[65][66][67]. Taken together, these findings suggest that the inability of the GV-CC cohort to appropriately respond to FSH/LH may, in turn, lead to failure of the CCs to effectively dampen the apoptotic signals, causing the COC to enter a stage of maturation arrest and follicular atresia.
Overall, biosynthesis was enriched in upregulated genes among the MII-CC cohort. Notably, several members of the CYP family, which were upregulated, and are involved in the biosynthesis of estrogen and androgens, are known to be regulated by LH and/or FSH [68][69][70]. Taken together, these findings suggest that the MII-CCs are responding adequately to gonadotropin administration.
Collectively, it appears that the GV-CC cohort, failed to adequately synthesize estrogen, despite exposure to gonadotropins and, thus, began to upregulate genes involved in apoptosis. This may be due to insufficient LH and/or FSH receptors on these COCs or due to another underlying malady.
Finally, we show that PDE3A, known to improve nuclear-cytoplasmic synchrony [71], is significantly upregulated in our MII-CC cohort. While this gene has not been studied in cumulus cells in the context of oocyte maturation in humans, it has been shown that an increase in oocyte PDE3A activity causes delayed spontaneous meiotic maturation, coupled with extended gap junctional communication between the CC and the oocyte. Such a delay has a positive effect on oocyte cytoplasmic maturation, thereby improving oocyte developmental potential [72]. The fact that upregulation of this gene was captured by our study design speaks once again to the strength of our study and to what it adds to current literature.
Methodological strengths of this study include (i) a sibling COC design allowing to minimize the biologic variability between cohorts, (ii) exploring transcriptomic dynamics in cumulus cells, which are considered valuable non-invasive markers for oocyte quality [73][74][75], and (iii) performing next generation sequencing (NGS), which is the most unbiased approach currently available for exploring transcriptomic signatures. A methodological weakness of this study is our inability to compare our findings with a third cohort of CCs encapsulating naïve GVs from the same stimulation cycle. Furthermore, to keep our sequencing sample size small, only 2 cumulus masses were selected from each patient; 1 to represent the immature GV COCs, and the other to represent the mature MII COCs. To lower the risk for selection bias due to this study design, we chose to perform validation studies that established the tight correlations within each cohort. Lastly, our small sample size and study design did not allow for tracking outcomes on an individual oocyte basis.

Conclusions
In conclusion, our findings enhance current literature on oocyte maturation by identifying CC genes not previously associated with oocyte maturation that may be involved in this process. Our novel list of genes can serve as a springboard for future studies. Our future studies will focus on determining the functional significance of these findings and on attempting to identify how different treatment options may favor a more synchronized mature/competent state. In addition, to further validate genes that are critical for oocyte maturation and competency, further large-scale studies correlating gene expression with clinical outcomes using a targeted transcriptome panel are needed.

Patient recruitment, data collection and cumulus cell isolation
Cumulus cell samples were collected from eighteen patients undergoing IVF-ICSI cycles at the CReATe Fertility Centre (Toronto, ON, Canada), between August 2016 and June 2017. Exclusion criteria were patients diagnosed with PCOS, as per Rotterdam criteria, as well as patients with endometriosis diagnosed by laparoscopy. Samples from eleven patients were used for RNA-seq (22 COCs, 11 mature (MII), and 11 immature (GV)), and samples from seven additional patients were used for qPCR validation of the findings (a total of 18 mature (MII) and 16 immature (GV) COCs). Patients were treated using a standard antagonist protocol, with initial gonadotropin dosing and subsequent adjustments at the discretion of the treating physician.
Ultrasound guided oocyte retrieval was performed 35-36 h post hCG injection. COCs were identified under a stereomicroscope and only COCs completely and tightly enclosed by compact CCs were used for this study to minimize the potential collection of contaminating granulosa cells. Selected COCs were serially washed three times in Quinn's Advantage Medium (Sage, USA) to remove cellular contaminants, and to further reduce the possibility of granulosa cell contamination. CCs were then mechanically separated from each oocyte individually in Quinn's Advantage Medium (Sage, USA), under paraffin oil by one experienced embryologist within 1 h of oocyte retrieval. The oocytes corresponding to individually collected CC were separately exposed to hyaluronidase (80 IU/ml) immediately after mechanical separation of CCs, washed in Quinn's Advantage Medium (Sage, USA). Maturational stage was assessed through the observation of the nucleus of the oocyte. Oocytes with an extruded polar body were deemed mature (MII), oocytes with an intact germinal vesicle was deemed immature (GV), oocytes without an observable germinal vesicle or an extruded polar body were deemed MI and excluded from further analysis. The CC were collected from single MII (n = 29) and GV oocytes (n = 27), frozen separately in 300ul of RNA lysis buffer RL (Norgen Biotek, Canada), and stored at − 80°C until RNA extraction. Clinical data including patient demographics, medical history, and ovarian stimulation related parameters, were collected for all patients enrolled in this study.

RNA extraction and NGS library preparation
RNA extraction, cDNA conversion and NGS library construction and normalization were conducted as previously described [76]. Briefly, for NGS one MII-CC and one GV-CC were chosen at random from each of the 11 patients for RNA extraction using the Total RNA Purification Kit Micro (Norgen Biotek, Canada). The quantity was assessed using Qubit RNA HS (ThermoFisher, USA) and RNA integrity assessed using 2100 Bioanalyser RNA 6000 Pico Total RNA Kit (Agilent Technologies, Canada). cDNA was synthesized using the SMART-seq v4 Ultra Low Input RNA Kit (Takara Bio Inc., Japan) according to the sample preparation guide and using 14 rounds of amplification. Sequencing libraries were constructed using Nextera XT (Illumina, USA) and 1 ng of amplified cDNA according to the sample preparation guide. Final sequencing libraries were assessed for quantity and quality using the KAPA Library Quantification Kit (KAPA Biosystems, Switzerland) and 2100 Bioanalyser High Sensitivity DNA Kit (Agilent Technologies), respectively. Normalized libraries were pooled, denatured, diluted to 1.4 pmol/l and loaded onto a High Output (300 cycle) flow cell (Illumina, USA) followed by sequencing (2 × 127 bp) on a NextSeq 550 (Illumina, USA).

Differential expression
The recommendations outlined by Ching et al. 2014 were followed when selecting the differential expression package, as well as using a paired-sample RNA-Seq design as suggested [17]. FASTQ files were generated using bcl2fastq2 (v2.17) and the read quality was assessed. Sequences were trimmed based on quality (Phred > 28). Raw trimmed reads were aligned to Human Genome Assembly 38 (hg38) using STAR (v2.5.3a) [77] and quantified to RefSeq (Release 84). Low expressed transcripts were excluded (maximum counts < 10) and differential expression (DE) was conducted on the remaining counts using DESeq2 (v3.5) [78]. Principal Component Analysis (PCA) and hierarchical clustering (HC) were conducted to assess the relationship between samples and determine covariates contributing to variation in the dataset. Principal component 1 (PC1) accounts for the largest proportion of the variability observed within the dataset, PC2 accounts for the second largest, and so on. Differentially expressed genes were identified by comparing all mature CC samples (MII-CC) to all immature CC samples (GV-CC), and were deemed to be differentially expressed if the gene had a Fold change (FC) of more than absolute value of 2, and a false discovery rate (FDR) < 0.05. This analysis was conducted in Partek Flow (version 8.0.19.0408).

Pathway analysis
Gene Set Enrichment Analysis (GSEA) was performed to determine the effect all differentially expressed genes have on cellular processes and functions [28]. The resulting pathway list was cross referenced with a custom gene set created and supported by the Bader Lab (University of Toronto) which is comprised of all GO database, KEGG, and Reactome gene sets (v2018-12-01) (http://download.baderlab.org/EM_Genesets/) [79]. Genes that could not be mapped to any gene-set term were excluded from the comparison. Gene sets with 10 or fewer genes and/or a q-value > 0.05 were excluded from further analysis. Following GSEA, leading edge analysis (LEA) was conducted to determine what genes were driving the gene set enrichment score, as well as to highlight genes that were shared between gene sets.
To further explore the impact FSH and/or LH may have on the transcriptome, we identified all differentially expressed genes that are known to be regulated by LH, FSH or both [80] and performed GSEA and LEA as described previously.

NGS validation by qPCR
Sixteen genes (and one reference gene) were chosen for validation from the list of differentially expressed genes. The choice of genes was based on previous annotations deeming these genes as biologically significant, identification by leading edge analysis, and/or participation in key ovarian gene pathways. Pre-designed and validated PrimeTime™ qPCR assays (IDT, USA) were used for validation of NGS results with RPLP0 as the reference gene. All targets were assayed in duplicate using PrimeTime™ Gene Expression MasterMix (IDT, USA) (polymerase activation at 95°C for 3 min; 45 cycles of 15 s denaturation at 95°C and 1 min annealing/extension at 60°C). Relative fold change (ΔΔCt) was employed to quantify gene expression [81]. Data analysis was performed using GraphPad Prism (version 5.02). The list of primers and probes used for validation are given in Supplemental  Table S1. qPCR for the above validated genes was also performed on multiple COCs of the same maturational stages on samples from 7 additional patients. This was carried out to ensure the validity of a random choice of a single COC as a representative of all COCs at the same maturational stage from the same patient. Similarities between individual CCs was measured by creating an overall Pearson correlation between the vectors of variables (16 targets genes) using SPSS Proximities. This created an overall Pearson R 2 . Paired t-test was used to compare correlations within GV-CC and within MII-CC groups, separately. A non-paired test was used to compare overall GV-CC and overall MII-CC correlations due to an unequal number of observations for some individuals.

Gene annotation and literature search
To determine the clinical significance of our bioinformatic findings, differentially expressed genes were cross referenced with available datasets in the literature by searching the PubMed database for previous studies assessing the transcriptome of human CC using NGS, Microarray, or qPCR. Differentially expressed genes were further reviewed in depth using the Ovarian Kaleidoscope Database [80] and GeneCards Human Gene databases (http://www.genecards.org/), to correlate our bioinformatic findings with hallmark physiological and pathological processes in the ovary.
Additional file 1 : Supplemental Table S1: List of primers and probes used for qPCR validation of RNAseq results. Supplemental Table S2: RNAseq quality control metrics for each cohort presented as mean ± SEM (Range). Supplemental Table S3: All differentially expressed genes between MII-CC and GV-CC cohorts. Supplemental Table S4: Differentially expressed genes known to be regulated by LH, FSH, or both with their associated fold change, and number of gene sets that are found in following leading edge analysis. Supplemental Table S5: Previously identified genes implicated in oocyte maturation in vitro that were not differentially expressed in this study. Supplemental Table S6: Top 50 genes identified by leading edge analysis ranked by the number of unique gene-sets the gene was identified in. Supplemental Table S7: Normalized read counts of all samples.