GCs collected from infertile women due to OEM were put in the OEM group, while GCs from women in couples with male factor (MF) infertility were regarded as control. The clinical characteristics of GC samples for RNA-seq, detailed in Table 1, showed that the only significantly different parameter between the two groups was the number of retrieved oocytes (p < 0.05).

The quality control statistics for the RNA-sequencing is presented in Supplementary Table S1. The RNA quality detection data of GCs for RNA-seq samples are presented in Supplementary Table S2 and Supplementary Figures S1A–C. The expression of cell-type-specific markers for GCs is shown in Supplementary Figures S1D–O. There were 15,890 mRNAs detected in total, of which 459 significantly DEGs (396 upregulated and 63 downregulated; Figure 1A and Supplementary Table S3) were found in the GCs by comparing the two groups. Z-score normalized fragments per kilobase of transcript per million mapped reads (FPKM) of differentially expressed genes were used as inputs for two-dimensional principal component analysis (PCA) by using R package Stats (Figure 1B). There was a clear separation between the distribution of the GC mRNA datasets from the OEM and MF group (p < 0.05).

To clarify the correlation between the expression of DEGs and the clinical characteristics of samples, we constructed their correlation heatmap (Figure 1C). The results showed that most of the DEGs were significantly correlated with the number of retrieved oocytes, followed by 2PN zygote numbers, cleavage stage embryo numbers, and AMH level.

Functional enrichment analysis of these DEGs helped to identify the pathways that affected the function of GCs from patients with OEM. GO analysis demonstrated that these DEGs were mainly involved in the following terms: cell junction, MAPK cascade, oocyte maturation, responses to FSH and LH, and cell cycle (Figure 2A). By using KEGG analysis, we found that these DEGs were involved in several pathways, including adhesion junction, Wnt signaling, MAPK signaling, FoxO signaling, cAMP signaling, and estrogen signaling pathway (Figure 2B).

In order to screen more potential DEGs for validation, we included the mRNA with threshold of |log₂(fold change (FC))| ≥ 1, p < 0.08, and FPKM >0.8 in all samples; then, 77 mRNA were selected. Among them, 60 genes were upregulated, and 17 genes were downregulated simultaneously. To validate the transcriptomic sequencing results, genes that were previously reported related to GCs/oocyte/follicular development, genes that intersected with reported original transcriptomic sequencing data for GCs (19) and OEM-affected oocytes (20), and genes that were involved in significantly enriched GO functions of this study were selected from the significant 77 DEGs and tested on more clinical GC samples (Figure 3A and Supplementary Table S4). These included up- and downregulated genes demonstrated in Figure 3B.

The clinical characteristics of patients’ cohort for validation are listed in Table 2. The average AMH, basal antral follicle count, retrieved oocytes number, growing follicle number, 2PN zygote number, cleavage stage embryo number, high qualified embryo number, and 2PN fertilization rates of OEM group were also significantly lower than those in MF group, while basal FSH level was also higher. These clinical manifestations are consistent with the phenomenon that OEM leads to the ovarian reserve dysfunction and compromised oocyte quality.

The quantitative PCR (qPCR) results of GCs for validation showed that the relative expressions of NR5A2 (Figure 4A), MAP3K5 (Figure 4B), PGRMC2 (Figure 4C), DEPTOR (Figure 4D), ITGAV (Figure 4E), KPNB1 (Figure 4F), PRKAR2A (Figure 4G), GPC6 (Figure 4H), EIF3A (Figure 4I), and SMC5 (Figure 4J) were indeed significantly higher in OEM group when compared with the MF group, while DUSP1 (Figure 4K) was lower expressed in OEM group.

