Age-related differential gene expression in granulosa cells and its effects on fertility using high-throughput transcriptomics

Abstract

More couples worldwide, delay their childbearing years. The increase in age causes a gradual decrease in female ovarian function and fertility, leading to an exponential decrease in women over 35 years of age having children. Although promising for some, assisted reproductive technology (ART) is not promising for older women. Decreased fertility in advanced age has become a growing concern in the field of reproduction. In this study, high-throughput transcriptome sequencing was used to identify the differentially expressed genes (DEGs) in the ovarian granulosa cells (GCs) of older women (aged 35–44) with infertility and younger women (aged 25–34). The enriched functions and signaling pathways of DEGs were analyzed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). The function of DEGs were analyzed and predicted combined with clinical ART data. Sequencing results were verified by quantitative reverse transcription-polymerase chain reaction. Retrospective clinical data and bioinformatics analyses revealed marked reductions in the retrieved oocyte, metaphase II oocyte, 2PN fertilization, and effective embryo numbers in older women. Although the clinical pregnancy and live birth rates did not differ notably between the groups, the miscarriage rate increased significantly in older women. In total, 620 DEGs were identified, of which 246 were upregulated, and 374 were downregulated in the older group. GO, and KEGG analyses indicated that the mechanism of fertility decline in older women was probably related to chronic inflammation, cytokine receptor interaction, and oxidative stress. In conclusion, combined with basic clinical ART data and pregnancy outcomes, we tried to provide a more intuitive and in-depth understanding of age-related reduction in ovarian function and pathogenesis of infertility with regard to chronic inflammation and oxidative stress.

Keywords:

• Age-related fertility decline
• RNA-sequencing
• ovarian granulosa cells
• differentially expressed genes

Introduction

Currently, the number of married couples who tend to delay their childbearing years has increased, which is potentially associated with social situations, such as improvement in education and income levels of women, changes in social fertility policies and concept, development of oral contraceptives, or even the occurrence of losing the only child at old age (Attali and Yogev 2021). This has subsequently increased the proportion of women with birth demand at advanced age worldwide. According to statistics, the average age of women in industrialized countries giving birth to live offspring has gradually increased from 23–25 years in the early 1970s to 27–29 years in 2017 (OECD Family Database 2019). In the USA, many women postponed their plans for pregnancy between 1991 and 2001 (Heffner 2004), and the proportion of women aged 35–39 and 40–44 years at first birth increased by 36% and 70%, respectively.

Female ovarian function and fertility decrease gradually with age; in particular, the exponential decrease in fertility in women over 35 years of age increases the incidence of infertility in this age group (Sun and Zhu 2014; Vannuccini et al. 2016; Frederiksen et al. 2018). A Chinese report indicated that about 15% of women with infertility were aged under 30 years old, 15% were aged 30–34 years, 30% were aged 35–39 years and 64% were aged 40–44 years (Sun and Zhu 2014). The same trend has been observed in the USA, and age-related infertility has become an important issue in recent decades (Chandra et al. 2013). However, more women are likely to become pregnant due to the emergence and continuous development of assisted reproductive technology (ART).

ART is not always successful in older women (Hansen et al. 2013; Qin et al. 2016; Qiao and Yang 2017); problems such as a reduction in oocyte number and quality, increased risk of embryo aneuploidy (Battaglia et al. 1996; Balasch and Gratacós 2012), the occurrence of lesions on the uterus and accessory structures (American College of Obstetricians and Gynecologists Committee on Gynecologic Practice and Practice Committee 2014), increased risk of pregnancy complications, and high pregnancy loss rates, cannot be easily addressed (Gill et al. 2012; Sauer 2015; Sheen et al. 2018). In a classic French study in which healthy women underwent artificial insemination because their husbands had azoospermia, the pregnancy rate gradually decreased with an increase in age (Schwartz and Mayaux 1982). Similarly, an analysis of cases undergoing in vitro fertilization-embryo transfer (IVF-ET) in the USA showed that the live birth rate in women over 35 years of age decreased significantly to less than 20%, while that in populations who received oocytes from young (women younger than age 35), healthy female donors reached 51% irrespective of the age of the recipients (Centers for Disease Control and Prevention, American Society for Reproductive Medicine Society for ART 2012). Currently, age is recognized as an independent factor that affects the final pregnancy outcomes in women (American College of Obstetricians and Gynecologists Committee on Gynecologic Practice and Practice Committee 2014; Reproductive Medicine Committee of Chinese Medical Doctor Association 2017).

