Next-Generation Sequencing Reveals Downregulation of the Wnt Signaling Pathway in Human Dysmature Cumulus Cells as a Hallmark for Evaluating Oocyte Quality

Abstract

Background: Dysmature cumulus cells are lower fertilization rates and abnormalities in embryonic development compared to maturation cumulus cells. Morphological evaluation of cumulus–oocyte complexes (COCs) considered the possibility that differences may also be found in gene expression. Purpose: To identify hallmarks for evaluating oocyte quality by investigating gene expression patterns in human cumulus cells surrounding oocytes. Methods: Cumulus cells were obtained from the cumulus–oocyte complex of infertile women treated with assisted reproductive technology. Based on maturity level, the cumulus cells were classified into two categories, i.e., dysmature cumulus cell (DCC) and maturation cumulus cell. DCCs were subjected to gene expression analysis using next-generation sequencing and compared with COCs that are in the process of maturation as controls. Results: The expression levels of genes involved in the Wnt signal/β-catenin pathway were significantly reduced in DCCs compared with those in control cells. Moreover, the expression levels of genes involved in multiple pathways associated with apoptosis were also significantly reduced compared with those in control cells. Conclusions: DCCs showed significant decreases in apoptosis- and Wnt/β-catenin signaling-associated gene expression. DCCs could be classified by morphological evaluation, and the method described in this study may be useful as an oocyte quality estimation tool.

1. Introduction

A human oocyte and its supportive cells, known as granulosa cells, are interdependent throughout ovarian development and ovulation for successful fertilization [1]. Oocytes isolated from primordial follicles fail to grow without the aid of granulosa cells [2,3]. Cumulus cells are the innermost subgroup of granulosa cells surrounding the oocyte and play a key role in oocyte maturation and fertilization. Cumulus cell function depends on the gap junctions between cumulus cells and the oocyte, which allow rapid transfer of molecules to the oocyte [4,5]; so, the quality of cumulus cells is thought to reflect oocyte quality. Therefore, transcriptomic profiling of cumulus cells is useful in identifying biomarkers indicating the developmental ability of oocyte maturation (Figure 1a).

The classification of mature cumulus–oocyte complexes (COCs), based on morphology, was first described in the 1980s [6] and has been used for many decades [7,8,9,10].

In the human menstrual cycle, only less than 1% of follicles lead to ovulation, and more than 99% of follicles lead to atretic follicles [11]. We have used these criteria for more than a decade and reported that the oocytes classified as “dysmature” exhibit a remarkably lower fertilization rate, compared to COCs that are in the process of maturation [12]. Our previous studies revealed significant differences in gene expression of DCCs compared to COCs that are in the process of maturation [13], and DCCs are reportedly associated with apoptosis [13] and autophagy-independent cell death [1,12].

To date, several studies using next-generation sequencing (NGS) to examine cumulus cells or granulosa cells have been conducted. Using NGS, it is possible to comprehensively investigate gene expression in samples, because it can determine millions of sequences at once. For example, Hu et al. [14] reported the cause of insulin resistance in polycystic ovarian syndrome (PCOS) analyzing cumulus cells. Xu et al. [15] searched the difference in gene expression of mouse cumulus cells with and without blastocyst formation. Fuchs et al. [16] studied the difference in the gene expression of human cumulus granulosa cells based on the types of drugs used to trigger IVF cycles. Although other studies using NGS for cumulus cells are also present [17,18], there have been a few studies on the gene expression related to COC maturity.

As described above, the evaluation of COC maturation is useful to estimate oocyte quality. Furthermore, the procedure used to assess cumulus cells is more straightforward than assessing the oocyte without damaging it. Nevertheless, there have been a few studies that estimate the difference in gene expression based on the maturity of COCs.

In the present study, we compared DCCs, which would become atretic follicles, with COCs that are in the process of maturation, and used NGS to search for candidate genes associated with follicular development, lower fertilization rate, apoptosis, and autophagy.

2. Materials and Methods

2.1. Patients