We analyzed the correlations of the qPCR value with the IVF clinical outcomes, including the numbers of retrieved oocytes, 2PN zygote, cleavage stage embryo on day 3 of culture (developed from 2PN zygotes), and growing follicles (diameter ≥ 14 mm on the day of hCG triggering), in validation cohort. Our results demonstrated that, among the validated DEGs, the number of retrieved oocytes was negatively correlated with the relative gene expression levels of NR5A2 (Figure 4L), PGRMC2 (Figure 4M), DEPTOR (Figure 4N), ITGAV (Figure 4O), KPNB1 (Figure 4P), PRKAR2A (Figure 4Q), GPC6 (Figure 4R), EIF3A (Figure 4S), and SMC5 (Figure 4T) and positively correlated with the relative expression of DUSP1 (Figure 4U). In addition, the number of 2PN zygote (Supplementary Figures S2A–J), cleavage stage embryo (Supplementary Figures S3A–J), and growing follicles (Supplementary Figures S4A–G) were correlated with the relative gene expression levels of the validated DEGs in this study.

Through RNA-seq, we detected 47,414 non-protein coding transcripts from these GC samples. Among these transcripts, there were 14,379 known lncRNA annotated in the genome, and 33,035 novel lncRNAs were predicted. Differential expression analysis demonstrated that 35 known lncRNAs were significantly upregulated including PTPRG-AS1, LINC01278, MEG8, LERFS, PRKCQ-AS1, and SNHG17, and 24 known lncRNAs were significantly downregulated including MALAT1, FTX, BDNF-AS, DANCR, and GASAL1 (Figure 5A and Supplementary Table S3). The correlation heatmap showed that most of the differentially expressed known lncRNA were also related to the number of retrieved oocytes, followed by 2PN zygote numbers and cleavage stage embryo numbers (Figure 5B). Interestingly, the correlation network diagram showed that most of significantly differentially expressed known lncRNAs were negatively correlated with these clinical indicators. For example, MRPL20-AS1 and MIF-AS1 were all negatively correlated with the number of retrieved oocytes, 2PN zygote numbers, and cleavage stage embryo numbers (Figure 5C).

For the predicted novel lncRNA, 360 of them were significantly upregulated, and 51 were significantly downregulated (Figure 5D and Supplementary Table S3). The correlation heatmap also showed that most of the differentially expressed novel lncRNAs were also related to the number of retrieved oocytes, followed by 2PN zygote numbers and cleavage stage embryo numbers (Figure 5E). Meanwhile, through screening, in the novel lncRNA whose average FPKM of all samples was ≥1, we found that most of these lncRNA were negatively correlated with these clinical characteristics. Novel lncRNA, such as MSTRG.47475.1, MSTRG.51840.1, MSTRG.47454.1, and MSTRG.14101.4 were all negatively expressed with the number of retrieved oocytes, 2PN zygote numbers, and cleavage stage embryo numbers (Figure 5F).

Establishing mRNA–lncRNA regulatory networks is useful for revealing interactions between mRNAs and co-expressed lncRNAs. We identified 210 significant mRNA–lncRNA(known) groups with correlation coefficient threshold |abs_rho| > 0.5 and p < 0.05. Most of them were (208/210) positively associated, and only two of the groups were negatively associated. The validated DEGs, including NR5A2, MAP3K5, PGRMC2, ITGAV, KPNB1, PRKAR2A, GPC6, EIF3A, and SMC5 co-expressed with known lncRNAs under the correlation coefficient threshold |abs_rho| > 0.9 and p < 0.05 (Figure 6A). Among them, EIF3A, PGRMC2, ITGAV, SMC5, and PRKAR2A had over 10 paired co-expressed known lncRNA (Supplementary Table S5).

Simultaneously, we identified 750 mRNA–lncRNA(novel) groups with correlation coefficient threshold |abs_rho| > 0.5 and p < 0.05 from the novel lncRNA with average FPKM value of all samples ≥ 1. Most of them were (744/750) positively associated, while only six of the groups were negatively associated. Identically, under the correlation coefficient threshold |abs_rho| > 0.9 and p < 0.05, validated DEGs including NR5A2, MAP3K5, PGRMC2, ITGAV, KPNB1, PRKAR2A, GPC6, EIF3A, and SMC5 had co-expressed novel lncRNAs, while DEPTOR and DUSP1 had not (Figure 6B). Among them, GPC6, NR5A2, EIF3A, and MAP3K5 had over 20 pairs of co-expressed novel lncRNA (Supplementary Table S5).