Age-related decrease in follicle number and oocyte quality is the main cause of female fertility decline, and its pathological mechanism is still unclear. Ovarian granulosa cells (GC) play a crucial role in follicular development by providing essential nutrients, growth factors, and steroid hormones to oocytes. This study used high-throughput transcriptome sequencing to screen differentially expressed genes (DEGs) in ovarian granulosa cells between older and younger women with infertility. Bioinformatics analysis and clinical ART data were combined to intuitively discuss the biological basis of age-related ovarian function decline, to probe into the mechanisms underlying the occurrence of infertility in older women, and to identify possible interventions for improving the fertility of aged women.

Results

General information

Detailed general information regarding the younger and older groups is presented in Table 1. In the older group, the age and basal follicle-stimulating hormone (FSH) levels were significantly higher than those in the younger group (p < 0.0001 and p = 0.0046, respectively), while the anti-Mullerian hormone level and the number of antral follicle count (AFC) were significantly lower (p < 0.0001). Statistically significant incidences of secondary infertility were observed in the older group (p = 0.0009). Infertility years, body mass index (BMI), basal estradiol level, and fertilization mode did not differ significantly between the two groups (P > 0.05).

Table 1. General characteristics of the two groups.

Term/Characteristics	Younger group (n = 42)	Older group (n = 47)	p Value
Age (years)	31.02 ± 2.34	40.04 ± 2.43	<0.0001
Infertility years	3.62 ± 2.01	3.62 ± 3.10	0.9971
Infertility type
Primary infertility (n%)	21 (50.00)	8 (17.02)	0.0009
Secondary infertility (n%)	21 (50.00)	39 (82.98)
BMI (kg/m^2)	23.16 ± 3.16	23.62 ± 2.92	0.4820
AMH (ng/mL)	3.01 ± 1.54	1.46 ± 0.97	<0.0001
AFC (n)	18.93 ± 4.57	11.64 ± 5.09	<0.0001
Basal FSH (IU/L)	7.68 ± 2.23	9.96 ± 4.78	0.0046
Basal E2 (pg/mL)	42.32 ± 15.44	43.48 ± 18.80	0.7544
Fertilization mode
IVF (n%)	33 (78.57)	35 (76.09)	0.7812
ICSI (n%)	9 (21.43)	11 (23.91)

BMI: body mass index; AMH: anti-Mullerian hormone; AFC: antral follicle count; FSH: follicle-stimulating hormone; E2: estradiol; IVF: in-vitro fertilization; ICSI: intracytoplasmic sperm injection.

Comparison of clinical outcomes between the younger and older groups

The number of retrieved oocytes (5.94 ± 3.77 vs. 10.07 ± 3.59%, p < 0.0001), metaphase II (MII) oocytes (4.64 ± 3.17 vs. 9.11 ± 3.52, P = 0.0265), 2PN fertilization (3.13 ± 2.36 vs. 5.38 ± 2.66, p = 0.0001), and competent embryos (1.81 ± 1.64 vs 3.07 ± 2.30, p = 0.0044) in the older age group were significantly lower than those in the younger age group, while the miscarriage rate (44.44% vs. 6.67%, p = 0.0274) in the older group was significantly higher, and all the differences were statistically significant. In addition, the number of high-quality embryos (0.51 ± 0.80 vs. 0.90 ± 1.16, p = 0.1169), the rate of clinical pregnancy rate (39.13% vs. 53.57%, p = 0.3039), and live birth (17.39% vs. 39.29%, p = 0.0877) in the older group were slightly lower, although the difference with the younger group was not statistically significant (Tables 2 and 3).

Table 2. Clinical outcomes of the two groups.

