DNA Methylation Changes in Endometrium and Correlation with Gene Expression During Ovarian Hyper-stimulation

Persistent supraphysiological serum estradiol (E2) levels during controlled ovarian hyper-stimulation (COH) have a detrimental effect on endometrial receptivity. In this study, we explored RNA expression and DNA methylation proles from patients’ endometrium. The patients were divided into two groups: the COH cycle (n=3, hCG+7) group and normal cycle group (n=3, LH+5). Quantitative RT-PCR was used to validate the expression of selected differentially expressed genes (DEGs). Comparing natural and stimulated endometrium transcriptome proles revealed 640 DEGs, with a > 2-fold change (FC) and p < 0.01. The DEGs were reported to be involved in endometrial receptivity and epithelial-mesenchymal transition (EMT). The expression of IGFBP-1, MMP9, FGF9, LIF, WNT4, HAND2, and immune system-related genes were signicantly up-regulated. By clustering and KEGG pathway analysis, molecules and pathways associated with endometrial receptivity (PI3K-Akt signaling pathway) were identied. DNA methylation was partially correlated to gene expression. In conclusion, RNA-seq COH affected endometrial receptivity and EMT/MET process by accelerated the decidualization process and broken the balance of estrogen and progesterone receptors expression. However, this was not associated with changes in DNA methylation.


Introduction
Controlled ovarian hyperstimulation (COH) is a technique widely used in vitro fertilization (IVF). Highquality embryos and a receptive endometrium are signi cant during COH. In recent years, several COH protocols have been adopted to optimize the quantity and quality of embryo transferred and improve on subsequent pregnancy success rate. However, even with a high-quality embryo, implantation failure can occur due to insu ciency of the endometrium de ciency and this is a major cause of pregnancy failure.
Previous human studies [1,2] have shown that the supraphysiological hormonal levels, especially serum estradiol (E2) impair endometrial receptivity or displace the window of implantation (WOI) affecting embryo adhesion [3]. In vitro studies have also con rmed that increasing E2 levels are deleterious to embryonic implantation because they directly affect the embryo [4] and endometrium [5].
Endometrial alterations in clinical practice are characterized by the endometrial thickness and morphological alterations which require further histological [6] and biochemical investigations [7].
Recently, researchers have focused on omics analysis of changes in the endometrium during COH. In stimulated cycles, progesterone receptor (PGR) and related pathways are reduced during WOI, while estrogen receptor (ER) and related pathway are increased accordingly [8][9][10][11]. In natural as well as stimulated cycles, the majority of the modulated genes during endometrial receptivity are reported to be up-regulated. Besides, the genes involved in the implantation process are absent in COH cycles [12][13][14].
Studies have reported that the endometrial genomic pro le in the presence of a GnRH antagonist protocol is more similar to the natural cycle compared to that treated with a GnRH agonist protocol [15][16][17].
Therefore, there is need to optimize COH protocols to improve endometrial receptively. DNA methylation (DNAm), as one of the most common forms of epigenetic modi cation. Methylation pro les in human endometrium also change across the menstrual cycle with thousands of genes differentially methylated between cycle stages [18]. A comparison of the changes in transcriptomes and corresponding DNAm from the same samples showed an association of DNA methylation and gene expression [18]. Genes involved in cell adhesion, immune response and developmental processes are the most affected by aberrant methylation changes during epithelial-mesenchymal transition (EMT)/mesenchymal-epithelial transition (MET) regulation processes [19,20]. Previous studies have shown that DNA hyper-methylation has a detrimental effect on endometrial receptivity in the presence of high progesterone (P4) levels during COH [21]. E2 and P4 induced decidualization of human endometrial stromal cell in vitro has been shown to exhibit DNAm changes [18,22]. However, the status of DNAm pro les in supraphysiological E2 levels during COH and their associations with changes in gene expression remain unknown.
In the present study, we explored the effect of supraphysiological levels of serum E2 on DNAm pro les during COH and correlate with gene expression pro les of endometrial receptivity.

Material And Methods
Patient characteristics and tissue collection. Endometrial samples (3 COH and 3 normal cycle specimens) were obtained from 6 healthy infertile patients. Inclusion criteria included age range between 25 -35 years, BMI between 18 -25kg/m 2 and a regular menstruation cycle of 21-35 days. In the NOR group, tissues were obtained from the patients at their mid secretory phase (LH+5) of the natural cycle and the patients had not been on hormonal treatment for at least three months before tissue sampling. For the COH group, gonadotrophin-releasing hormone agonist (GnRH-a) and Gn were used for ovarian stimulation. Oocyte retrieval was performed 36-37h after HCG administration by the transvaginal ultrasound approach. Insemination was carried out 4 to 5 hours after ovum pick-up (OPU). About 16 to 18 hours after insemination, the presence of pro-nuclei was detected and freezing of all embryos performed. Patients were diagnosed with ovarian hyperstimulation syndrome (OHSS) and after 5 days of OPU (hCG+7), their endometrium was collected.
Endometrial biopsies were collected with Pipelle catheters under sterile conditions from the uterine fundus. A portion of the collected endometrial tissue was snap-frozen in liquid nitrogen and stored at -80°C in RNAlater (Ambion Inc., Austin, TX, USA) until further processing. Tissue collection was approved by the Ethics Review Committee of the University of Zhengzhou, and written informed consent was obtained from all participants.