The Venn diagrams demonstrated that MIF-AS1, LERFS, AC123912, LINC00662, MUC20-OT1, and FP671120 were co-expressed known lncRNAs, which both correlated with clinical characteristics and validated DEGs (Figure 6C). Novel lncRNAs, such as MSTRG.15872.5, MSTRG.47475.1, MSTRG.19364.1, MSTRG.14101.4, and MSTRG.42454.1, were all associated with clinical characteristics and validated DEGs (Figure 6D).

The use of human GC samples for this study was approved by the Institutional Review Board of Sir Run Run Shaw Hospital, Zhejiang University School of Medicine (ethical approval no. 20200423-44), and signed informed consent was obtained from all participants before inclusion.

GCs were collected from infertile women who underwent IVF due to unilateral or bilateral OEM or from women in couples with MF infertility (henceforth “MF group”) from May 2018 to April 2019. The inclusion criteria for both groups of women included an age of 20–42 years, regular menstrual cycles of 25–35 days, body mass index (BMI) of 18–25 kg/m², and normal hormone profile in the early proliferative phase (days 2–5 of the menstrual cycle). Furthermore, patients from the MF group did not show signs of endometriosis or any other female infertility problems. For these patients, the main reason of infertility was diagnosed low sperm quantity and quality, abnormal chromosome, or sexual dysfunction of their husband. Simultaneously, patients from the OEM group underwent previous laparoscopic surgery for unilateral or bilateral ovarian endometriomas, with clear pathological diagnosis, and were classified as stage II–IV according to the World Endometriosis Society consensus on the classification of endometriosis (2). The detailed information about the distribution (bilaterally or unilaterally), location, and size of lesions of all OEM patients included in this study is listed in Supplementary Table S6. In addition, patients with chromosomal abnormalities, acute inflammation, malignant tumors of the reproductive system, previous endometrial neoplasias, or other anatomical endometrial pathologies were excluded.

All patients included in this study underwent controlled ovarian stimulation with gonadotropins. Oocyte retrieval was scheduled 36 h later after human chorionic gonadotropin (hCG, Livzon, Guangdong, China) was administered. When preovulatory follicles were aspirated, the oocyte–cumulus complexes were pulled from the follicle wall with the fluid. Under a microscope, all the oocyte–cumulus complexes from one single patient were identified and picked out and then washed in Gemmate medium (Cook Medical, Indianapolis, IN, USA). Before oocytes were removed for further culture prior to IVF, oocytes and cumulus cells were separated from the oocyte–cumulus complexes with syringe needles. The rest of the cumulus cells, as a subtype of GCs, were pooled and collected for RNA extraction.

Total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s procedure. The total RNA quantity and purity were analyzed using an Agilent Bioanalyzer 2100 and RNA 1000 Nano LabChip kits (Agilent Technologies, Santa Clara, CA, USA) with an RNA integrity number >7.0. An RNA-seq library was constructed following the manufacturer’s instructions using SMARTer^® Stranded Total RNA-Seq kits v. 2 (Takara Bio USA, Inc., Mountain View, CA, USA). The first-strand cDNA was synthesized from 10 ng of diluted RNA. Then, Illumina adapters and barcodes were added to the cDNA fragments. After purification of the RNA-Seq library with AMPure beads and depletion of ribosomal cDNA with ZapR v. 2 and R-Probes v. 2, the final RNA-Seq library was amplificated and enriched. The mean insert size for the paired-end libraries was 300 bp (± 50 bp).

Paired-end sequencing was performed on an Illumina X10 sequencing system with a constructed human GC RNA-seq library at LC-Bio Technology Co. Ltd. (Hangzhou, China). The average sequencing depth is 10.98 G. After quality control for reads, 5.81 G was left.