Term/Characteristics	Younger group (n = 42)	Older group (n = 47)	p Value
Retrieved oocyte number	10.07 ± 3.59	5.94 ± 3.77	<0.0001
MII oocyte number	9.11 ± 3.52	4.64 ± 3.17	0.0265
2PN fertilization number	5.38 ± 2.66	3.13 ± 2.36	0.0001
Competent embryo number	3.07 ± 2.30	1.81 ± 1.64	0.0044
High-quality embryo number	0.90 ± 1.16	0.51 ± 0.80	0.1169

MII: metaphase II; 2PN: two pronuclei.

Table 3. Pregnancy outcomes between the two groups.

Term/Characteristics	Younger group (n = 42)	Older group (n = 47)	p Value
Clinical pregnancy rate	53.57 (15/28)	39.13 (9/23)	0.3039
Miscarriage rate	6.67 (1/15)	44.44 (4/9)	0.0274
Live birth rate	39.29 (11/28)	17.39 (4/23)	0.0877

Analysis of DEGs

In total, 620 statistically significant DEGs (|log2fc| >1, p < 0.05) were screened from the two groups, among which 246 genes were up-regulated (FC> 2) and 374 genes were down-regulated (FC < 0.5) in the older group. The overall distribution of DEGs is depicted in a volcano plot (Figure 1A). Fragments per kilobase of transcript per million mapped read (FPKM) values were normalized by Z-score to represent gene expression. Significant differences between the two groups were noted following a clustering analysis (Figure 1B).

Figure 1. Volcano plot and hierarchical cluster analysis for the differentially expressed genes (DEGs) in the two groups. (A) Volcano plot for the differentially expressed genes (DEGs). The X-axis represents the fold change in gene expression in different samples and the Y-axis represents the -log10(P-value) as per the figure; red, green, and blue represent the DEGs that are significantly upregulated, down-regulated, and normally expressed, respectively. (B) Hierarchical cluster analysis of differentially expressed gene (DEG) levels. The X-axis represents samples (A_granule1 represents sample 1 in older group, and Y_granule1 represents sample 1 in younger group, etc.) and the Y-axis represents genes; upregulated genes are shown in red and downregulated genes in dark blue.

Gene ontology (GO) annotation of DEGs

Gene functions were annotated using GO enrichment analysis (Figure 2A). Among the upregulated DEGs in older group, the most enriched terms were signal transduction, inflammatory response, and immune response (GO:0007165, 0006954, 0006955) in biological process (BP), and membrane and cytoplasm (GO:0016020, 0005737) in cellular component (CC) (Figure 2B). Among the down-regulated DEGs, the most enriched terms in older groups were cell cycle and cell differentiation (GO:0007049, 0030154) in BP, and membrane and cytoplasm (GO:0016020, 0005737) in CC (Figure 2C).

Figure 2. Scatter plot and histogram of GO enrichment for differentially expressed genes (DEGs). (A) Scatter plot of Gene Ontology (GO) annotation for differentially expressed genes (DEGs). Rich Factor on the X-axis refers to the ratio of DEGs enriched in the GO term to the total number of DEGs in the GO analysis. The larger Rich Factor corresponds to higher richness of the GO term. (B) Histogram of GO enrichment for upregulated DEGs. Distribution of DEGs enriched in GO terms for biological processes (BP), cellular component (CC), and molecular function (MF). (C) Histogram of GO enrichment for down-regulated DEGs. Distribution of DEGs enriched in GO terms for biological processes (BP), cellular component (CC), and molecular function (MF).

Enrichment analysis of Kyoto encyclopedia of genes and genomes (KEGG) pathways for DEGs

KEGG analysis identified 33 biological signaling pathways that were most abundant in the DEGs, including cytokine-cytokine receptor interaction and NF-κB signaling pathway (p < 0.05) (Figure 3A), among which the upregulated DEGs in the older group were concentrated in KEGG Subclasses, including the immune system, metabolism of terpenoids and polyketides, and infectious and immune diseases (Figure 3B), and the down-regulated genes were primarily enriched in KEGG main classes, including environmental information processing, signaling molecules and interaction, and signal transduction (Figure 3C).

