Construction of Competing Endogenous RNA Networks Incorporating Transcription Factors to Reveal Differences in Granulosa Cells from Patients with Endometriosis

Abstract

Purpose:

This study was designed to reveal the molecular differences between granulosa cells (GCs) from patients with endometriosis and normal controls.

Methods:

RNA sequencing was performed on GCs from patients with EM-related infertility (n = 3) and controls (n = 3). Differentially expressed long noncoding RNAs [differentially expressed lncRNAs (DELs), |log2 FC|>4, false discovery rate (FDR) <0.05] and genes [differentially expressed genes (DEGs), |log2 FC|>1.4, FDR <0.05] in patients with EM-related infertility and controls were screened. Protein-protein interaction (PPI) networks of the DEGs were constructed. Then, mRNA-miRNA-lncRNA pairs based on DEGs and DELs were constructed by comprehensive bioinformatics analyses. In addition, overlapping genes identified from both the PPI and the mRNA-miRNA-lncRNA pairs were selected. Finally, a competing endogenous RNA (ceRNA) network incorporating transcription factors (TFs) was constructed.

Results:

A total of 25,806 lncRNAs and 19,684 mRNAs were detected, and 7 DELs and 46 DEGs were identified. Five hub genes from the PPI network were also identified. A single overlapping gene, NR4A2, from both the PPI network and mRNA-miRNA-lncRNA pairs was identified. Finally, a ceRNA network incorporating TFs, including one mRNA (NR4A2), one miRNA (hsa-miR-217), three lncRNAs (XIST, MCM3AP-AS1, and C17orf51), and five TFs (SRF, POLR2A, NRF1, MNT, and TCF7L2), was successfully constructed.

Conclusions:

The proposed ceRNA network and the prediction of TFs in GCs from EM-related infertility revealed differences in GCs from patients with EM. Importantly, the novel TFs, lncRNAs, miRNAs, and mRNAs involved in the ceRNA network might provide new insights into the underlying molecular mechanisms of EM-related infertility.

Introduction

Endometriosis (EM) is a benign gynecological disease caused by the appearance of endometrial tissue (glands and stroma) anywhere outside the uterine cavity. About 50% of patients with EM have infertility (Eisenberg et al., 2018). The mechanisms underlying EM-related infertility have not been fully elucidated (Vercellini et al., 2014). EM is characterized by “omni-directional” interference in the reproductive process, including endometrial receptivity deficiency; influence of pelvic inflammatory microenvironment on the quality of sperm, oocytes, and embryos; abnormal follicular development; epigenetic changes; influence of EM after surgery on ovarian function and pelvic adhesion; and so on (Pallacks et al., 2017; Lessey et al., 2018). Among these, the decrease in oocyte quality caused by the change in the microenvironment in the follicular fluid of patients with EM is one of the important factors (Dumesic et al., 2015). Compared with tubal factor infertility (TFI), EM-related infertility shows slower growth of follicles, decreased dominance of follicles, low oocyte retrieval rate, relatively increased number of empty follicles and atresia follicles, poor oocyte quality, and reduced fertilization rate.

Granulosa cells (GCs), as the main cells for follicular development, can synthesize and secrete estrogen and progesterone, promote follicular development, and provide favorable conditions for follicle maturation. At the same time, they can provide nutrients and energy for the development, maturation, and normal functioning of oocytes. The proliferation, differentiation, and apoptosis of GCs are also key factors that affect follicular development and oocytes (Dumesic et al., 2015). Therefore, its normal function is essential to the growth, development, and maturation of oocytes and follicles. Follicle maturation disorders and poor oocyte quality may be closely related to GCs in the follicle in patients with EM-related infertility. However, the molecular biological mechanism of GCs in EM-related infertility has not been reported.

Therefore, this study aimed to explore the possible causes of EM-related infertility by studying ovarian GCs in patients with EM-related infertility and controls. GCs are the main somatic components in the follicle, and they form a functional whole through a wide gap connection with the oocytes. Therefore, changes in the function of GCs play a decisive role in the growth and development of follicles and the maturation of oocytes.