We aligned reads to the UCSC Homo sapiens reference genome GRCh38 (http://genome.ucsc.edu/) using the HISAT package (http://daehwankimlab.github.io/hisat2/). HISAT allows multiple alignments per read (up to 20 by default) and a maximum of two mismatches when mapping the reads to the reference genome.

The mapped reads of each sample were assembled using StringTie (https://ccb.jhu.edu/software/stringtie/). Then, all transcriptomes from the samples were merged to reconstruct a comprehensive transcriptome using gffcompare (github.com/gpertea/gffcompare/). StringTie (66) was used to calculate the expression levels for mRNAs and lncRNAs by calculating FPKM. The differentially expressed mRNAs and lncRNAs were selected with |log₂(FC)| ≥ 1 and with statistical significance (p < 0.05) using R package edgeR (67).

To identify novel lncRNA, transcripts that overlapped with known mRNAs and transcripts shorter than 200 bp were discarded. Then, we utilized CPC0.9-r2 (68) (cpc2.cbi.pku.edu.cn/) with default parameters (cpc2 -i novel.fa -o cpc2.out) and CNCI2.0 (69) (www.bioinfo.org/software/cnci) with default parameters (CNCI.py -f novel.fa -o CNCI.result -p 1 -m ve -g novel.gtf -d genome.fa) to predict the coding potential for transcripts. All transcripts with CPC score <−1 and CNCI score <0 were removed, and the remained transcripts were considered as novel lncRNAs.

All the significantly DEGs, including 396 upregulated and 63 downregulated genes, were subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment and Gene Ontology (GO) evaluation. “S” was the number of genes with significant difference in specific GO/KEGG categories, “TS” represented the number of genes with significant difference, “B” was the number of genes for a specific GO/KEGG, and “TB” was the total number of genes. KEGG pathways enriched and presented downregulated and upregulated genes with similar biological function. The enrichment of down- and upregulated genes for the same GO function was also demonstrated. The higher rich factor (S gene number divided by B gene number) represents the higher enrichment degree. The calculation formula of p-values of enrichment analysis was: P = 1 − ∑ i = 0 S − 1 ( i B ) ( T S − i T B − B ) ( T S T B ) . The R package Ggplot2 (70) was used to visualize the data of analysis.

Selected DEGs were corroborated through real-time quantitative polymerase chain reaction (qPCR) as described previously (6). Primers are shown in Supplementary Table S7. The validation of primers for qPCR was estimated through the melting curve of qPCR procedure. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as the endogenous control. For each sample, the average cycle threshold (Ct) was calculated from triplicate wells that varied by <0.5 Ct. All Ct values were normalized against the average Ct values of the GC samples from the MF patient group.

The Weighted Gene Co-expression Network Analysis (WGCNA) (71) was used to obtain the network relationship. According to the relationship, Cytoscape was used to draw the network diagram. Correlation network was performed using the OmicStudio tools at https://www.omicstudio.cn/tool. A Pearson’s correlation analysis was subsequently performed to cluster samples and detect outliers. The threshold of co-expression network between known/novel lncRNAs and clinical characteristics was with |correlation coefficient| > 0.5, p < 0.05. The threshold of co-expression network between known/novel lncRNAs and validated DEGs was with |correlation coefficient| > 0.9, p < 0.05.

Linear regression was applied to assess the correlation between validated DEGs and IVF clinical outcomes. Spearman correlation analysis was used to construct correlation heatmap between clinical characteristics and the expression of mRNAs and lncRNAs. Pearson correlation analysis was utilized to analyze co-expression relationships between mRNAs and lncRNAs. p < 0.05 is considered statistically significant. The data were analyzed using IBM SPSS Statistics v. 20.0 for Mac (IBM Corp., Armonk, NY, USA). Venn’s diagrams were schematized by using https://bioinfogp.cnb.csic.es/tools/venny/index.html.