Figure 3. Scatter plot and histogram of KEGG enrichment for differentially expressed genes (DEGs). (A) Scatter plot of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for differentially expressed genes (DEGs). The top 20 statistically significant KEGG pathways are shown. Rich Factor on the X-axis refers to the ratio of DEGs enriched in the KEGG pathway to the total number of DEGs in the KEGG analysis. Larger Rich Factor corresponds to higher richness of the KEGG pathway. (B) Subclass histogram for the enrichment of upregulated DEGs. (C) Subclass histogram for the enrichment of downregulated DEGs.

Validation of DEGs

Quantitative reverse transcription-polymerase chain reaction (qRT-PCR) was used to validate DEGs. On the basis of sequencing results and literature, we selected four DEGs with significant differential expression for validation. As shown in Figure 4, a statistically significant increase in the IL1B and CXCL8, decrease in the GSTA1 and SEMA3A expression levels were noted in the older group compared to that in the younger subjects, which was consistent with the results of sequencing.

Figure 4. The validation of four selected differentially expressed genes between the two groups using quantitative reverse transcription-polymerase chain reaction. The Y-axis represents the relative gene expression. * refers to statistically significant differences (* for p < 0.05, ** for p < 0.01, and *** for p < 0.001).

Discussion

Decreasing fertility in older women is associated with a reduction in follicle number and oocyte quality (Perheentupa and Huhtaniemi 2009), although the specific biological mechanism remains elusive. This study performed a retrospective analysis involving both clinical data and bioinformatics analytical data, which revealed a significant decrease in the number of retrieved oocytes, MII oocytes, 2PN fertilization, and competent embryos. Although the clinical pregnancy rate and the live birth rate did not differ significantly between the two groups, we observed an increase in the miscarriage rate in women over 35 years old, indicating that the incidence of adverse pregnancy outcomes increased with age. Transcriptome sequencing was used for ovarian granulosa cells of different age groups, and 620 DEGs were identified, among which 246 genes were upregulated and 374 genes were down-regulated in the older group. GO and KEGG enrichment analysis suggested that the reduction in fertility in older women may be related to chronic inflammation, cytokine receptor interaction, and immune response.

Generally, chronic low-grade inflammation is associated with the incidence of multiple age-dependent diseases (Howcroft et al. 2013; Camell et al. 2019; Messerschmidt et al. 2019), primarily characterized by the accumulation of immune cells and production of proinflammatory cytokines (Watanabe et al. 2017). The term inflammaging has been used to describe the relationship between aging, inflammation, and diseases (Franceschi and Campisi 2014; Pence and Yarbro 2019). As we know, aging is associated with a reduced immune response, which is accompanied by multiple changes in the expression of immune-related cells and pathways. This phenomenon is known as immunosenescence (Shen-Orr et al. 2016; Schmitt et al. 2013; Lim et al. 2014). The link between low-grade chronic inflammation and age-related disease may be related to a cross-talk between an aging immune system and damaged tissue. During the aging process, the balance between proinflammatory cytokines and anti-inflammatory cytokines is shifted to the proinflammatory state (Sparkman et al. 2019), where the level of the inflammatory cytokines such as TNF-α, IL-6, IL-1β and inflammatory pathways such as NF-κB show an evident age-related increase (Furman et al. 2013; Chung et al. 2015; Bellaver et al. 2017; Yo and Rünger 2018; Kim et al. 2019; Wang et al. 2019), which is closely related to the occurrence of cardiovascular, cerebrovascular, hypofunctional, and degenerative diseases (Jiang and Cadenas 2014; Furman et al. 2017; Yao and Zhao 2018; Prince et al. 2019; Shi et al. 2019; Liao et al. 2020; Furtado et al. 2020). This is consistent with our finding that the terms associated with inflammatory response and immune response are enriched in upregulated genes in ovarian granulosa cells of advanced women.

