RNA sequencing-based transcriptome analysis of granulosa cells from follicular fluid: Genes involved in embryo quality during in vitro fertilization and embryo transfer.

The ability to identify oocyte quality remains one of the most significant challenges in assisted reproductive technology. The quality of oocyte is an important predictor of implantation and live birth [ 1 , 2 ]. The presence of a mature and high-quality oocyte plays an essential role in the development of a high-quality embryo [ 3 ]. This means that the selection of high-quality embryos begins at the time of oocyte selection. Embryo quality is a strong indicator for the success rate of in vitro fertilization (IVF) program, as the live birth rate increases when good quality embryos are transferred [ 4 ]. The incidence of only poor-quality embryos is usual in IVF cycles of patients with advanced maternal age and in low responding patients. According to Semondade et al. [ 5 ], the incidence rate of only poor embryos is approximately 10% at the first IVF cycle, and the recurrence rate is 3%. Morphological features and development rates are key indicators for oocyte and embryo selection during IVF. These methods are familiar for both clinician and embryologist and are used as the standard method in embryos selection step. However, the usefulness of their assessment is being questioned because of personal bias from the embryologist [ 6 ]. In recent years, techniques for embryo selection that provide chromosomal analysis to improve clinical pregnancy have recently been developed, such as preimplantation genetic testing (PGT) [ 7 ]. Although PGT is a strong predictor for implantation, it is expensive and requires invasive embryo biopsy, which involves technical expertise [ 8 ]. Therefore, there is a high demand for a non-invasive and easy-to-perform screening tools to improve the selection of the most pregnancy-competent embryo. Good-quality and mature oocyte in IVF are important for fertilization and embryo development [ 9 ]. Oocyte maturation occurred during folliculogenesis through orchestrated cross-talk between the oocyte and granulosa cells (GCs). Therefore, it is advocated that the functions of GCs indirectly reflect oocyte developmental competence. GCs can be easily recovered in large quantities during oocyte collection. Thus, the gene expression analysis in GCs could provide a non-invasive assessment for identifying the most competent oocytes and embryos. However, the transcriptomic analysis of GC used to identify embryo quality is still controversial. While some studies have identified candidate genes expressed in GC that could be expected as biomarkers of oocyte and embryo quality [ 10 , 11 ] and successful clinical outcomes [ 12 – 14 ], others have reported that there are no significant differences in gene expression between embryos that did or did not successfully implant [ 15 , 16 ]. Therefore, the present study aims to investigate and compare the GC transcriptomic obtained from subjects producing good-quality embryos and those producing only poor quality embryos in human IVF. We carried out RNA sequencing (RNA-seq) of GCs isolated from follicular fluid to identify novel gene transcription factors correlated with embryo quality among the embryos with good quality.