RNA extraction and expression analysis.
Total RNA was extracted from endometrial tissues using the TRIzol reagent (Invitrogen, Carlsbad, CA) according to the manufacturer's instructions. The RNA quality was determined using a NanoDrop 2000 Spectrophotometer (ThermoFisher, USA). Total RNA extracted was used for cDNA library construction as described in Supplemental Figure 1 and sequencing performed using Illumina HiSeq TM 2500 (Gene Denovo Biotechnology Co., Guangzhou, China). All the high-quality reads were mapped using Bowtie2 [23] to ribosomal RNA (rRNA). All rRNA reads were mapped to the human genome (GRCh38/hg38) using TopHat2 (version 2.0.3.12) [24]. Cu inks [25] and Cuffmerge were used to reconstruct and merge the transcripts.
Principal component analysis (PCA) was performed with R package gmodels (http://www.r-project.org/). Differentially expressed genes (DEGs) were identi ed using edgeR with a fold change (FC) >2 and a false discovery rate (FDR) <0.05. The DEGs were then subjected to enrichment analysis including Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis.
Quantitative real-time PCR was used to validate the RNA sequencing (RNA-seq) results. Synthesis of cDNA from extracted total RNA (2ug) was performed with a Fast Quant RT kit (TIANGEN, China) according to the manufacturer's instructions. Real-time PCR was performed using the SYBR Green master mix (Vazyme, China) on a Step One Plus Real-Time PCR System (Applied Biosystems, USA). All quantitative PCR reactions were run in triplicate. Primer sequences are presented in Supplemental Table  1.

DNA extraction and methylation analysis.
Genomic DNA was isolated from endometrial tissue using Qiagen DNeasy Kits and quanti ed using Pico-Green assay. DNA concentration and integrity were assessed using a Nano Drop spectrophotometer and Agarose Gel Electrophoresis, respectively. DNA bisul te sequencing libraries were prepared as described in Supplemental Figure 1. After ltering, the clean reads were mapped to the species reference genome using BSMAP software (version: 2.90). A custom Perl script was used to call methylated cytosine's which were tested with the correction algorithm as described in previous studies [26].
To identify DMRs between two samples, the minimum read coverage to call a methylation status for a base was set to 4. DMRs within a 200bp window for each sequence context was determined according to different criteria: 1) For CG and CHG, the numbers in each window ≥ 5, the absolute value of the difference in methylation ratio ≥ 0.25, and q ≤ 0.05; 2) For CHH, the numbers in a window ≥ 15, the absolute value of the difference in methylation ratio ≥ 0.15, and q ≤ 0.05; 3) For all C, the numbers in a window ≥ 20, the absolute value of the difference in methylation ratio ≥ 0.2, and q ≤ 0.05.

Statistical analysis
The data for patients in the COH and NOR groups were expressed as the mean ± standard deviation or number (percentage). The Mann-Whitney U test was used to study the continuous variable and non-normally distributed data were tested by Wilcoxon's rank-sum test. For the categorical variable data, we performed chi-squared test analysis or 2-way ANOVA was performed using SPSS software (version 21.0). *P<0.05 was considered to be statistically signi cant.