Previous studies have shown that miRNAs can reduce mRNA stability or inhibit translation by binding to miRNA recognition element (MRE) (Pritchard et al., 2012). A single miRNA can modulate multiple targets containing miRNA-specific MRE, as well as a single RNA containing multiple MREs regulated by multiple miRNAs (Allegra et al., 2020). At the same time, lncRNAs have been reported to act as miRNA sponges that compete with mRNAs to attract miRNAs for interactions and influence the expression of mRNAs. This phenomenon, in which certain genes targeted by common miRNAs “compete” for these miRNAs, thereby regulating each other by making others free from miRNA regulation, is called competing endogenous RNA network (ceRNA network). However, the biological characteristics and clinical relevance of pseudogenes that function as ceRNAs in EM-related infertility remain unclear, especially in GCs from EM-related infertility. In addition, the eukaryotic transcription initiation process is very complicated and often requires the assistance of multiple protein factors. Transcription factors (TFs) and RNA polymerase II form a transcription initiation complex and participate in transcription initiation together.

In this study, total RNA of ovarian GCs was extracted from patients with EM-related infertility and controls. Next, a ceRNA network incorporating TFs was constructed to reveal differences in GCs from patients with EM by RNA sequencing and comprehensive bioinformatics methods. Finally, the novel TFs, lncRNAs, miRNAs, and mRNAs involved in the ceRNA network incorporating TFs in this study might provide new insight into the underlying molecular mechanisms of EM-related infertility.

Methods

Patient history

The study was approved by the Ethics Committee of the First Affiliated Hospital of Xiamen University (KY2017-054). Also, written informed consent was obtained from patients who met the eligibility criteria before oocyte retrievals. The study enrolled six infertile patients (three infertile patients with ovarian EM and three normo-ovulatory patients with TFI), who underwent ovarian stimulation for IVF (in vitro fertilization) at the ART Center (the First Affiliated Hospital of Xiamen University, Xiamen) between November 2018 and December 2018. The patients were diagnosed with EM either by laparoscopy or pathological examination. The stages of patients with ovarian EM were categorized according to the revised American Society of Reproductive Medicine classification (1997) (stage III, n = 3). The inclusion criteria for the ovarian EM and control groups were as follows: (1) a history of infertility more than 1 year; (2) age 20-35 years; (3) normal liver and kidney function, without gynecology and other systemic disease; and (4) basal serum FSH levels <10 IU/L, before controlled ovarian hyperstimulation. The exclusion criteria were as follows: (1) patients combined with ovulation dysfunction, such as polycystic ovary syndrome, POI, and hyperprolactinemia; (2) patients also having serious diseases such as cardiovascular system, liver, kidney, and hematopoietic system; (3) autoimmune disease; (4) uterine fibroids, endometritis, nonvegetative ovarian cysts, ovarian malignancies, and internal genital tuberculosis; and (5) smoking, alcoholism, or drug addiction. No patients with ovarian EM received any hormones within 6 months.

Ovarian stimulation protocols

All patients with EM and controls underwent an improved long protocol; they were injected with long-acting GnRH agonist (Duffelin, Ipson Biotechnology, France) on the second to the fourth day of the menstrual cycle. The dose of GnRH agonist might differ depending on the individual differences. After 28-35 days, the hormone level and ultrasonographic findings were measured, and Gn was started after reaching the downregulation standard. The successful downregulation standard was as follows: serum FSH <5 U/L, LH <5 U/L, serum E2 < 183.5 pM, and progesterone (P) <3.18 nM. The ultrasonography showed endometrial thickness <5 mm and follicle diameter <10 mm in both ovaries. In the process of controlled ovarian stimulation, the starting dose (150-300 U) of Gn (rhFSH, Merck Serono, Germany) was determined according to the ovarian responsiveness, Body Mass Index, and previous ovulation induction of the patients. When the follicle diameter ≥18 mm dominated the follicle by 50-60%, follicular maturation was triggered through an injection of 5000-10,000 IU of hCG (Lizhu Pharmaceutical, China), and the oocytes were retrieved after 38 h.

Collection and purification of GCs

GCs were collected and purified as previously described (Jin et al., 2018). Briefly, GCs were collected from the patient follicular fluid on the oocyte retrieval day. Only clear follicular fluid without blood or flushing solution was collected. Each follicular fluid sample in a 15-mL disposable sterile tube was immediately centrifuged at 1000 rpm for 10 min. After removing the supernatant, the cellular precipitate was mixed with PBS at a ratio of 1:1, slowly added to the upper layer of 45% Percoll (GE Healthcare, Sweden) separation solution, and centrifuged at 1000 rpm for 20 min. Then, the middle ring-like layer (GCs) was collected. GCs from each participant were stored at −80°C for RNA extraction.

RNA extraction, strand-specific library construction, and sequencing