Inflammatory factors in the ovarian microenvironment affect follicular development, oocyte maturation, and ovarian reserve (Fréour et al. 2012; Boots and Jungheim 2015; Zhang et al. 2020). The increase in the levels of proinflammatory cytokines may lead to deterioration in oocyte quality and eventually result in poor outcomes in IVF (Winger et al. 2008; Revelli et al. 2009). A study in mice has revealed considerable changes in the ovarian tissue microenvironment during reproductive aging, particularly with regard to proinflammatory phenotypes (Lliberos et al. 2021). Proinflammatory cytokines such as IL-1B and TNF-α become increasingly active in the ovaries of aged mice (Briley et al. 2016). IL1B is an immunoregulatory gene of the IL-1 family with two main monomeric isoforms, α and β, of which IL-1β is a proinflammatory cytokine responsible for the expression of IL-1B. IL-1B, primarily produced by monocytes and activated macrophages, can induce the expression of various genes involved in inflammatory reactions, along with the generation of free radicals that trigger a series of inflammatory reactions (Li 2015). A study reported that the expression of IL1B mRNA in the ovary increased with age, and chronic and mild inflammatory reaction is one of the characteristics of the aging process. It was also noted that IL-1 deficiency due to drug intervention could delay, but not reverse, age-related ovarian failure (Uri-Belapolsky et al. 2014; Dou et al. 2017). Another study (Nakamura et al. 2021) on poor fertility induced by tubal disorders showed that IL-1B regulated vascular endothelial growth factors in human fallopian tubes and that low IL-1B levels in the local environment might contribute to increased tubal vascular permeability during peri-implantation (Nasu et al. 2006). In addition, oviduct epithelial cells (OECs) isolated from aged cows were found to express higher levels of inflammatory factors, including IL-1A and IL-1B (Tanaka et al. 2016; Nakamura et al. 2019), which induced OECs to produce and release inflammatory cytokines via a feedback mechanism, induce neutrophil accumulation via CCL2, and regulate the migration and adhesion of polymorphonuclear neutrophils. Consequently, the production of inflammatory cytokines is accelerated (Gillies et al. 2018), resulting in a vicious cycle of inflammatory reactions in the fallopian tube, which is not conducive for sperm-egg binding and pregnancy. This could also be used to interpret low pregnancy and successful ART rates in older women.

Glutathione S-transferases (GSTs) are a family of multifunctional isoenzymes found in various organisms. By catalyzing the binding between reduced glutathione and various endogenous or exogenous electrophilic chemicals, GSTs can increase the hydrophilicity of toxic substances, which are then pass through the cell membrane and are excreted after catabolism, contributing to detoxification without oxidative stress damage to cellular macromolecules (Hayes and Pulford 1995; Sweeney et al. 2000; Strange et al. 2001; Polimanti et al. 2013; Charo et al. 2021). Age is an important factor that affects the expression of GSTs (Xu et al. 2018). In an animal experiment, a significant reduction in the level of the GSTA1 mRNA was observed in the heart, liver, and brain tissues of HtrA2mnd2 mice with aging phenotypes compared to that in wild type mice, indicating a reduction in anti-oxidative stress in the three organs to different extents (Zhong 2019). Another study showed that the level of the GSTA1 mRNA in the liver of 12-month-old mice was significantly higher than that in 1-month-old mice, which suggests that the body has to catabolize more toxic substances with aging. In this process, the oxidative stress load increases, and more GST enzymes are required to maintain the balance of oxidation/anti-oxidation in the body to protect it from oxidative stress damage, while high levels of GSTA1 can ensure the production of sufficient amounts of GSTs (Jiang et al. 2013; Xu et al. 2018). Our sequencing and analytical results suggested that the expression level of GSTA1 in older women with infertility was lower than that in younger individuals; hence, we rationally hypothesized that with the increase in demand for anti-oxidation with age, the level of GSTA1 may increase progressively, but eventually decreases with age, resulting in a gradual reduction in the body's anti-oxidation ability and ovarian function. The differences in the observations across studies may be related to the small sample sizes or differences in the definition of reproductive age among different species.

In summary, we investigated the DEGs associated with reduced fertility in advanced age in this study, which may bring new avenues for future research regarding the biological mechanism of age-related ovarian dysfunction and reduced reproductive outcomes. However, the limitations of this study include small sample size and the fact that the mechanisms by which the identified DEGs influence fertility were not investigated; thus, in future studies, the sample size should be increased for the further validation of the fertility-associated DEGs and the investigation of the mechanisms underlying the action of different genes should be ascertained using animal experiments.