Baseline parameters
The baseline characteristics of the six patients are listed in Supplemental Table 2. Age, duration of fertility, body mass index, type of infertility and antral follicle count were comparable between the two groups. No statistically signi cant difference was reported in the basal hormone levels.  Figure 1B shows volcano plots providing an overview of the differential status of gene expression.
GO and KEGG pathway analysis was performed on the DEGs. Enriched Go terms are showed in Figure 2A B and arranged and displayed according to biological processes, cellular components, and molecular functions. For the cellular component annotation classi cation, both up-and down-regulated genes were mainly localized in the membranes, endomembrane system, vesicle, nucleus, and extracellular region. Concerning the molecular function, both up-and down-regulated genes were mainly associated with binding, catalytic activity, and signal transduction. The enriched biological processes for the up-and down-regulated genes were mainly involved in biological regulation, metabolic processes, response to stimulus and multicellular organismal process. KEGG pathway analysis showed that the top 20 KEGG pathways were signi cantly enriched and are shown in Figure 2C. Among these pathways, the immune system (complement and coagulation cascades) and signal transduction (PI3K-Akt signaling pathway) were found to be the most enriched. Besides, the Ras signaling pathway, MAPK signaling pathway, and WNT signaling pathway were also enriched.
To further investigate the DEGs relationships, a co-expression network of 640 DEGs (p-value < 0.01) based on protein-protein interaction (PPI) is shown in Figure 3A. Four key genes, including C3, CXCL8, CXCL12, and PENK were considered to have important regulation and control ability. A total of 221 DEGs related to EMT, 83 DEGs involved in endometrial receptivity and 57 DEGs associated with both functions were screened ( Figure 3B, Supplementary Data 2). The six up-regulated transcripts with the highest FC in EMT were MMP1, SOHLH2, ELOVL3, LBP, IL1B, and IGFBP1. The six down-regulated transcripts identi ed as receptivity-associated genes were NLRP5, SLC7A4, SEPRINB9, CD36, CREB3L1, and ATP6V0A4.

DNAm pro ling analysis
Genome-wide DNAm pro les were studied by reduced-representation bisul te sequencing (RRBS-Seq). The mean genome-wide methylation levels of the six samples in C, CG, CHG, and CHH are shown in Supplementary Table 3. CG sites (mCG) modi ed by methylation were found to be stable and fully methylated. The methylated CHH sites (mCHH) were unstable and were in a hemi-methylated or hypomethylated state.
The DMRs statistics between NOR and COH groups are shown in Figure 4A. The violin plots showing the distribution of methylation in each locus on 10 genes revealed signi cant differences in methylation ( Figure 4B). The region-level analysis of all CG revealed 38 signi cant DMRs, of which 20 exhibited increased (associated with 10 genes) and 18 decreased (associated with 10 genes) methylation (Supplementary Table 4). The most signi cant DMRs included CG in the upstream region of RAP1B and the gene-body of NPAS3. RAP1B was also one of the genes highlighted in the site-level analysis (Supplementary Data 3).
GO and KEGG pathway analyses were performed to characterize the genes annotated to DMRs. In sitelevel analyses, no GO terms were signi cantly represented (Supplementary Data 4). This is typically mirrored by region-level analyses (Supplementary Data 5). Pathway enrichment analyses for the same gene lists showed enrichment in 12 pathways in site-level analysis, including insulin signaling, Rap1 signaling, renal cell carcinoma, neurotrophic signaling, N-Glycan biosynthesis, measles, Jak-STAT signaling, transcriptional misregulation in cancers, endocytosis, cytokine-cytokine receptor interaction, HTLV-I infection, PI3K-Akt signaling pathway. No enrichment was seen in region-level analysis; however, genes showing a correlation between methylation and gene expression were enriched for the Wnt signaling pathway.

Correlation analysis
Screening for gene function and expression and selected TSPAN4, NPAS3, RAP1B, and TGFB3 for qPCR was performed to evaluate the correlation between DNAm and transcriptome (Supplementary Figure 2). In this study, there was a limited correlation between methylation and transcriptome during COH.

Discussion
The quality of endometrial receptivity during the WOI or the change of endometrial development time is highly implicated in implantation rates. Over the last decades, endometrial receptivity and EMT/MET processes have been widely studied. However, the detailed molecular mechanism during COH cycles remains unclear. Supraphysiological hormonal treatment in IVF cycles has shown a considerable impact on the molecular alterations in the endometrium. The main changes of molecules are focused on the endometrial glandular and stromal cells secreting and differentiating [28]. Elevated serum E2 levels in COH cycles are associated with increased glandular-stromal dyssynchrony and defective induction of PGR and a considerable decline in the number of cytosolic PGRs [3]. Besides, an increase in E2 and P4 causes endometrium transfer from the proliferative phase to the secretory phase and this leads to the WOI opening. Stromal cell decidualization reduces endometrial receptivity. High E2 levels limit the development of pectin but stimulate hyperplasia of endometrial collagen bers, which can cause receptivity decline and implantation failure. However, with omics detection technology, it is possible to accurately and reliably assess the endometrium in clinical routine practice. Therefore, this study investigated the RNA expression and DNAm pro le in patients who underwent the COH cycle in comparison with the natural cycle.
In this study, ER1 was reported to be signi cantly decreased and ER2 increased in the COH group compared with the natural group. E2 mediated epithelial proliferation occurred through the nuclear-and membrane-initiated steroid signaling. In the nuclear pathway, the down-regulation of ER1 was predominant. In the membrane pathway, E2 exerted its actions through a subpopulation of ER at the plasma membrane (mER) and GPER. Upon activation of these receptors, various signaling pathways were rapidly activated which ultimately in uenced downstream transcription factors.