There were 141 women undergoing IVF and embryo transfers (ETs), which included 48 women who had good-quality (GQ) ETs, 79 with only poor-quality (PQ) ETs, and 14 with either GQ or PQ that were excluded from this study ( Fig 1 ). Most patients (n = 83) underwent double ETs, 59 samples came from a single ET, and 3 from triple ETs. ETs for 91 samples occurred on day 3 or 4, and for 50 on day 5, no significant difference between the GQ and PQ groups according to the day of ET. Table 1 presented the characteristics of patients between the GQ and the PQ groups. There were no significant differences in age, body mass index, infertility duration, and anti-Mullerian hormone level. However, there were more ETs in the PQ group than in the GQ group. The clinical pregnancy rate was higher in the GQ (52.0%) than in the PQ embryo transfer group (40.5%); however, the difference was statistically insignificant ( p   =  0.09), and the ongoing pregnancy rate of the GQ group was significantly higher than that of the PQ group (45.8% vs. 24%, respectively, p   =  0.02). BMI, body mass index; IVF, in vitro fertilization; AMH, anti-Mullerian hormone; E2, estradiol; FSH, follicle-stimulating hormone; LH, luteinizing hormone; P4, Progesteron 4; hCG, human chorionic gonadotropin Three GCs samples from each group were selected and used for RNA-seq analysis. Six hundred and forty genes were differentially expressed in our cohort. We analyzed DEGs between the two groups. Compared with that in the GQ embryo group, 18 genes were upregulated, whereas 208 genes were downregulated in the PQ embryo group ( p 1) ( Fig 2A ). As shown in Fig 2B , we developed unsupervised hierarchical clustering of DEGs. A detailed list of 226 DEGs that were upregulated and downregulated in each group is presented in S2 Table . (A) A Venn diagram depicting the distribution of significant differentially expressed genes (DEGs) (fold change > 2.5, q <0.05) in our study. Numbers represent the number of different genes showing upregulated or downregulated expression in PQ embryo samples compared with GQ embryo samples with red representing upregulated expression and blue indicating downregulated expression. Numbers in parenthesis represent the number of DEGs in only PQ or GQ group. (B) Heatmap of the DEGs between the GQ and PQ embryo groups. The color key from blue to red represents the relative gene expression level from low to high, respectively. (C) Volcano plot showing DEG analysis between GQ embryo and PQ embryo groups using DESeq2; 640 genes were differentially expressed (210 downregulated (FC < -2 and FDR 2 and FDR < 0.05). The upregulated genes are represented using red dots; downregulated genes are denoted using blue dots; and the black dots indicate genes with no significant changes. To further extend the molecular properties of the 226 DEGs, we performed a GO analysis of the up- or down-regulated groups ( Fig 3 ). The DEGs were divided into three categories: biological processes, molecular functions, and cellular components. In the GO category biological process, DEGs were enriched in the response to stimulus, metabolic process, biological regulation, immune system process, binding and metabolic process, cellular process, single-organism process, cell and cell part, developmental process, cellular component organization or biogenesis, and the reproductive process. The top 20 genes displayed in the different embryo qualities are presented in Tables 2 and 3 . KEGG pathway enrichment analysis showed that the DEGs in our study participated in 13 pathways ( Fig 4 ). Several of these are related to reproduction, such as the chemokine signaling pathway, phagosome, cytokine-cytokine receptor interaction, cell adhesion molecules, and NF-kB signaling pathway. Genes were classified into three GO domains: biological process, cellular component, and molecular function. The left y-axis shows the gene count in each category. The solid columns indicate DEGs. Fold change was calculated using the 2 ΔΔCt method between the good-quality (GQ) and poor-quality (PQ) embryos. * p <0.05. (A) Six genes associated with embryo quality both in our and previous studies. (B) Seven genes randomly selected from the NF-kB pathway. When comparing our DEGs to those previously reported in the literature, Table 4 summarizes the findings of the detailed list of 68 genes that show genes positively associated with embryo quality. The GC gene biomarker studies have varied depending on the endpoints chosen, including ovarian development, follicular development, folliculogenesis, oocyte and embryo quality. GC, granulosa cells; RT-qPCR, reverse transcription- quantitative polymerase chain reaction; RNA seq, ribonucleic acid sequencing To examine the reliability of the RNA-seq data, we selected 13 genes to verify their expression in GC using RT-qPCR ( Fig 4 ). When comparing our DEGs to those reported in previous studies, we found six genes that were related to embryo quality both in our study and previous published literature, including colony stimulating factor receptor 1 (CSFR1) [ 38 ], cathepsin H (CTSH) [ 37 ], cadherin 13 (CDH13) [ 35 ], serpin family A member 1 (SERPINA1) [ 39 ], cytochrome P450 family 27 subfamily A member 1 (CYP27A1) [ 40 ], and integrin subunit beta 2 (ITGB2) [ 16 ]. Among the significantly enriched pathways, we randomly selected seven genes from the NF-kB pathway: four genes that have refuting findings between different studies (interleukin 1 beta (IL1β), tumor necrosis factor (TNF), BCL2 related protein A1 (BCL2A1), C-C motif chemokine ligand 4 (CCL4)), and three genes that did not support an association with pregnancy or embryo quality (MYD88 innate immune signal transduction adaptor(MYD88), mitogen-activated protein kinase kinase kinase 7(MAP3K7), and TGF-beta activated kinase 1 binding protein 1(TAB1)). Thirteen candidate genes, comprising eleven downregulated genes in GCs from only PQ embryo groups, CSF1R, CTSH, IL1β, TNF, BCL2A1, CCL4, SERPINA1, CYP27A1, ITGB2, MYD88, MAP3K7, and two upregulated genes in GC from PQ embryo groups, TAB1, and CDH13, were selected and analyzed using RT-qPCR. The RT-qPCR results were consistent with the RNA-seq data, which means that the RNA-seq results were dependable and that it could be used to perform accurate differential expression analysis of mRNA.

We found that the expression profiles of specific genes in GCs was associated with embryo quality during IVF and suggested the need to develop treatment strategies that can compensate for the poor quality due to the deficiency of specific genes during follicular development. Furthermore, we reported the published data that support or confute granulosa expressed genes as biomarkers for identifying oocyte or embryo quality. This study improves our understanding of reproductive function of GCs, which could be helpful for more targeted studies aiming to improve oocyte and embryo competence in the future.

This prospective cohort study was conducted between January 2019 and February 2021 at the CHA Fertility Center Gangnam. The study was approved by the CHA University Gangnam CHA hospital institutional review board (GCI-19-10, May 15, 2019), Republic of Korea. All women gave written informed consent to provide material for this study. All procedures followed the rules for studies with human-origin materials established by the IRB. Of 141 embryos transfer (ETs), 48 were good quality (GQ) and 79 were poor quality (PQ) ETs; 14 cases of GQ and PQ ETs cases with endometriosis, endometrial pathologies, uterine fibroids or hydrosalpinx were excluded. Couples with severe male factor infertility as defined based on severe oligozoospermia (<5 million sperm/mL) or a history of testicular biopsy were excluded. Patients underwent controlled ovarian stimulation using either the midluteal long gonadotropin-releasing hormone (GnRH) agonist protocol or GnRH antagonist protocol. Gonadotropin doses were individualized according to patients’ age, anti-Mullerian hormone level, antral follicle count, and previous response to stimulation. Cycle monitoring with transvaginal ultrasonography and serum estradiol measurement, was continued until hCG administration. When at least two follicles reached 18 mm in diameter, 250–500 μg of recombinant hCG was administered for final oocyte maturation. Transvaginal ultrasound-guided oocyte-retrieval under conscious sedation was carried out 34–36 hours after hCG administration. The mature oocytes were inseminated using intracytoplasmic sperm insemination (ICSI). One or two embryos were transferred to each patient, and ETs occurred on day 3, 4, or 5. All patients received luteal phase support with progesterone after ET until 8–10 weeks after pregnancy. Follicular fluids were aspirated and pooled for each patient during oocyte-retrieval. Follicle aspirates, which were not clear and were contaminated with endometriosis cysts, were discarded. Granulosa cells recovered from one woman were used as one sample. The 127 samples were divided into two groups according to the patient’s embryo quality; 48 patients with good quality blastocysts having a grade of at least 3BB and 97 patients with only poor quality blastocysts. The follicular fluid was centrifuged at 1000 x g for 20 min to separate erythrocytes, leukocytes and GCs. The cell pellets were washed in RBC lysis buffer (Roche Diagnostics, Basel, Switzerland) for RBC removal and centrifuged at 500g for 10 min and separated enzymatically with instigation at 37°C for 30 min in the enzyme solution [Hanks’ balanced salt solution (HBSS, Gibco, Grand Island, NY) containing 0.5 mg/ml collagenase (type IV; Gibco) and 0.25mg/ml dispase II (neutral protease, grade III, Roche)]. The suspension of GCs was washed and filtered through a 40μm mesh (BD, Franklin Lakes, NJ). GCs using TRIzol Reagent (Life Technologies, San Diego, CA, USA) were stored -75°C until analysis. Embryo grading was evaluated by two embryologists and was assessed by another senior embryologist with over at least ten years of work experience before embryo transfer. For standardization, all embryologists at our center were trained in embryo grading by the same laboratory director. Embryo quality assessment was carried out according to previously described protocol [ 17 ]. Briefly, good quality embryos were identified using the following characteristics. Cleavage embryos were classified as GQ according to Cummins criteria [ 18 , 19 ]. Blastocyst quality was assessed regarding to the degree of blastocoel expansion, inner cell mass (ICM) and trophectoderm (TE) morphology on day 5 or 6 [ 20 ]. Total RNA was isolated from cells using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. RNA integrity was checked by using an Agilent 2100 BioAnalyzer (Agilent, CA, USA) with an RNA integrity number value greater than 6. Only qualified samples were used for RNA library constitution. The libraries were arranged for 151 bp paired-end sequencing by the TruSeq stranded mRNA sample preparation kit (Illumina, CA, USA). After the sequential process of end repair, A-tailing, and adapter ligation, cDNA libraries were amplified by Polymerase Chain Reaction (PCR). They were quantified using the KAPA library quantification kit (Kapa Biosystems, MA, USA) following the manufacturer’s quantification protocol. After cluster amplification of denatured templates, sequencing was carried out as paired-end (2×151bp) using IlluminaNovaSeq6000 (Illumina, CA, USA). Based on the estimated read counts in the previous step, DEGs were screened using the R package TCCv.1.26.0 [ 21 ]. The TCC package uses well-set normalization strategies to compare the tagcount data. Normalization factors were estimated by the iterative DESeq2 [ 22 ] /edgeR method [ 23 ]. The Q-value was assessed based on the p-value using the p.adjust function of the R package by default parameter programs. The DEGs were evaluated with a Q- value threshold of < 0.05 for fixing errors due to multiple testing. The GO database indicates a set of hierarchically controlled vocabulary classified into three categories: biological process, cellular component, and molecular function. For functional characterization of the DEGs, GO-based trend test was assessed by the R package called GOseq [ 24 ] through the Wallenius non-central hypergeometric distribution. Selected genes with P-values < 0.05 following the test were regarded as statistically significant. Approximately 1 μg RNA from each sample was reverse transcribed to cDNA using a SensiFAST™ cDNA Synthesis Kit (Meridian Bioscience, Tennessee, USA) following the manufacturer’s instructions. Thirteen genes were selected for validation from the list of DEGs. Primers were designed using the Primer3 primer-design program (version 0.4.0; http://bioinfo.ut.ee/primer3-0.4.0/ ) ( S1 Table ). The qPCR reactions were performed on the CFX Connect Real-Time PCR detection System (Bio-Rad Laboratories, Inc., USA) using iQ™ SYBR® Green Supermix (Bio-Rad Laboratories, Inc., USA). Relative quantification analyses were carried out using the comparative CT method, and relative gene expression levels were calculated using the 2 −ΔΔCT method. The primary outcomes of the study were clinical pregnancy outcomes, including clinical pregnancy rate, and ongoing pregnancy rate. Clinical pregnancy was defined as the presence of at least one gestational sac with a fetal heartbeat on ultrasonography. Ongoing pregnancy was defined when a positive heartbeat was at 12 weeks or more of gestation on ultrasonography. All statistical analyses were conducted using SPSS (version 25.0; Chicago, IL, USA) software. Categorical variables were performed using the chi‐square and Fisher’s exact tests. Continuous variables were assessed using Student’s t test. A probability ( p ) value <0.05 was considered to determine statistical significance. A p value below 0.05 was considered statistically significant.

(DOCX) Click here for additional data file. (XLSX) Click here for additional data file.