This study was reviewed and approved by the ethical committee of the National Center for Child Health and Development (# 581, 3 July 2012). Women presenting unexplained infertility aged 33–42 years who underwent assisted reproductive technology (ART) were included. Briefly, cumulus cells were collected from patients who consented to participate in this study. Based on the morphological criteria, DCCs were selected, and COCs that are in the process of maturation were categorized as control (Supplemental Table S1). Twenty samples of cumulus cells were collected from patients who provided consent for their samples to be used for research purposes.

The patients underwent standard in vitro fertilization (IVF) treatment, which included gonadotropin administration to induce multiple follicular development, and transvaginal oocyte aspiration. COC retrieval was performed by vaginal puncture under ultrasound echo-guidance 36 h after human chorionic gonadotropin (hCG) administration (10,000 IU). The classification of COCs was based on the following morphological criteria: We defined DCCs as the presence of dark spots in the layers of cumulus cells, with or without complete expansion with a visible halo of the cumulus layer, to distinguish them from maturation cumulus cells (Figure 1b). The cumulus cells were mechanically removed at 40 h after hCG administration, before conventional IVF or intracytoplasmic sperm injection, washed in Sydney IVF fertilization (COOK Medical, Brisbane, Australia) medium, and frozen individually before total RNA extraction. The sample collection flow is shown in Figure 2.

2.2. RNA Extraction and NGS

Total RNA of cumulus cells from oocytes were extracted using the RNeasy mini kit (Qiagen, Valencia, CA, USA) following the manufacturer’s protocol. The extracted RNA quality and quantity were assessed on a 2100 Bioanalyzer, using the RNA 6000 Nano kit (Agilent Technologies, Palo Alto, CA, USA). Libraries were prepared using a TruSeq RNA Sample Preparation Kit (Illumina Inc., San Diego, CA, USA) to generate clusters on the Illumina cBot (Illumina). Sequencing was performed on the Hiseq 2000/2500 (Illumina) as 50-bp reads, while image analysis and base calling were performed with CASAVA ver.1.7.0/1.8.2 (Illumina), following the manufacturer’s instructions. High-quality sequences, which passed the default quality-filtering parameters in the Illumina pipeline GERALD stage were retained and exported into Strand NGS (Agilent Technologies) for subsequent data analysis. The next-generation Web sequence data were deposited in the Japanese Genotype-phenotype Archive (JGA) [Accession #: JGAS00000000212].

Sequences were aligned with the human genome sequence (Feb. 2009/hg19), using the gene annotation file downloaded from the UCSC sequence and annotation database (http://hgdownload.cse.ucsc.edu/downloads.html#human), and Strand NGS (Agilent Technologies). The downloaded UCSC annotation database included 40,437 genes and transcripts consisting of 19,805 protein-coding genes, 12,221 pseudogenes, and 7365 non-coding RNA (Supplemental Table S2).

2.3. Gene Expression Profiling Analysis

Uniquely aligned reads were quantified, and both gene and exon-level reads per kilobase per million mapped read (RPKM) values were calculated. The RPKM values represent counts of reads mapping a feature (gene, or exon), normalized to both overall sequencing coverage and the size of the feature [19].

Before gene expression analysis and the use of log-transformed and base-lined counts, we performed principal component analysis (PCA) and hierarchical clustering to overview the expression profiles of all samples. The values were normalized to mean 0 and variance 1 before applying log transformation.

Differentially expressed genes in DCCs were identified based on a 2-fold change compared to the mean value of the control, and statistical significance based on p-value < 0.05, using one-way analysis of variance (ANOVA), followed by post hoc analysis using Student–Newman–Keuls. The list of differentially expressed genes was utilized for gene ontology (GO) analysis for biological/molecular function, using Strand NGS with a false discovery rate-corrected p-value < 0.01. Furthermore, the list of differentially expressed genes was analyzed for pathway enrichment analysis using WikiPathway, and statistical significance was set at a p-value < 0.05.

2.4. Immunofluorescence

Cumulus cells prepared as described above were permeabilized with 0.5% Triton X-100 in phosphate-buffered saline (PBS) for 5 min at 25 °C and blocked with 1% BSA and 0.1% Tween 20 in PBS for 1 h at 25 °C. The samples were incubated overnight at 4 °C with monoclonal anti-β-Catenin-Cy3 conjugate (C7738, Sigma, St Louis, MO, USA), prepared at a 1:500 dilution in Dako antibody diluents (Dako, Carpinteria, CA, USA). The complexes were then washed three times with 0.1% Tween 20 in PBS. The nuclei were then counterstained with 4′,6-diamidino-2-phenylindole (DAPI) at a final concentration of 10 μg/mL in 0.1% Tween 20 in PBS for 30 min at 25 °C. Thereafter, samples were washed three times with 0.1% Tween 20 in PBS, immersed in Vectashield (Vector Laboratories Inc., Burlingame, CA, USA), and covered with coverslips. Microphotographs were captured using a confocal laser scanning microscope LSM510 Meta (Carl Zeiss MicroImaging, Inc., Thornwood, NY, USA). Photographs captured from three different areas of each sample were analyzed using Olympus DP Controller and DP Manager (Olympus, Tokyo, Japan). In each experiment, at least 300 cells were measured in randomly selected microscopic fields.

3. Results

3.1. Sequencing and Mapping Summary

Sequencing data for each sample resulted in approximately 10–27 M reads (Supplemental Table S3). Among the raw sequence reads, 77–85.5% aligned with the human genome sequence (Supplemental Table S3). Read counts were converted to RPKM [19] to assess global gene expression profiles. Based on the gene expression patterns, PCA was carried out, followed plotting 20 samples that were divided into DCCs and control. The results showed that samples collected from DCCs were differentially located compared to the control (Supplemental Figure S1a).

3.2. Differentially Expressed Genes

To determine the statistical difference in gene expression between DCCs and the control, we performed a one-way ANOVA followed by pair-wise comparisons. Among genes with a statistical significance, we selected those with two-fold differential expression. The Venn diagrams show genes with both two-fold differential expression and a statistical significance threshold, analyzed using a one-way ANOVA. Among these, 165 genes were upregulated, while 363 were downregulated (Supplemental Figure S1b, Supplemental Table S4).

3.3. GO and Pathway Analysis

We performed GO analysis to examine the enrichment of biological processes, molecular function, and cellular components and identify biological characteristic for differentially expressed gene list. Upregulated genes were enriched in four functions: single-organism processes, cellular processes, metabolic processes, and response to stimuli. Whereas in downregulated genes, only the plasma membrane portion was enriched in cellular components (Supplemental Table S5).

In upregulated genes, biological processes such as mitotic prometaphase, mitotic metaphase and anaphase, cell cycle checkpoints, kinesins, and mitotic G1-G1-S phases were affected significantly (all, p < 0.0001) (Supplemental Table S6). In contrast, the Wnt signaling pathway and pluripotency, the integrated pancreatic cancer pathway, platelet homeostasis, arachidonic acid metabolism, and the netrin-1 signaling pathway were significantly affected in downregulated genes (all, p < 0.01) (

Table 1. Downregulated pathways in DCCs.

Pathway	p-Value	Matched Entities	Pathway Entities of Experiment Type
Hs_Wnt_Signaling_Pathway_and_Pluripotency_ WP399_74897	0.0000586	5	100
Hs_Integrated_Pancreatic_Cancer_Pathway_ WP2377_71228	0.009160941	4	200
Hs_Platelet_homeostasis_WP1885_77031	0.005834554	2	31
Hs_Arachidonic_acid_metabolism_ WP2650_76814	0.000182	3	27
Hs_Netrin-1_signaling_WP1868_76847	0.000635	3	42
Hs_Integration_of_energy_metabolism_ WP1831_77011	0.004354545	3	81
Hs_GPCR_downstream_signaling_ WP1824_76910	0.005903522	6	406

Nine DCC samples were analyzed, and downregulated pathways were presented.

).

Among downregulated genes, we focused on those associated with Wnt signaling, since they were reportedly associated with the state of oocytes and cumulus cells [20,21]. When relative values of gene expressions were estimated, the expressions of five genes were significantly reduced in DCCs. Two cytokines (WNT9B and WNT10B) were significantly downregulated in DCCs (4.095 vs. 1.0, and 2.444 vs. 1.0, DCCs and the control, respectively; p = 0.0007 and p = 0.0001). Frizzled-7 (FZD7), one of the Wnt receptors, was expressed in cumulus cells, and its expression was significantly reduced in DCCs (2.052 vs. 1.0; p < 0.0001). AXIN2 binds to various proteins constituting the Wnt signaling pathway, and its expression was also significantly reduced in DCCs (2.150 vs. 1.0; p < 0.0001). Besides, the expression of Cyclin D1 (CCND1), a representative gene among target genes associated with cell proliferation and the cell cycle, was also reduced in DCCs (2.835 vs. 1.0; p = 0.0001) (Figure 3a). As depicted in Figure 3b, we presumed that these reduced gene expressions impair the transduction of the Wnt signaling pathway in the DCCs.

Based on the above results, we considered the possibility that the Wnt signaling pathway might not be acting only in DCCs compared to the control. Therefore, immunostaining was performed using β-catenin in DCCs and control. The cells were immunostained with β-catenin, and the nuclei were stained with DAPI. Of these cells, DCCs and control exhibited remarkable immunoreaction with β-catenin, but a difference in fluorescence patterns between these cells was also noted (Figure 4a).

Although β-catenin was diffusely localized in the cytoplasm of DCCs and control, it also appeared to be localized in intracellular structures, such as the nuclear membrane in DCCs (Figure 4b,c). The number of β-catenin-expressing cells was higher in DCCs than in control (Figure 4d). Our results suggested that excess β-catenin expression was similar in DCCs and control, but subcellular β-catenin localization was distinct between these cells, implying that the progression in CC maturation may result from altered localization β-catenin. Therefore, we assumed that these cells might be classified differently in different situations.

Dysfunction of the Wnt/β-catenin signaling pathway may be a factor, resulting in follicle atresia, lower fertilization rate, or abnormal embryogenesis, compared to other cumulus cells. These results suggest that DCCs are indeed present in atretic follicles and can be distinguished from other cumulus cells in the maturation process by morphological evaluation.

4. Discussion

In this study, we found that DCCs showed the most variations in the activation of the Wnt signaling pathway compared to the control. It has been reported that the importance of the Wnt/β-catenin signaling pathway in human cumulus cell function is related to follicle development and folliculogenesis [20,21,22,23]. DCCs were suggested to have different follicular development from that of control. The Wnt signaling pathway has two activation mechanisms, depending on β-catenin or not. Cytoplasmic β-catenin is kept at low levels through continuous proteasome-mediated degradation, which is controlled by the GSK-3/APC/AXIN complex. When ligands bind to receptors FZD and lipoprotein receptor-related proteins 5 and 6 (LRP5/6), the degradation of β-catenin is inhibited due to dissociation of β-catenin degradation complex, and consequently, β-catenin accumulates in the cytoplasm and nucleus. Nuclear β-catenin interacts with transcription factors, such as lymphoid enhancer-binding factor 1/T cell-specific transcription factor (LEF/TCF) to affect transcription [20,24].

Among the 19 Wnt ligands identified in humans, DCC showed a remarkable decrease in the expression of two ligands, WNT10B and WNT9B. WNT10B reportedly functions in the Wnt/β-catenin signaling pathway [21]. While there are no reports directly associating WNT9B with Wnt/β-catenin signaling, WNT9B and WNT3 or WNT10B and WNT1 are closely located on the same chromosome of 17q21 or 12q13. WNT1, 2, 3, and 3a [25] ligands have strong transformation abilities, leading to the activation of the Wnt/β-catenin signaling pathway, and conceivably, it may be associated with β-catenin activation. Especially, WNT1 strongly transforms C57MC mammary epithelial cells [25] and participates in the activation of the Wnt/β-catenin signaling pathway. Accordingly, WNT9B and WNT10B likely disturb the activation of the Wnt/β-catenin signaling pathway. In immunostaining, localization of β-catenin in DCCs was different from that in the control, and the number of β-catenin-expressing cells in the nucleus was also significantly lower. In addition, in the control, the nuclei are round and regular, while the nuclei of DCC shows irregular form, likely nuclear fragmentation, implying that DCC may be a state of apoptosis during atresia. These results suggested that the Wnt/β-catenin signaling pathway is inactivated in DCCs. Ten FZD members are capable of functioning as Wnt receptors [26]. There is no published report describing a direct relationship between FZD7 and cumulus cells. Conceivably, FZD7 expression is considered to be reduced, because the Wnt/β-catenin signaling pathway is not active. When the Wnt/β-catenin signaling pathway is inactive, AXIN forms a complex with APC and GSK-3β to act on β-catenin degradation [23]. DCCs also showed a decrease in AXIN expression, possibly due to inadequate transcription of AXIN, without the complex formation in the Wnt/β-catenin signaling pathway. Furthermore, due to lower levels of intracellular β-catenin, the target genes, such as CCND1, could not be transcribed, resulting in a decreased expression.

Besides the genes related to Wnt signaling pathway and pluripotency, we found that six pathways were significantly downregulated (Table 1). Among them, G protein-coupled receptors (GPCRs) and Netrin-1 are reportedly associated with the Wnt signaling pathway. GPCRs comprise a very large family of proteins that play an essential role in cell proliferation, differentiation, neurotransmission, physiological processes, development, and apoptosis [27]. FZD family, the receptors for Wnt signaling, include GPCR [28], and the present results suggested that downregulation of FZD7 led to a decrease in GPCR activity. Netrin-1 is the first axon guidance molecule found in vertebrates that has potent chemotaxis functions for axon guidance, cell migration, morphogenesis, and angiogenesis [29]. Netrin-1 has been reported to significantly promote β-catenin, which is vital for Wnt signaling after spinal cord injury [30], and reduced Netrin-1 activity in cumulus cells may negatively affect Wnt signaling activity. Moreover, regarding the other four pathways, the pancreatic cancer pathway is comprehensive, including the cell cycle and apoptosis. It is reportedly involved in the inhibition of pancreatic cancer via apoptosis [31] and may also be involved in cumulus cell apoptosis. Platelet homeostasis pathway reportedly plays an important role in the regulation of cell coagulation and hemostasis, as well as being involved in the regulation of angiogenesis and apoptosis [32]. Cumulus cells may be similarly affected by the decline in the platelet homeostasis pathway. Since arachidonic acid metabolism by cyclooxygenase (COX) is a major pathway for platelet activation [33], a decrease in arachidonic acid metabolic activity may have also reduced platelet homeostasis activity and led to apoptosis. In the integration of the energy metabolism, mitochondria play a central role, and they may be associated with apoptosis regulation [34,35]. Although these four pathways induce apoptosis as a result of their decreased activity, they are not associated with follicular development. Therefore, we did not consider them important and concentrated on the Wnt signaling pathway.

In conclusion, the DCC developed via controlled ovarian stimulation during the oocyte retrieval cycle suggests that dysfunction of the Wnt/β-catenin signaling pathway may be a factor, resulting in follicle atresia, lower fertilization rate, or abnormal embryogenesis compared to other cumulus cells. These results suggest that DCCs are indeed present in atretic follicles and can be distinguished from other cumulus cells in the maturation process by morphological evaluation. Gene expression analysis of cumulus cells is a useful and safe method that evaluates the condition of oocytes. Based on the gene expression patterns of cumulus cells obtained in this study, there is evidence that they largely contribute to the development of future ART.