PI3K/AKT signaling pathway
The PI3K/AKT pathway was signi cantly enriched in various signal transduction pathways. The binding of growth factors (GF) to tyrosine kinase (RTK) receptors and G protein-coupled receptors (GPCR) stimulated PI3K class IA and IB isoforms, respectively.PI3K catalyzed the production of phosphatidylinositol-3,4,5-triphosphate (PIP3) in the cell membrane. PIP3 served as a second messenger which activated AKT. AKT upregulated key cellular processes by phosphorylating substrates involved in protein synthesis, metabolism, and cell cycle.

RAS/MAPK signaling pathway
E2 not only induces epithelial proliferation in the endometrium but also induces LIF in the glandular epithelium [33]. LIF action on the luminal epithelium regulates WNT/ß-catenin signaling, IGF1 signaling and FGF signaling [33]. LIF was signi cantly up-regulated in our results. FGF family are induced by E2 in the endometrial stroma [34]. As a potential paracrine mediator, FGF induces the RAS/MAPK and AKT pathway. In this study, part of the FGF family, particularly FGF9 was found to be up-regulated. Most of the E2-related pathways were overexpressed but the level of ER remained low. These ndings suggested that proliferation-associated pathways were induced by E2. However, these ndings should be veri ed using a larger sample size.

WNT/ß-catenin signaling pathway
The molecular protagonists in the decidualization route are P4 and cAMP [29]. The PGR-regulated genes and pathways include HOXA10, HAND2, IGFBP-1, and BMP2; ERK/MAPK, WNT/ß-catenin, and IHH-COUPTFII pathways [35]. HOXA10 is a PGR target in the endometrium and has been reported to in uence WNT4 expression around the implantation site in mice [36]. However, HOXA10 was not found in the DEGs, but the HOA family, HOXA5/6/7, HOXD8/9/10/11 were signi cantly up-regulated. Inhibitor DKK1 was found to be upregulated but did not block the WNT canonical pathway and this contributed to endometrial proliferation but not decidualization. Therefore, it is suggested that the decidualization process was partly impeded by the PGR-regulated pathway. The imbalance between ER and PGR related genes and pathways might impair endometrial receptivity.

Complement and coagulation cascades
Previous studies have con rmed that PRL, IGFBP-1, MMP9, and IL-6 are markers of decidualization [37]. Increased MMP secretion promotes endometrial stromal cells to support the decidualization process [38,39]. In this study, both MMP9 and MMP2 were up-regulated suggesting promotion of the MET process. Decidualization is an example of the MET process and was found to be correlated with altered immune modulation and the complement cascade pathway [10]. Excessive E2 levels may be one of the reasons for maternal immune rejection due to the activation of the complement system. Furthermore, we established a co-expression network using the DEGs based on their signi cance and correlation. Both CXCL12 and its receptor CXCR4 were found to be upregulated and this con rmed that the proin ammatory response existed despite none embryo implantation. A large number of in ammatory responses inhibited endometrial decidualization and embryo implantation.
Changes in DNAm during COH revealed 59 signi cant DMRs, out of which 23 exhibited increased and 36 decreased methylation between the two groups. In this study, no signi cant correlation between DNAm and transcriptome was reported. However, the methylation ndings con rmed the dominance of the E2regulated pathway caused by COH. Hence, it changed the endometrial receptivity and EMT process which mainly acted by P4.

Conclusion
In conclusion, RNA-seq COH had a deleterious effect on endometrial receptivity and EMT/MET process by advancing the decidualization process and affecting the balance between ER and PGR. However, these effects were not associated with changes in DNA methylation. This study has the following limitations; small sample size, single-center study, and retrospective design. Therefore, future prospective trials or meta-analyses are needed to con rm these ndings.

Clinical Perspectives
Supraphysiological serum E2 levels during COH is associated with embryo implantation and have an adversed effects on endometrial receptivity. However, the status of DNAm pro les in supraphysiological E2 levels during COH and their associations with changes in gene expression remain unknown.
We con rmed that COH had a deleterious effect on endometrial receptivity and EMT/MET process, but these effects were not associated with changes in DNA methylation.
Our results showing that the gene and DNAm change of endometrial during COH and to provide a theoretical basis of improving the endometrial receptively during IVF/ICSI clinical practices.

Declarations
Ethics approval and consent to participate This study design was reviewed and approved by the Research Ethics Committee of the First A liated Hospital of Zhengzhou University.

Consent for publication
Not applicable.

Availability of data and materials
All data generated or analysed during this study are included in this published article and its supplementary information les. Figure 1 A)Transcriptome analysis based on RNA-seq and PCA showed the clustering of the samples in each group B) volcano plots providing an overview of the differential status of gene expression.