Conclusion

We identified some DEGs associated with the decrease in fertility in older women using high-throughput sequencing in this study. Combined with clinical ART data and pregnancy outcomes, we tried to provide a more intuitive and in-depth understanding of age-related reduction in fertility with regard to chronic inflammation and oxidative stress.

Materials and methods

General information

A retrospective analysis was performed using the data of 89 women with infertility who underwent ART at the Affiliated Hospital of Shandong University of Traditional Chinese Medicine (TCM), China, because of fallopian tube-associated factors or male factors from January 2018 to June 2020. Infertility was diagnosed if the patients could not have children for more than one year despite living together, leading normal sexual life, and not using contraception measures. Cases of infertility caused by other diseases, such as polycystic ovarian syndrome, premature ovarian failure, endometriosis and adenomyosis, multiple endometrial polyps, hyperprolactinemia, thyroid dysfunction, and chromosomal abnormalities were excluded. For the clinical outcome study, subjects were assigned to either the older (35–44 years old, n = 47) or the younger group (25–34 years old, n = 42) based on the age of the women. Five samples from each group were selected randomly for sequencing analysis.

Induction of ovulation

The Gonadotropin-releasing hormone antagonist (GnRH-ant) protocol was scheduled for controlled ovarian hyperstimulation in both groups. Basic endocrine and ultrasound examinations were performed on the third day of menstruation. Recombinant human follitropin alpha solution for injection (rhFSH; GONAL-f, Merck Serono, Switzerland, S2011005) was administered at an initial dose of 150–250 IU, according to the BMI of the subjects. Cetrorelix acetate powder for injection (0.25 mg/d) (GnRH-ant, Cetrotide, Merck Serono, Switzerland, H20140476) was administered when B-ultrasound revealed dominant follicles of 13–14 mm in diameter until the trigger day, followed by injection of 250 IU recombinant human choriogonadotropin alfa solution (rhCG; OVIDREL, Merck Serono, Switzerland, S20130091) when diameter >17 mm was observed in three or more dominant follicles. Transvaginal ultrasound-guided puncture for egg extraction was performed after 36 h.

Extraction of granulosa cells

Follicular fluid was collected following oocyte retrieval and transferred to a centrifuge tube for centrifugation at 2,000 rpm for 5 min. The supernatant was discarded, and the precipitate was harvested in a round-bottom tube, which was followed by the addition of

3 mL phosphate-buffered saline (PBS) and mixing. The mixture was added slowly to the upper layer of the lymphocyte separation solution (5 mL) in a 15 mL centrifuge tube and centrifuged at 1,500 rpm for 20 min. The cells at the middle interface were aspirated into a new 15 mL centrifuge tube (containing 5 mL PBS), mixed well, and centrifuged for 5 min at 2,000 rpm. The supernatant was discarded, and the granulosa cells obtained were collected in a 1.5-mL centrifuge tube for further centrifugation at 1,500 rpm for 2 min. The supernatant was discarded, and the final products were stored at −80 °C for later use.

RNA isolation, purification, and quantification

TRIzol reagent (Invitrogen, Carlsbad, CA, USA) was used to isolate and purify total RNA according to the manufacturer's instructions. A NanoDrop ND-1000 device (NanoDrop, Wilmington, DE, USA) was used to quantify the amount and purity of each sample. An Bioanalyzer 2100 (Agilent, CA, USA) was used to assess the integrity of RNA, and samples with RIN number >7.0 were used further.

cDNA library construction