After total RNA was extracted using a TRIzol Reagent Kit (Invitrogen, Carlsbad, CA) following the manufacturer's protocol, ribosome RNAs (rRNAs) were removed to retain mRNAs and ncRNAs. The enriched mRNAs and ncRNAs were fragmented into short pieces using fragmentation buffer and reverse transcribed into cDNA with random primers. Second-strand cDNA was synthesized by DNA polymerase I, RNase H, dNTP (dUTP instead of dTTP), and buffer. Next, the cDNA fragments were purified using a QiaQuick PCR Extraction Kit (Qiagen, Venlo, The Netherlands), end-repaired, mixed with poly(A), and ligated to Illumina sequencing adapters. Then, uracil-N-glycosylase was used to digest the second-strand cDNA. The digested products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced using Illumina HiSeq 4000 (or other platforms) by Gene Denovo Biotechnology Co. (Guangzhou, China). Finally, the raw reads were obtained.

Data preprocessing

First, raw reads obtained from the sequencing machines were further filtered using fastp (Chen et al., 2018b) (version 0.18.0) to obtain high-quality clean reads. After data filtering, the composition and mass distribution of bases were analyzed to visualize the quality of the data. The more balanced the base composition, the higher the quality, and the more accurate the subsequent analysis. Then, the short reads alignment tool Bowtie2 (Langmead and Salzberg, 2012) (version 2.2.8) was used for mapping reads to the rRNA database. The rRNA mapped reads were then removed. The remaining reads were further used in the assembly and analysis of the transcriptome. In addition, an index of the reference genome was built, and paired-end clean reads were mapped to the reference genome using HISAT 2 (Kim et al., 2015) (version 2.1.0) with “-rna-strandness RF” and other parameters were set as a default. Based on the results of all reads (Total_Mapped reads) that could be mapped to the genome, the distribution of reads in the reference genome was calculated. Finally, genetic analysis (clean reads), including gene coverage, randomness analysis, and sequencing saturation analysis, was performed to ensure the accuracy of subsequent analysis.

lncRNA analysis

The transcripts were reconstructed with software Stringtie (Pertea et al., 2016) (version 1.3.4) and HISAT2, which allowed biologists to identify new genes and new splice variants of known ones. The novel lncRNAs were predicted using CNCI (Sun et al., 2013) (version 2) and CPC (Kong et al., 2007) (version 0.9-r2), and their types were analyzed. In addition, lncRNA expression (FPKM) was obtained after correcting sequencing depth and transcript length. The FPKM values of the lncRNAs were used for subsequent analysis. Then, principal component analysis (PCA) (ggplot2 R package) and Pearson correlation coefficient analysis between samples were used to understand the repeatability between samples based on the expression level of lncRNA in each sample, which helped exclude outliers. Finally, the normalized data of lncRNAs were used for differential analysis using the DESeq2 (Love et al., 2014) package, and the differentially expressed lncRNAs (DELs) were identified by fold-change (FC) screening at a threshold of 16.0-fold or greater and false discovery rate (FDR) <0.05. DELs were used for subsequent analysis.

mRNA analysis

Based on the comparison results of HISAT2, the transcripts were reconstructed and the expression levels of all genes in each sample were calculated. The FPKM value of the gene was obtained after correcting sequencing depth and transcript length to ensure the accuracy of subsequent analysis. The FPKM value of the gene was used for subsequent analysis. Differentially expressed genes (DEGs) with the parameter of FDR <0.05 and absolute log2-FC ≥1.4 between two different groups (controls vs. EM-related infertility) were screened out. In this study, DEGs were used for subsequent analysis.

Protein-protein interaction network analysis and screening of hub genes

The Search Tool for the Retrieval of Interacting Genes (STRING, Version: 11.0) (Szklarczyk et al., 2017) can be used to provide information regarding predicted and experimental interactions of proteins. The prediction method involved co-occurrence, neighborhood, coexpression experiments, gene fusion, text mining, and databases. The DEGs were mapped into protein-protein interactions (PPIs), and a combined score of >0.4 was set as a threshold value. PPI networks were constructed and visualized with Cytoscape software (Version: 3.7.2). In addition, the nodes with higher degrees of interaction from the PPI network were considered as hub nodes. The cytoHubba (Version: 0.1) is a tool for screening hub genes in the Cytoscape software. A recent study suggested that the newly proposed method MCC was better than the other 11 previously reported methods in cytoHubba plug-in. The MCC method was used to screen hub genes.

Construction of mRNA-miRNA-lncRNA pairs