Purified poly(A) RNA was obtained after two rounds of purification using Dynabeads Oligo (dT)25-61005 (Thermo Fisher, CA, USA). Next, the poly(A) RNA was fragmented into small pieces using Magnesium RNA Fragmentation Module (NEB, cat.36150, USA) at 94 °C for 5–7 min. The cleaved RNA fragments were reverse-transcribed to generate cDNA by SuperScript™ II Reverse Transcriptase (Invitrogen, cat. 1896649, USA), which were then used to synthesize U-labeled double-stranded DNA with E. coli DNA polymerase I (NEB, cat.m0209, USA), RNase H (NEB, cat.m0297, USA) and dUTP Solution (Thermo Fisher, cat.R0133, USA). A-bases were added to the blunt ends of each strand for ligation to the indexed adapters, which harbor T-base overhangs. Single- or dual-index adapters were ligated to the fragments, and size selection was performed with AMPureXP beads (A63881, Beckman Coulter Inc, America). Following heat-labile UDG enzyme (NEB, cat.m0280, USA) treatment of the U-labeled double-stranded DNAs, the ligated products were amplified with PCR under the following conditions: initial denaturation at 95 °C for 3 min; 8 cycles of denaturation at 98 °C for 15 s, annealing at 60 °C for 15 s, and extension at 72 °C for 30 s; a final extension at 72 °C for 5 min. The average insert size of the final cDNA library was 300 bp (±50 bp).

RNA-seq and identification of DEGs

We performed 2 × 150 bp paired-end sequencing (PE150) on an Illumina NovaSeq™ 6000 (LC-Bio Technology CO., Ltd., Hangzhou, China) following the manufacturer's recommended protocol. Gene expression was represented by FPKM, and the expression of known genes in different samples was recorded. Differential expression analysis was conducted for genes that had been StringTie assembled and quantified using the R package edgeR; DEGs were defined as those with |log2 (fold change) | ≥1 and p < 0.05.

GO and KEGG enrichment analysis of DEGs