The lncRNA-miRNA pairs and miRNA-mRNA pairs were constructed based on miRcode (Version 11), miRDB (Version 7.0), miRTarBase, and TargetScan (Version 7.2) databases using the ggalluvial R package (Version: 1.0.0) to construct mRNA-miRNA-lncRNA pairs (Rosvall and Bergstrom, 2010). MiRcode provides “whole transcriptome” human microRNA target predictions based on the comprehensive GENCODE gene annotation, including >10,000 long noncoding RNA genes. Coding genes were also covered, including atypical regions such as 5′-UTRs and CDS. miRDB is an online database for miRNA target prediction and functional annotations. All the targets in miRDB were predicted using a bioinformatics tool, MirTarget, which was developed by analyzing thousands of miRNA-target interactions from high-throughput sequencing experiments. MiRTarBase is a database of experimentally validated microRNA targets. TargetScan predicts the biological targets of miRNAs by searching for the presence of conserved 8mer, 7mer, and 6mer sites that matched the seed region of each miRNA.

In addition, genes present in all three databases were regarded as the target genes of these miRNAs. Comparing predicted target genes with DEGs, only the remaining overlapped genes and their interaction pairs were used for constructing mRNA-miRNA-lncRNA pairs.

mRNA-miRNA-lncRNA pair topology analysis and stability analysis

The network topology properties were analyzed using the plug-in Network Analyzer Tool Kit in Cytoscape. Network Analyzer network analysis consisted mainly of network diameter, number of connected components, average clustering coefficient, distribution map, and so on. The connection nodes, the length of the path, and the closeness centrality of the node were calculated.

A ceRNA network incorporating TF construction

The overlapped genes of PPI network and mRNA-miRNA-lncRNA pairs were selected. Venn diagrams online software was used to visualize the overlapped genes. At the same time, a ceRNA based on these overlapped genes was constructed. In addition, iRegulon software was used to predict the TFs of mRNA and miRNA in ceRNA. Then, a ceRNA network incorporating TFs was constructed successfully.

Results

Data preprocessing

After screening using rigorous criteria, adapter reads, reads with the N ratio >10%, all A base reads, and low-quality clean reads (the number of bases with the quality value Q ≤ 20 accounting for more than 50% of the entire reads) were removed. The base composition in the data appeared more balanced and of higher quality. This indicated that the analysis was more accurate. In addition, the statistical results of the comparison region showed that most of the sequencing reads were distributed in the exon region. Furthermore, the randomness analysis and sequencing saturation analysis results showed that reads were evenly distributed in all parts of the genes, and the number of detected genes tended to be saturated. These results facilitated an accurate analysis (Supplementary Figs. S1 and S2A-L).

lncRNA analysis

The statistics of new lncRNA transcript types after reconstructing the transcripts are shown in Supplementary Figure S2M. A total of 25,806 lncRNAs (FPKM) were detected. The expression distribution results showed no significant difference between any two samples (Supplementary Fig. S3A). The results of PCA and Pearson correlation coefficient analysis showed that each group had high similarity and low variation among its samples (Supplementary Fig. S3B and Fig. 1B). In addition, seven DELs with symbols were screened (Fig. 2). The heatmap of the lncRNAs showed that EM-related infertility clustered separately from the controls (Fig. 3).

FIG. 1.

Workflow of this study and data processing. (A) Workflow of this study. (B, C) Sample correlation heatmap.

FIG. 2.

Volcano plot. Volcano plots of differentially expressed mRNA and lncRNA between patients with EM-related infertility and controls. The values shown on the x-axis and y-axis in the scatter plot are the normalized signal values of each sample (log2 scale). mRNAs; EM, endometriosis; FC, fold-change.

FIG. 3.

Hierarchical clustering heatmap. Hierarchical clustering heatmap of DEGs and DELs.

mRNA analysis

In this study, 19,684 mRNAs were detected; the expression distribution results showed no significant difference between any two samples (Supplementary Fig. S3C). Each group had high similarity and low variation among its samples (Supplementary Fig. S3D and Fig. 1C). Finally, 46 DEGs were screened through differential analysis (Fig. 2). The heatmap of the mRNAs showed that EM-related infertility clustered separately from the controls (Fig. 3).

Construction of the PPI network

A total of 18 nodes and 21 edges were mapped in the PPI network, and the top 5 genes, with the highest interaction degrees, including N4RA2, N4RA1, N4RA3, KDR, and FOSB, were identified (Fig. 4A).

FIG. 4.

PPI network and ceRNA network construction. (A) Construction and analysis of the PPI network. (B) Sankey diagram for the ceRNA network. Each rectangle represents a gene, miRNA, or lncRNA, and the connection degree of each RNA is visualized based on the size of the rectangle. ceRNA, competing endogenous RNA; PPI, protein-protein interaction.

Construction of mRNA-miRNA-lncRNA pairs

In this study, mRNA-miRNA-lncRNA pairs, including three lncRNAs (XIST, MCM3AP-AS1, and C17orf51), seven miRNAs (hsa-miR-125a-5p, hsa-miR-125b-5p, hsa-miR-17-5p, hsa-miR-20b-5p, hsa-miR-217, hsa-miR-363-3p, and hsa-miR-449c-5p), and five mRNAs (BTG2, SCAMP5, KIAA1191, SCAMP5, and NR4A2), were successfully built (Fig. 4B).

Topology and stability analyses of mRNA-miRNA-lncRNA pairs

The topological analysis was conducted on the constructed mRNA-miRNA-lncRNA pairs. Figure 5A shows the node degree distribution. It was obvious that the maximum degree of the node was 3, and the lowest was 1. The closeness (centrality, CC) of the node was computed as the number of connection lengths between the specified node and other nodes; more concentrated nodes had higher scores and, therefore, the closeness centrality meant the shortest path. Figure 5B shows some nodes with the same number of neighbors, and they were relatively concentrated. In contrast, some nodes with higher connectivity were relatively scattered.

FIG. 5.

Topology analysis and stability analysis of mRNA-miRNA-lncRNA pairs. (A) Node degree distribution of the ceRNA network. (B) Closeness centrality distribution map of the ceRNA network. (C) Shortest path distribution map of the ceRNA network. (D) Degree distribution density map of nodes. (E) Venn diagram. (F) A ceRNA network combining TFs. TFs, transcription factors.

The path can reflect the combination of all nodes in the network. Figure 5C shows the shortest path length distribution of the network and demonstrates that the path length distribution was more concentrated. The extreme values at both ends were less. The upper limit was 4. The lower limit was 1, indicating that most nodes in the network could be connected through a shorter path correlation.

The degree distribution density map of nodes is shown in Figure 5D. It was obvious that as the degree of nodes increased, the number of nodes became increasingly less, suggesting that most nodes in the network showed isolation. In the course of the occurrence of a disease, it might be that a small number of key nodes changed and their interaction affected neighboring nodes, which generated coexpression, in turn affecting downstream BP.

Construction of a ceRNA network incorporating TF

An intersection existed between the mRNA-miRNA-lncRNA pairs and PPI network, which was NR4A2 (Fig. 5E). In addition, it was obvious that the TFs with the top 5 of NES, including SRF, POLR2A, NRF1, MNT, and TCF7L2, all regulated NR4A2 in the predicted results of TFs of mRNAs and miRNAs. Finally, a ceRNA network incorporating TFs (Fig. 5F), including one mRNA (NR4A2), one miRNA (hsa-miR-217), three lncRNAs (XIST, MCM3AP-AS1, and C17orf51), and five TFs (SRF, POLR2A, NRF1, MNT, and TCF7L2), was constructed. In this study, NR4A2 and MCM3AP-AS1 were downregulated, whereas XIST and C17orf51 were upregulated in EM-related infertility.

Discussion