GO describes the attributes of genes and gene products in organisms with respect to three ontologies: MF, CC, and BP. In this analysis, all DEGs were mapped to each term in the GO database (http://www.geneontology.org), and the most enriched GO terms were defined by the number of DEGs in each term scaled by the hypergeometric test. The p-value of GO enrichment can be calculated as follows: $$P=1-\mathrm{ }\sum _{i=0}^{S-1}\frac{\left(\begin{array}{c}B\\ i\end{array}\right)\left(\begin{array}{c}TB-B\\ TS-i\end{array}\right)}{\left(\begin{array}{c}TB\\ TS\end{array}\right)}$$

KEGG is a database (www.genome.jp/kegg) providing pathway information that helps in understanding the high-level biological functions of genes. Significantly activated pathways of DEGs were identified using a hypergeometric test. The p-value of KEGG enrichment can be calculated as follows: $$P=1-\mathrm{ }\sum _{i=0}^{S-1}\frac{\left(\begin{array}{c}B\\ i\end{array}\right)\left(\begin{array}{c}TB-B\\ TS-i\end{array}\right)}{\left(\begin{array}{c}TB\\ TS\end{array}\right)}$$ p < 0.05 was considered significant enrichment of GO terms or KEGG pathways.

Validation of DEGs using quantitative reverse transcription-polymerase chain reaction (qRT-PCR)

qRT-PCR was performed to confirm the reliability of the sequencing result, and four significant DEGs (including two upregulated and two downregulated) were selected for validation. SYBR Green I was used to detect the expression of target genes in the samples. The absorbance of RNA was measured using an ultra-micro nucleic acid protein analyzer (Scandrop100, Germany) at 260 and 280 nm. The cDNA was synthesized in a 20 µL reaction system using the AidLab reverse transcription kit (TUREscript 1^st strand cDNA synthesis kit). The sequences of the primers used for PCR are shown in Table 4.

Table 4. Sequences of qRT-PCR primers.

Gene	Primer F/R	Sequence (5′ to 3′)	Melting temperature (Tm)
GAPDH	F	TGCACCACCAACTGCTTAGC	60
GAPDH	R	GGCATGGACTGTGGTCATGAG
CXCL8	F	CTTTCAGAGACAGCAGAG	60
CXCL8	R	ACAGTGAGATGGTTCCTT
GSTA1	F	ACCTCCAGGAGACTGCTAT	60
GSTA1	R	CGGGTGGACTCCATTCTG
IL1B	F	CCTGTACGATCACTGAACTG	60
IL1B	R	CACCACTTGTTGCTCCAT
SEMA3A	F	ACCTCCAGGAGACTGCTAT	60
SEMA3A	R	CGGGTGGACTCCATTCTG

Real-time PCR was performed using an analytic-jena-qTOWER 2.2 fluorescent quantitative PCR instrument (Germany) under the following reaction conditions:

• Step 1: 95 °C, 3 min

• Step 2: 95 °C, 10 s

• Step 3: 60 °C, 30 s + plate read

• Go to step 2, 39 cycles

Melt curve analysis (60 − 95 °C, +1 °C/cycle, holding time 4 s)

The samples were centrifuged at 6000 rpm, 4 °C, for 30 s and amplification was performed using an analytic-jena-Easycycler (Germany) quantitative PCR instrument. The relative gene expression was calculated by the qPCRsoft3.2 software automatically. Gene glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was set as the internal control, and the 2^-ΔΔCt method was used to calculate the relative gene expression. Pfaffl method was used by the software, as shown below: $$\mathit{Ratio}=\frac{{\left(1+E_{\mathit{target}}\right)}^{\mathrm{\Delta }{\mathrm{Ct}}_{\mathrm{target}\mathrm{ }(\mathrm{control}-\mathrm{expt})}}}{{\left(1+E_{\mathit{target}}\right)}^{\mathrm{\Delta }{\mathrm{Ct}}_{\mathrm{reference}\mathrm{ }(\mathrm{control}-\mathrm{expt})}}}$$

Statistical analysis

Statistical analysis was performed using the SAS Studio platform. T-tests were performed for data following normal distribution, and the mean value and standard error of the mean (SEM) were reported. For non-normally-distributed data, the Wilcoxon rank-sum test was performed, and the frequency and percentage of each variable were reported as results. Statistical significance was set at p < 0.05 for all tests.

Ethics approval

The study was approved by the Reproductive Medicine Ethics Committee of the Affiliated Hospital of Shandong TCM University (approval number CA1801(007)) and was performed in accordance with the ethical standards as laid down in the 1964 Declaration of Helsinki and its later amendments or comparable ethical standards. The study protocol was registered at the Chinese Clinical Trial Registry, registration number ChiCTR1800019604. Written informed consent for participation in this study was obtained from each patient.

Authors’ contributions

Developed the project: YY, FL; collected the data: LD, LZ, YQ; managed the data: YY, JY; analyzed the data: MJZ, SX; prepared the manuscript: CFW; Wrote the manuscript: YY; Proofread the manuscript: JYS.

Abbreviations
2PN	=	two pronuclei
ART	=	assisted reproductive technology
AFC	=	antral follicle count
BMI	=	body mass index
BP	=	biological process
CC	=	cellular component
CXCL8	=	chemokine (C-X-C Motif) Ligand
DEGs	=	differentially expressed genes
Fc	=	fold change
FPKM	=	fragments per kilobase of transcript per million mapped read
FSH	=	follicle-stimulating hormone
GAPDH	=	glyceraldehyde-3-phosphate dehydrogenase
GnRH-ant	=	gonadotropin-releasing hormone antagonist
GO	=	gene ontology
GSTA1	=	glutathione S-transferase alpha 1
IL-1	=	interleukin-1
IVF-ET	=	in vitro fertilization-embryo transfer
KEGG	=	kyoto encyclopedia of genes and genomes
MII	=	metaphase II
MF	=	molecular function
OECs	=	oviduct epithelial cells
PBS	=	phosphate-buffered saline
SEMA3A	=	semaphorin 3A
TCM	=	traditional Chinese medicine

Declaration of interests

The authors report no conflict of interest.

Data availability statement

The data that support the findings of this study are available from the corresponding author, FL, upon reasonable request. Sequence data are available in figshare database, DOI: 10.6084/m9.figshare.17827307.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China under Grant 81774355; Development Pan of Medical and Health Science and Technology in Shandong Province under Grant 2017WS260; and Shandong Traditional Chinese Medicine Science and Technology Project under Grant 2020Q013.

This work was supported by the National Natural Science Foundation of China under Grant 81774355; Development Pan of Medical and Health Science and Technology in Shandong Province under Grant 2017WS260; and Shandong Traditional Chinese Medicine Science and Technology Project under Grant 2020Q013.