EM is a female hormone-dependent disease characterized by the ectopic placement of EM outside the uterine lumen. About 50% of patients with EM have infertility (Khine et al., 2016; Eisenberg et al., 2018). The mechanisms underlying infertility in EM have not been fully elucidated, but they might be associated with follicular development (D'Hooghe et al., 2019). GCs, as the main cells of follicle development, can synthesize and secrete estrogen and progesterone, promote follicle development, and provide favorable conditions for follicle maturation (Amar et al., 2020; Chou et al,. 2020). At the same time, they can also provide nutrients and energy to oocytes, promote the development and maturation of oocytes, and maintain the normal functions of oocytes. The proliferation, differentiation, and apoptosis of GCs are also the key factors affecting follicle development and oocytes. Therefore, the normal function of GCs is crucial to the growth, development, and maturation of oocytes (Yamamoto et al., 2015; Guo et al., 2016). Existing studies showed that the follicular development potentially decreased, the apoptotic rate of GCs increased, and the changes in follicular microenvironment and endocrine occurred in EM-related infertility (Gupta et al., 2008; Behera et al., 2016). Poor quality of oocytes might be closely related to GCs in follicles in patients with EM-related infertility.

Compared with controls, patients with EM-related infertility showed decreased follicular growth rate, decreased number of dominant follicles, low oocyte yield, relatively increased number of empty follicles and atresia follicles, poor oocyte quality, and decreased fertilization rate (Ojima et al., 2019). The DNA analysis of follicular GCs in infertile patients by flow cytometry showed that more GCs in patients with EM were in the apoptotic phase and S phase. In contrast, the number in the G2/M phase decreased compared with that in infertile patients with other factors (Toya et al., 2000). In addition, infertility caused by EM is also associated with increased oxidative stress in GCs (Augoulea et al., 2009). These changes result in the disruption of the cell cycle of normal GCs, which subsequently results in empty follicles, lower egg quality, and lower fertilization rates (Sun et al., 2018). Therefore, this study attempted to reveal the differences in GCs from patients with EM.

In this study, total RNA of ovarian GCs was extracted from patients with EM-related infertility and controls. Then, a specific ceRNA regulatory network combining TFs was constructed by RNA sequencing and comprehensive bioinformatics methods. A total of 25,806 lncRNAs and 19,684 mRNAs were detected, and 7 DELs and 46 DEGs were screened out. The overlapped gene of mRNA-miRNA-lncRNA pairs and PPI network was NR4A2. Finally, a ceRNA network combining TFs, including one mRNA (NR4A2), one miRNA (hsa-miR-217), three lncRNAs (XIST, MCM3AP-AS1, and C17orf51), and five TFs (SRF, POLR2A, NRF1, MNT, and TCF7L2), was constructed. Importantly, NR4A2 and MCM3AP-AS1 were downregulated in EM-related infertility, whereas XIST and C17orf51 were upregulated in EM-related infertility. Among the ceRNA regulatory network incorporating TFs, three lncRNAs (XIST, MCM3AP-AS1, and C17orf51) were associated with one miRNA (hsa-miR-217). It was proposed that the possible competition of XIST, MCM3AP-AS1, and C17orf51 in binding to hsa-miR-217 influenced the downstream regulation of NR4A2. It was worth noting that NR4A2 was regulated by 5 TFs.

NR4A2 (also called Nurr1, NOT, TINUR, or NGFI-Bβ), an orphan nuclear receptor, is rapidly and strongly induced by oxidized lipids, cytokines, and pathogen-associated molecular patterns (Mahajan et al., 2015). Abnormal expression of NR4A2 has a correlation with many inflammatory disorders (Chen et al., 2018a; Liu et al., 2018a; Wirth et al., 2020). Increased expression and function of NR4A2 have been reported to be related to a variety of signaling pathways, such as activation of the p53 pathway, which decreased NR4A2 expression; the p38 MAPK pathway, which could modulate the activity of Nurr1; and the Wnt signaling pathway inhibited by Nurr1 and PI3K-Akt-mTOR pathway (Beard et al., 2016). EM is a disease associated with chronic inflammation; the levels of inflammatory factors, such as IL-1, IL-6, IL-8, IL-10, and TNF-α, significantly increased in the follicular fluid of infertile patients with EM (Moberg et al., 2015; De Andrade et al., 2017). However, high levels of inflammatory cytokines, IL-1, IL-6, IL-8, and TNF-α, affect GCs and their function (Li et al., 2019). These inflammatory factors can inhibit the proliferation of GCs and induce the apoptosis of GCs by regulating the expression of p53, Bax, caspase, Bcl2, and survivin (Liu et al., 2018b; Haraguchi et al., 2019). These findings indicated that NR4A2 might be related to infertility caused by EM. However, no relevant report on the relationship between the two is available.

Previous studies showed that the occurrence and development of EM were related to the imbalance of miRNA. miR-217 found in the constructed ceRNA network combining TFs was a newly identified miRNA molecule with low expression in many cancers; it affected the ability of cell proliferation, migration, and invasion (Li et al., 2013; Xia et al., 2013; Li et al., 2016; He et al., 2019). However, whether miR-217 is involved in EM-induced infertility is not known and needs to be explored.

The mechanism underlying the occurrence of EM-related infertility still needs further elucidation. The findings of the present study improved the understanding and provided a broader prospect for improving the internal environment and treating EM-related infertility.

Conclusions

The proposition of a ceRNA network and the prediction of TFs in GCs from patients with EM-related infertility revealed differences in GCs. Importantly, the novel TFs, lncRNAs, miRNAs, and mRNAs involved in the ceRNA network incorporating TFs in this study might provide new insights into the underlying molecular mechanisms of EM-related infertility.

Ethics Approval

All procedures were performed according to standard protocols or manufacturers' instructions, and the study was approved by the Ethics Committees of the First Affiliated Hospital of Xiamen University.

Consent to Participate

Written informed consent was obtained from individual or guardian participants.

Consent for Publication

All authors agreed to publish the article.

Availability of Data and Materials

The expression data associated with this study were sequenced by endometrial transcriptome microarray in six volunteers.

Code Availability

Raw data and codes are available on request.

Footnotes

Authors' Contributions

R.F.W. and J.Z.L. contributed equally to this study. They mainly achieved all the experimentation. J.J.L. and N.Q.Z. were involved in sample collection and supervision. L.L.R. and W.D.Z. participated in data analysis and discussion. Q.H.C. and Y.Z.L. organized and supervised the project. All authors read and approved the final article.

Author Disclosure Statement

No competing financial interests exist.

Funding Information

This study was financially supported by the National Science Foundation of China (No. 81701419 and No. 81871145), the Natural Science Foundation of Fujian Province (2019J01566 and 2019J01565), and the Key Medical and Health Program of Xiamen (3502Z20209001).

Supplementary Material

References

Allegra

, Musolino

, Tonacci

, et al. (2020) Interactions between the microRNAs and microbiota in cancer development: roles and therapeutic opportunities. Cancers (Basel), 12:805.

Amar

, Binet

, Teteau

, et al. (2020) Bisphenol S impaired human granulosa cell steroidogenesis in vitro. Int J Mol Sci, 21:1821.

Augoulea

, Mastorakos

, Lambrinoudaki

, et al. (2009) The role of the oxidative-stress in the endometriosis-related infertility. Gynecol Endocrinol, 25:75-81.

Beard

, Tenga

, Hills

, et al. (2016) The orphan nuclear receptor NR4A2 is part of a p53-microRNA-34 network. Sci Rep, 6:25108.

Behera

, Sharma

, Singh

, et al. (2016) Relationship between endometritis and oxidative stress in the follicular fluid and luteal function in the buffalo. Reprod Domest Anim, 51:844-847.

Chen

, Zhang

, Liang

, et al. (2018a) Regulatory effects of Nr4a2 on Th2 cells from patients with pemphigus vulgaris. Oncotarget, 9:11258-11267.

Chen

, Zhou

, Chen

, et al. (2018b) fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 34:i884-i890.

Chou

, Chen

, et al. (2020) Exposure to mono-n-butyl phthalate in women with endometriosis and its association with the biological effects on human granulosa cells. Int J Mol Sci, 21:1794.

De Andrade

, Nacul

, Dos Santos

, et al. (2017) Circulating and peritoneal fluid interleukin-6 levels and gene expression in pelvic endometriosis. Exp Ther Med, 14:2317-2322.

10.

D'Hooghe

, Fassbender A

F O

D, Vanhie

. (2019) Endometriosis biomarkers: will codevelopment in academia-industry partnerships result in new and robust noninvasive diagnostic tests?. Biol Reprod, 101:1140-1145.

11.

Dumesic

, Meldrum

, Katz-Jaffe

, et al. (2015) Oocyte environment: follicular fluid and cumulus cells are critical for oocyte health. Fertil Steril, 103:303-316.

12.

Eisenberg

, Weil

, Chodick

, et al. (2018) Epidemiology of endometriosis: a large population-based database study from a healthcare provider with 2 million members. BJOG, 125:55-62.

13.

Guo

, Yang

, Liu

, et al. (2016) Effects of cumulus cell removal time during in vitro fertilization on embryo quality and pregnancy outcomes: a prospective randomized sibling-oocyte study. Reprod Biol Endocrinol, 14:18.

14.

Gupta

, Goldberg

, Aziz

, et al. (2008) Pathogenic mechanisms in endometriosis-associated infertility. Fertil Steril, 90:247-257.

15.

Haraguchi

, Hirota

, Saito-Fujita

, et al. (2019) Mdm2-p53-SF1 pathway in ovarian granulosa cells directs ovulation and fertilization by conditioning oocyte quality. FASEB J, 33:2610-2620.

16.

, Wang

, Tang

, et al. (2019) miR-217 inhibits proliferation, migration, and invasion by targeting SIRT1 in osteosarcoma. Cancer Biother Radiopharm, 34:264-270.

17.

Jin

, Yang

, Zhou

, et al. (2018) Profiles for long non-coding RNAs in ovarian granulosa cells from women with PCOS with or without hyperandrogenism. Reprod Biomed Online, 37:613-623.

18.

Khine

, Taniguchi

, Harada

(2016) Clinical management of endometriosis-associated infertility. Reprod Med Biol, 15:217-225.

19.

Kim

, Langmead

, Salzberg

(2015) HISAT: a fast spliced aligner with low memory requirements. Nat Methods, 12:357-360.

20.

Kong

, Zhang

, Ye

, et al. (2007. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res, 35:W345-W349.

21.

Langmead

, Salzberg

(2012) Fast gapped-read alignment with Bowtie 2. Nat Methods, 9:357-359.

22.

Lessey

, Gordts

, Donnez

, et al. (2018) Ovarian endometriosis and infertility: in vitro fertilization (IVF) or surgery as the first approach?. Fertil Steril, 110:1218-1226.

23.

, Zhao

, Zhang

, et al. (2013) MicroRNA-217, down-regulated in clear cell renal cell carcinoma and associated with lower survival, suppresses cell proliferation and migration. Neoplasma, 60:511-515.

24.

, Li

, Zhang

(2016) Tumor suppressor role of miR-217 in human epithelial ovarian cancer by targeting IGF1R. Oncol Rep, 35:1671-1679.

25.

, Li

, Ouyang

, et al. (2019) Investigating the impact of local inflammation on granulosa cells and follicular development in women with ovarian endometriosis. Fertil Steril, 112:882-891 e881.

26.

Liu

, Liu

, Shi

, et al. (2018a) NR4A2 protects cardiomyocytes against myocardial infarction injury by promoting autophagy. Cell Death Discov, 4:27.

27.

Liu

, Zuo

, Wang

, et al. (2018b) Ethanol promotes apoptosis in rat ovarian granulosa cells via the Bcl-2 family dependent intrinsic apoptotic pathway. Cell Mol Biol (Noisy-le-grand), 64:118-125.

28.

Love

, Huber

, Anders

(2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol, 15:550.

29.

Mahajan

, Saini

, Chandra

, et al. (2015) Nuclear receptor Nr4a2 promotes alternative polarization of macrophages and confers protection in sepsis. J Biol Chem, 290:18304-18314.

30.

Moberg

, Bourlev

, Ilyasova

, et al. (2015) Endometrial expression of LIF and its receptor and peritoneal fluid levels of IL-1alpha and IL-6 in women with endometriosis are associated with the probability of pregnancy. Arch Gynecol Obstet, 292:429-437.

31.

Ojima

, Saito

, Tsuchiya

, et al. (2019) Runx3 regulates folliculogenesis and steroidogenesis in granulosa cells of immature mice. Cell Tissue Res, 375:743-754.

32.

Pallacks

, Hirchenhain

, Krussel

, et al. (2017) Endometriosis doubles odds for miscarriage in patients undergoing IVF or ICSI. Eur J Obstet Gynecol Reprod Biol, 213:33-38.

33.

Pertea

, Kim

, Pertea

, et al. (2016) Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc, 11:1650-1667.

34.

Pritchard

, Cheng

, Tewari

(2012) MicroRNA profiling: approaches and considerations. Nat Rev Genet, 13:358-369.

35.

Rosvall

, Bergstrom

(2010) Mapping change in large networks. PLoS One, 5:e8694.

36.

Sun

, Luo

, Bu

, et al. (2013) Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res, 41:e166.

37.

Sun

, Zhang

, Wang

, et al. (2018) TMCO1 is essential for ovarian follicle development by regulating ER Ca(2+) store of granulosa cells. Cell Death Differ, 25:1686-1701.

38.

Szklarczyk

, Morris

, Cook

, et al. (2017) The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res, 45:D362-D368.

39.

Toya

, Saito

, Ohta

, et al. (2000) Moderate and severe endometriosis is associated with alterations in the cell cycle of granulosa cells in patients undergoing in vitro fertilization and embryo transfer. Fertil Steril, 73:344-350.

40.

Vercellini

, Vigano

, Somigliana

, et al. (2014) Endometriosis: pathogenesis and treatment. Nat Rev Endocrinol, 10:261-275.

41.

Wirth

, Mariani

, Bergant

, et al. (2020) Loss-of-function mutations in NR4A2 cause dopa-responsive dystonia Parkinsonism. Mov Disord, 35:880-885.

42.

Xia

, Ooi

, Hui

(2013) MicroRNA-216a/217-induced epithelial-mesenchymal transition targets PTEN and SMAD7 to promote drug resistance and recurrence of liver cancer. Hepatology, 58:629-641.

43.

Yamamoto

, Kuwahara

, Taniguchi

, et al. (2015) Tumor necrosis factor alpha inhibits ovulation and induces granulosa cell death in rat ovaries. Reprod Med Biol, 14:107-115.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

1.20 MB

1.11 MB

0.38 MB