Skip to main content

Expression of long noncoding RNAs in the ovarian granulosa cells of women with diminished ovarian reserve using high-throughput sequencing

Abstract

Background

Infertility is a global reproductive-health problem, and diminished ovarian reserve (DOR) is one of the common causes of female infertility. Long noncoding RNAs (lncRNAs) are crucial regulators of numerous physiological and pathological processes in humans. However, whether lncRNAs are involved in the development of DOR remains to be elucidated.

Methods

Ovarian granulosa cells (OGCs) extracted from infertile women with DOR and from women with normal ovarian reserve (NOR) were subjected to high-throughput sequencing. Comprehensive bioinformatics analysis was conducted to identify the differential expression of messenger RNAs (mRNAs) and lncRNAs. Sequencing results were validated by the selection of lncRNAs and mRNAs using real-time reverse transcription-quantitative polymerase chain reaction (RT-qPCR).

Results

Compared with the NOR group, a total of 244 lncRNAs were upregulated (53 known and 191 novel), and 222 lncRNAs were downregulated (36 known and 186 novel) in the DOR group. Similarly, 457 mRNAs had differential expression between the two groups. Of these, 169 were upregulated and 288 were downregulated. Bioinformatics analysis revealed that the differentially expressed genes of mRNA and lncRNAs were considerably enriched in “cell adhesion and apoptosis”, “steroid biosynthesis”, and “immune system”. A co-expression network comprising lncRNAs and their predicted target genes revealed the possible involvement of the “thyroid hormone signaling pathway” and “protein binding, digestion and absorption” in DOR pathogenesis. The expression of SLC16A10 was positively regulated by multiple lncRNAs. After RT-qPCR validation of seven differentially expressed lncRNAs and mRNAs, respectively, the expression of lncRNA NEAT1, GNG12, ZEB2-AS1, and mRNA FN1, HAS3, RGS4, SUOX were in accordance with RNA-sequencing.

Conclusions

We presented the first data showing that the expression profiles of lncRNA and mRNA in OGCs between NOR and DOR patients using RNA sequencing. The lncRNAs and mRNAs that we identified may serve as novel diagnostic biomarkers for patients with DOR.

Background

Infertility is a worldwide problem involving couples of reproductive ages. It will become the third major health threat after cancer and cardiovascular diseases [1]. As a serious medical problem, infertility can have destructive effects on the quality of life and marital satisfaction. As economic and social pressures increase, the age at which women have children continues to increase. The term of diminished ovarian reserve (DOR) was first described by Navot and colleagues in 1987 [2]. It is characterized by a decline in the quantity and/or quality of oocytes, with increased serum levels of follicle-stimulating hormone (FSH) and decreased serum levels of anti-Müllerian hormone (AMH) [3]. The main clinical manifestations are irregular menstruation, perimenopause syndrome, and infertility. The United States Assisted Reproductive Technology Population Survey revealed that the prevalence of DOR has increased from 19% in 2004 up to 26% in 2011 [3], and is gradually increasing as a proportion of infertile women. Studies have demonstrated that DOR typically develops into premature ovarian failure (POF) in 1–6 years [4].

Women with DOR usually result in poor pregnancy outcomes and present major challenges in reproductive medicine, particularly receiving in vitro fertilization/intracytoplasmic sperm injection (IVF/ICSI) treatment [5]. However, the etiology and pathogenesis of DOR are not clear. The key etiological factors are age, genetics, iatrogenicity (surgery, radiotherapy, or chemotherapy), and immune and environmental factors [6]. Therefore, to improve pregnancy outcomes in women with DOR, the clinicians should focus on the molecular mechanisms underlying DOR occurrence. In this way, timely prediction of the diagnosis, necessary interventions, and therapeutic measures through identification of clinical biomarkers can be undertaken.

Long noncoding RNAs (lncRNAs) are functional RNA molecules of > 200 nucleotides in length that cannot be translated into proteins. Recently, studies have revealed that lncRNAs regulate the expression of genes in terms of transcriptional, posttranscriptional, and epigenetic modifications [7]. These actions may affect the growth, development, aging, and disease progression of organisms [8]. Additionally, lncRNAs are involved in the regulation of stem cell differentiation, embryonic development, cell proliferation, apoptosis, cell metabolism, and immune reaction [9, 10]. Recent studies have shown differential expression of lncRNAs in several gynecological disorders, such as endometriosis [11], premature ovarian insufficiency (POI) [12], and polycystic ovary syndrome (PCOS) [13]. Although several lncRNAs have been found to be involved in the abovementioned ovarian diseases, the association between lncRNAs and DOR is largely unknown.

Technologies based on transcription, gene chips, and high-throughput sequencing, have found that various lncRNAs are related to the occurrence, development, and treatment of diseases [14]. The use of biological methods to study differences in the expression of lncRNAs between patients and healthy individuals and prescreening lncRNAs by high-throughput sequencing may aid in the discovery of new targets related to many diseases. In-depth study of the biomolecular mechanisms of lncRNAs can help explain the complexity and pathological states of diseases and support the diagnosis and treatment of refractory diseases.

In this study, we aimed to employ transcriptome analysis to investigate the biological role of lncRNAs during DOR progression. Ovarian granulosa cells (OGCs) from women with DOR and those from women with normal ovarian reserve (NOR) were obtained to perform RNA-sequencing (RNA-seq) and analyze the expression profiles between two groups. Bioinformatics analysis was undertaken to identify differentially expressed (DE) lncRNAs and messenger (m) RNAs and their associated functions and signaling pathways. In addition, we applied real-time reverse transcription-quantitative polymerase chain reaction (RT–qPCR) to verify the differential expression of seven lncRNAs and mRNAs. These aberrantly expressed lncRNAs and mRNAs may serve as potential biomarkers for the development of DOR and may introduce new routes for the diagnosis of women with POF.

Results

Baseline characteristics of patients used for sequencing

As presented in Table 1, patients in the DOR group had low AMH levels, low antral follicle count (AFC), low number of oocytes obtained and high basal FSH levels compared to the NOR group, the differences in the above characteristics being indicators of ovarian reserve (P < 0.001). No other features were significantly different (P > 0.05).

Table 1 Comparison of the baseline characteristics between the NOR group and DOR group used for sequencing

Sequencing of complementary (c) DNA libraries and transcriptome profiles of OGCs

A total of 999,067,624 raw reads were acquired using the NovaSeq™ 6000 platform (Illumina, San Diego, CA, USA). After excluding adapter contamination, undetermined bases, and low-quality bases, 902,773,866 clean reads were obtained, accounting for 135.42 GB. Although the GC content of the clean data was ~ 42%, the quality scores of clean reads were > 99.9 and > 98.83% for Q20 and Q30, respectively. In summary, these results indicated that the reliability and quality of the sequencing data were sufficient for further analyses (Table 2). Over 96% of clean reads were mapped to the reference human genome using HISAT (http://daehwankimlab.github.io/hisat2/), including uniquely mapped reads from 84.10 to 86.76% (Supplemental Table 1).

Table 2 Statistical data of the reads for 12 cDNA libraries

Identification of lncRNAs and mRNAs

To further recognize the protein-coding ability of unknown transcripts, Coding Potential Calculator (CPC; http://cpc2.gao-lab.org/) and Coding-Non-Coding Index (CNCI; https://github.com/www-bioinfo-org/CNCI/) were designed to exclude those with coding potential. As a result, a total of 45,333 novel lncRNAs were detected from 12 cDNA libraries (Supplemental Table 2). These lncRNAs were dispersed evenly across the 46 human chromosomes (Fig. 1a). With regard to the genomic location of the lncRNAs, the types of lncRNAs were class_code “i” (intraintronic transcript, 72.24%), class_code “u” (intergenic transcript, 18.37%), class_code “o” (sense transcript, 3.55%), class_code “j” (bidirectional transcript, 3.25%), and class_code “x” (antisense transcript, 2.59%). These lncRNAs showed no apparent bias for genomic location (Fig. 1b). Additionally, a total of 47,520 known lncRNAs were identified between the two groups.

Fig. 1
figure 1

Comparison of the properties of mRNA and lncRNA from women with NOR or DOR. a Genomic mapping of lncRNAs in different samples. b Location of lncRNA types on 46 chromosomes. Each circle represents a type of lncRNA and corresponds to “i”, “j”, “o”, “u”, or “x” from the outer circle to the inner circle. c The expression and amount of lncRNAs and mRNAs. d Distribution of the transcript length of lncRNAs and mRNAs. e Distribution of the exon number of lncRNAs and mRNAs. f Distribution of the ORF length of lncRNAs and mRNAs

We compared lncRNAs with mRNAs to illustrate the structural features of lncRNAs. The expression and number of lncRNAs were comparable to those of mRNAs (Fig. 1c). Furthermore, the median length of lncRNA transcripts was 1504 bp, which was shorter than the median length of 2317 bp for mRNA transcripts. These findings suggest that these lncRNAs are shorter than the mRNAs (Fig. 1d). There were fewer exons in lncRNAs (mean, 3) than in mRNAs (mean, 9). In total, 68.86% of mRNAs had ≥5 exons, whereas 77.02% of lncRNAs had ≤3 exons (Fig. 1e). Moreover, in OGCs, the predicted open reading frame (ORF) length of lncRNAs was shorter than that of mRNAs (Fig. 1f).

Expression profile of lncRNAs

We discovered that 466 lncRNAs had differential expression between the DOR group and NOR group. Of these, 244 were upregulated (53 known and 191 novel), and 222 were downregulated (36 known and 186 novel) (Supplemental Table 3, Fig. 2a-b). We also calculated the average expression of lncRNAs between the NOR and DOR groups separately, filtered out genes with fragments per kilobase of exon model per million mapped reads (FPKM) values below 0.1, and then intersected the two groups to plot a Venn diagram to show the unique and common lncRNAs between NOR and DOR (Fig. 2c). To further probe the functions of these lncRNAs, we forecasted the cis-regulated target genes of the DE lncRNAs between the DOR group and NOR group. A total of 52 lncRNAs were found to have target genes if 100 kilobase pairs (kbp) was used as the cutoff. Some lncRNAs had 2–3 target genes, and 57 probable lncRNA target genes were identified (Supplemental Table 4). Based on these cis-regulated target genes, analyses of functional enrichment were performed using the Gene Ontology (GO) (http://geneontology.org/) database. A total of 76 GO terms were enriched significantly (P < 0.05). For biological process (BP), 57 terms were found, including “stem cell proliferation (GO:0072089). For “cellular component” (CC), 6 terms were found, including “proteasome activator complex” (GO:0008537). For molecular function (MF), 13 terms were found, including “CXCR3 chemokine receptor binding” (GO:0048248) (Fig. 2d-e).

Fig. 2
figure 2

Identification and characterization of the DE lncRNAs between women in the DOR group and NOR group. a Number of DE lncRNAs with upregulated and downregulated expression. b Volcano plot of the DE lncRNAs. c Venn diagram showing the number of unique and common lncRNAs between the two groups. d Histogram of the GO-enriched terms for the target genes of the DE lncRNAs. e Scatter plot of the GO enrichment for the target genes of the DE lncRNAs. f Scatter plot of the signaling pathway enrichment for the target genes of the DE lncRNAs using the KEGG database

We used the Kyoto Encyclopedia of Genes and Genomes (KEGG) (www.genome.jp/kegg/) database to identify signaling pathways significantly enriched for lncRNA target genes. Three significantly enriched signaling pathways related to the lncRNA target genes (P < 0.05) were found, including “ferroptosis”, “fatty acid biosynthesis”, and the “JAK-STAT signaling pathway” (Fig. 2f).

We discovered the predicted outcomes of the DE lncRNAs with cis-regulated genes. We listed the first-six and last-four lncRNA–gene pairs according to the Pearson correlation coefficient. The first-six lncRNA–gene pairs were regulated in the same direction, whereas the last-four lncRNA–gene pairs were in the opposite direction (Table 3).

Table 3 Differentially expressed lncRNA–gene pairs between the DOR group and NOR group

Expression profile of mRNAs

A total of 457 DE mRNA genes were found between the OGCs of the two groups. Of these genes, 169 genes were upregulated, and 288 genes were downregulated (Fig. 3a-b). We similarly plotted Venn diagram to demonstrate the unique and common mRNAs between NOR and DOR (Fig. 3c). To further explore the functions of these DE mRNA genes, analyses of functional enrichment were performed using the GO database. Among the 457 DE mRNA genes, 3404 GO terms with functional information were enriched (Fig. 3d-e). A total of 542 GO terms were enriched significantly: 359 BP terms (e.g., ‘cell adhesion’, GO:0007155; ‘positive regulation of apoptotic process’, GO:0043065 and ‘steroid biosynthetic process’, GO:0006694), 137 MF terms (e.g., ‘oxidoreductase activity’, GO:0016491), and 42 CC terms (e.g., ‘endoplasmic reticulum lumen’, GO:0005788) (Supplemental Table 5). Analyses of signaling pathway enrichment using the KEGG database revealed that the DE mRNA genes were significantly enriched in nineteen signaling pathways, including “steroid biosynthesis”, “focal adhesion”, “ECM–receptor interaction”, “PPAR signaling pathway”, and “PI3K-AKT signaling pathway” (Fig. 3f).

Fig. 3
figure 3

Identification and characterization of the DE mRNA genes between women in the DOR group and NOR group. a Number of DE mRNA genes showing upregulated and downregulated expression. b Volcano plot of the DE mRNA genes. c Venn diagram showing the number of unique and common mRNAs between the two groups. d Histogram of the GO-enriched terms of DE mRNA genes. e Scatter plot of the GO-enriched DE mRNA genes. f Scatter plot of the signaling pathway enrichment of DE mRNA genes using the KEGG database

Coenriched GO terms of the DE lncRNAs and mRNAs

We then aimed to identify the key pathways that regulate the ovarian reserve. Hence, we identified five significantly enriched GO terms in the enrichment of the DE mRNAs and the enrichment of the lncRNA target genes. The significantly enriched GO terms were “CXCR chemokine receptor binding”, “CXCR3 chemokine receptor binding”, “positive regulation of transforming growth factor-beta production”, “regulation of bile acid biosynthetic processes”, and “cellular response to organonitrogen compound”. Three pathways were involved in BP. and the other two pathways were involved in MF (Table 4).

Table 4 Coenriched GO terms of DE lncRNAs and mRNAs

Construction of a co-expression network of DE lncRNAs and mRNAs

To further explore the potential regulatory mechanisms, we constructed lncRNA-mRNA co-expression networks and performed comprehensive analyses. After cis-regulated target genes were predicted for DE lncRNAs, 57 target genes were identified that were cis-regulated by 52 lncRNAs, both of which were DE between DOR and NOR groups. Twenty-four positively correlated lncRNA–mRNA co-expression pairs containing 24 DE lncRNAs and 19 DE mRNAs were acquired at the verge of a Pearson correlation coefficient (r > 0.40) (Supplemental Table 6). Then, a lncRNA–mRNA regulatory network comprising dysregulated lncRNAs and their cis-target mRNAs was visualized using Cytoscape v3.7.1 (https://cytoscape.org/). Solute carrier family 16 member 10 (SLC16A10) was included in this network and was regulated by five lncRNAs: MSTRG.66173, MSTRG.66167, MSTRG.66177, MSTRG.66174, and MSTRG.66170 (Fig. 4). Among them, SLC16A10 was enriched in several GO terms, including “protein binding” (GO:0005515), “thyroid hormone generation” (GO:0006590), and “cell junction” (GO:0030054). SLC16A10 was also enriched in the “thyroid hormone signaling pathway” and “protein digestion and absorption”.

Fig. 4
figure 4

The lncRNA–mRNA regulatory network. Each lncRNA is a red round rectangle. Each mRNA is a green ellipse

Verification of the DE lncRNAs and mRNAs by RT–qPCR

To verify the RNA-seq results, RT–qPCR was performed on the seven DE lncRNAs and mRNAs respectively identified between the OGCs from the DOR group and NOR group. Table 5 displayed the baseline characteristics of the case used for RT-qPCR validation, with patients in the DOR group not significantly different from those in the NOR group (P > 0.05), except for indicators related to ovarian reserve such as AMH level, AFC, and number of oocytes retrieved (P < 0.05). These seven lncRNAs include nuclear rich abundant transcript 1 (NEAT1), Ras-related protein in brain 5A (RAB5A), maternally expressed gene 3 (MEG3), guanine nucleotide-binding protein subunit gamma-12 (GNG12), macrophage stimulating 1 like (MST1L), general control nonderepressible 1 (GCN1), zinc-finger E-box binding homeobox 2-AS1(ZEB2-AS1). These seven mRNAs include fibronectin 1 (FN1), hyaluronan synthase 3 (HAS3), matrix metalloproteinase 15 (MMP15), glyoxalase 1 (GLO1), regulator of G-protein signaling 4 (RGS4), sulfite oxidase (SUOX), Semaphorin 3A (SEMA3A). Table 6 presented information on the relevant sequencing results for the above lncRNAs and mRNAs.

Table 5 Comparison of the baseline characteristics between the NOR group and DOR group used for RT-qPCR
Table 6 Sequencing information of lncRNAs and mRNAs used for RT-qPCR verification

Among the lncRNAs, the expression of NEAT1, GNG12, and ZEB2-AS1 was consistent with the RNA-seq results. The expression of the remaining 4 lncRNAs was not significantly different (Fig. 5). Among the mRNAs, the expression of FN1, HAS3, RGS4, and SUOX was consistent with the RNA-seq results. The expression of the remaining 3 mRNAs did not differ significantly (Fig. 6). These results indicate that the seq-RNA data were reliable.

Fig. 5
figure 5

Expression of key lncRNAs as determined by RT–qPCR. ns, P > 0.05; *P < 0.05; **P < 0.01

Fig. 6
figure 6

Expression of key mRNAs as determined by RT–qPCR. ns, P > 0.05; *P < 0.05; **P < 0.01

Discussion

Follicular reserve and ovarian function are measures of reproductive endocrine functionality and fertility potential [15]. OGCs play important roles in the regulation of follicle activation, growth, development, and atresia [16]. Recently, several studies targeting specific genes and regulating various signaling pathways have shown that lncRNAs have vital roles in the formation and development of ovarian follicles [17,18,19].

We report, for the first time, the expression profiles of lncRNAs and mRNAs in OGCs of DOR patients and NOR patients, as identified by high-throughput sequencing. We identified 466 DE lncRNAs and 457 mRNAs associated with ovarian reserve. We showed that these dysregulated genes were involved in multiple biological processes, including cell adhesion and apoptosis, the immune system, and signaling pathways such as steroid biosynthesis. Notably, multiple lncRNAs could contribute to DOR by targeting and regulating signaling pathways mediated by SLC16A10, such as thyroid hormone synthesis and protein digestion and absorption. Overall, our results indicate the potential regulatory role of related lncRNAs and mRNAs in DOR.

In recent years, several studies have identified an important role for lncRNAs in the regulation of follicle development [20, 21]. Studies have identified differential expression of lncRNAs associated with reproductive function by using RNA-seq of OGCs obtained from humans [22], with most studies focusing on PCOS [23], endometriosis [24, 25], and POF [26, 27]. The lncRNAs we identified had fewer exons, shorter ORFs, and shorter transcript lengths compared with those in protein-coding transcripts, and these characteristics were similar to those from porcine OGCs [28] and human placentas at term [29]. These data demonstrate that the lncRNA outcomes we obtained were reliable.

Several research teams have reported aberrant lncRNA expression in POI and identified hundreds of lncRNAs associated with POI. Some of these lncRNAs could be used as biomarkers for the diagnosis or prognosis of DOR [18, 30]. The discrepancies between the DE lncRNAs identified in our study and those in other studies are due (at least in part) to differences in sample selection and sequencing methods. Moreover, several researchers have conducted functional studies on these DE lncRNAs [31,32,33] and reported different underlying molecular mechanisms. However, most of those studies concentrated on a particular gene or signaling pathway and did not analyze them from a holistic perspective. Additionally, most of the identified lncRNAs and mRNAs were associated with proliferative and apoptotic processes. The DE genes identified in our study were mostly involved in the signaling pathways related to the adhesion, proliferation, and apoptosis of cells, steroid biosynthesis and metabolism, and the immune response. All of these are biological processes involved in follicular atresia induced by OGC apoptosis [34, 35].

Folliculogenesis is regulated by multiple mechanisms involving endocrine and intraovarian signaling pathways [36, 37], and changes caused by the differential expression of genes may further affect ovarian reserve. We revealed several DE lncRNAs and mRNAs that could be associated with ovarian reserve. The lncRNAs validated by RT-qPCR were NEAT1, GNG12, and ZEB2-AS1; the mRNAs were FN1, HAS3, RGS4, and SUOX, and we will focus on the above genes. NEAT1 is a relatively well-studied lncRNA involved in the formation and maintenance of a “nuclear paraspeckle” (a nuclear body with multiple gene-expression functions). NEAT1 can promote disease through regulation of mitochondrial function [38]. NEAT1 is expressed abundantly in oocytes as well as in primordial, primary, and small antral follicles. It is involved mainly in the expression of genes related to apoptosis and extracellular matrix-related functions to achieve epigenetic control of follicular development [39, 40]. It was found that NEAT1 was DE in patients with PCOS [41] and POF [42], and that overexpression or interference with NEAT1 could improve the pathological changes in rat ovarian tissue by affecting the apoptosis of OGCs. GNG12, a particular G protein-coupled receptor, which plays crucial roles in the proliferation of tumor cells as a transducer and transmembrane signaling regulator. ZEB2-AS1 [43, 44] and GNG12 [45] can regulate cell proliferation, apoptosis, and migration processes and thus promote tumor cell growth of by affecting various pathways such as phosphoinositide 3-kinase/protein kinase B (PI3K/AKT), Wnt/β-catenin, etc. However, they have not been adequately investigated in the context of follicle growth, and the lonely study has confirmed that ZEB2-AS1 could enhance the activity of trophoblast cells and prevents the development of recurrent spontaneous abortion [46]. In the present study, we observed lncRNA GNG12 and ZEB2-AS1 expression was downregulated in the OGCs from DOR patients, suggesting it may contribute to DOR development by increasing apoptosis and inhibiting cell proliferation in OGCs.

Among the DE mRNAs, FN1 is a glycoprotein component of the extracellular matrix that is widely involved in processes such as cell migration, adhesion, and proliferation. FN1 was found to be low in PCOS patients with decreased angiogenic capacity [47]. Combined with the finding in this study that FN1 expression was decreased in women with DOR, we hypothesized that compromised vascularization around the ovarian follicles may lead to follicular growth arrest. SUOX regulated sulfite oxidase converts sulfite to sulfate in vivo, which has been found to be associated with PCOS [48] and the GO term “mitochondrion” [49]. Our study found low expression of SUOX in DOR patients, which may contribute to the development of DOR by affecting mitochondrial metabolic pathways. As for HAS3, its regulated synthesis of hyaluronan has the property of maintaining tissue homeostasis. A study has shown that there was a lower level of HAS3 in PCOS endometrium compared to women with regular menstrual cycles in the proliferative phase [50]. Similar to the findings of our study, when goat oocytes were cultured in vitro, the levels of HAS3 showed a likewise upward trend as the oocytes aged, and the number of apoptotic cells increased [51]. The variability in the above findings may be due (at least in part) to differences in species, sample selection and disease. This also suggests that HAS3 plays an important role in reproductive disorders, but the exact mechanism of action is unclear. We look forward to further studies on the role of HAS3 in influencing follicle development and oocyte quality in the future. RGS4, a GTPase-activating protein, changes in its expression levels can affect apoptosis, invasion and migration capacity and is linked to the development of diseases such as cancer [52]. It has been less studied in the reproductive system, where it has been found to accelerate the kinetics of K(+) channels in the Xenopus oocyte system [53]. This study found that RGS4 expression was down-regulated in DOR patients, this result provides a direction for future research, the exact action mechanism still needs to be further clarified.

We assessed the enrichment of functions and signaling pathways using the GO database and KEGG database, respectively, to reveal the biological functions of the DE lncRNA target genes and mRNAs. The DE genes were involved mainly in cell adhesion, steroid biosynthesis, and pathways related to immunity and apoptosis, which are the functional pathways of OGCs [34, 54, 55]. Cell adhesion is involved in the in vitro culture of porcine [56] and human [57] OGCs and can affect the proliferative potential and survival capability. Pashaiasl and collaborators measured the mRNA expression profile of the OGCs of DOR patients. Similar to our findings, they identified DE genes involved in focal adhesion [58], as well as DE genes involved in ossification, ovarian-follicle development, vasculogenesis, and sequence-specific activity of DNA-binding transcription factors. OGCs are the most important somatic cells for the synthesis of steroid hormones. We found that the DE genes could affect steroid biosynthesis and thus regulate the proliferation and apoptosis of OGCs through multiple pathways. These findings have been documented in studies in sheep [59], mice [60] and women with POF [61].

Autoimmunity is responsible for approximately 4–30% of POI cases [62]. Women suffering from POI carry a high risk of having autoimmune disorders [63]. Chemokines play an important role in reproductive immunology [64] and regulate leukocyte migration by attracting cells that express their cognate ligands. C-X-C motif chemokine ligand 12 (CXCL12), C-X-C motif chemokine receptor 4 (CXCR4), and CXCR7 can influence the development of endometriosis, Asherman’s syndrome, endometrial cancer, and ovarian cancer [65]. CXCR3, which was identified in the present study, has been investigated less thoroughly. Therefore, women should be screened for common autoantibodies (e.g., steroid cell, anti-ovarian, and anti-thyroid autoantibodies), which can aid in the early diagnosis and treatment of DOR and the prevention of POF. A prospective cohort study involving 164 women found a higher prevalence of OGC apoptosis, worse ovarian response, lower oocyte count, lower embryo count, and poor pregnancy outcomes in patients with DOR [66].

We found that the PI3K-AKT and Janus kinase/signal transducer and activator of transcription (JAK-STAT) pathways enriched by the DE genes were associated with apoptotic processes. Several studies have demonstrated that the inhibition of the PI3K-AKT pathway leads to apoptosis and autophagy in OGCs [60, 61]. The importance of the JAK/STAT signaling pathway in ovarian development and folliculogenesis in horses [67] and humans [68] has been demonstrated. Polyfluoroalkyl substances reduce the ovarian reserve and decrease endogenous hormone synthesis by activating the peroxisome proliferator-activated receptor (PPAR) pathway [69], which has been indicated as a new biomarker of follicular capacity [70]. By constructing a lncRNA–mRNA co-expression network, we found that SLC16A10 was positively regulated by multiple lncRNAs. Additionally, the function of the target gene SLC16A10 was enriched mainly in thyroid hormone synthesis as well as the binding, digestion, and absorption processes of proteins. One retrospective study found that increased levels of thyrotropic hormone may be associated with decreased serum levels of AMH [71], whereas the ovarian reserve was not associated with anti-thyroid peroxidase antibodies or anti-thyroglobulin antibodies. Transplantation of human amniotic epithelial cells has been shown to restore ovarian function and improve ovarian reserve and fertility in cyclophosphamide-induced POI rats by affecting protein digestion/absorption, and steroid-biosynthesis pathways [72].

Conclusions

The expression profiles of lncRNAs and mRNAs in OGCs from women with DOR and NOR were mined using transcriptome sequencing. A total of 466 lncRNAs and 457 mRNAs were identified to have differential expression in the transcripts. Bioinformatic analysis revealed that most of these lncRNAs and mRNAs were involved in important pathways and biological processes related to cell adhesion, apoptosis, steroid biosynthesis, and the immune system. Dysregulated lncRNAs and mRNAs could serve as diagnostic biomarkers. Further studies should concentrate on the elucidation of specific molecular mechanisms to provide a new understanding for the diagnosis and treatment of DOR.

Methods

Ethical approval of the study protocol

The research protocol was approved by the Reproductive Medicine Ethics Committee of the Affiliated Hospital of Shandong University of Traditional Chinese Medicine (Ref: SZ2018090103). All patients provided written informed consent before sample collection.

Patients and sample collection

Samples and clinical data were collected at the Center for Integrative Reproduction and Genetics, Affiliated Hospital of Shandong University of Traditional Chinese Medicine, from September 2018 to October 2020. Twenty-four patients were recruited in our study. These patients underwent IVF/ICSI using a gonadotrophin-releasing hormone antagonist for controlled ovarian hyperstimulation. All patients satisfied the following criteria: age of 25–40 years; body mass index (BMI) of 18–25 kg/m2; and no endocrine disease, uterine anomaly, endometriosis, or chromosomal abnormalities.

The study group (n = 12) comprised infertile women with DOR. The diagnostic criteria adopted an amended version of the Bologna standard [73], with participants fulfilling a minimum of two of the following conditions: (i) bilateral AFC < 6; (ii) AMH < 1.10 ng/mL; (iii) menstrual basal FSH level of 10–40 mIU/mL. The control group (n = 12) was women with NOR and pure male factor (azoospermia or severe oligospermia/aspermia) infertility. Six patient samples from each group were selected randomly for lncRNA sequencing. The other six samples underwent RT–qPCR.

Oocyte retrieval was performed 36 h after injection of human chorionic gonadotropin once follicles > 17 mm in diameter. Follicular fluid was obtained from the first aspirated follicle of each ovary during ultrasound-guided transvaginal ovarian puncture for oocyte retrieval. OGCs were collected and purified from the follicular fluid by density gradient centrifugation. Follicular fluid was centrifuged at 380×g for 5 min at 23 °C and the supernatant was removed. Phosphate-buffered saline (5 mL) was added to the precipitate and mixed. Then, 5 mL of Ficoll-Paque (GE Healthcare, Chicago, IL, USA) was added to a 15-mL centrifuge tube. The suspension was added slowly to the upper layer of Ficoll-Paque, and centrifugation at 380×g for 20 min at 23 °C was carried out. The OGC layer was aspirated and transferred to a 1.5-mL centrifuge tube, and repeated pipetting was undertaken followed by centrifugation at 380×g for 3 min at 23 °C. Finally, the supernatant was removed and the OGCs were stored in a − 80 °C for further analyses.

Preparation and sequencing of RNA libraries

TRIzol® Reagent (Invitrogen, Carlsbad, CA, USA) was employed to isolate and purify total RNA according to the manufacturer’s instructions. The amount and purity of total RNA in each sample was analyzed using a bioanalyzer (2100 series; Agilent Technologies, Santa Clara, CA, USA) and an RNA Nano1000 Assay Kit (Agilent Technologies) so that the RNA integrity number (RIN) > 7.0. Ribosomal RNA was removed according to the instructions of the Ribo-Zero™ rRNA Removal Kit (Illumina). At 95 °C, the remaining RNA was broken into short fragments using divalent cations. Subsequently, cDNA was created by reverse transcription using fragmented RNA as a template, followed by addition of deoxynucleoside triphosphate, RNase H, and Escherichia coli DNA polymerase I to synthesize the second strand. After purification of AMPure XP beads (Beckman Coulton, Fullerton, CA, USA) and end repair, poly(A) was added, and sequencing connectors were attached. Then, cDNA of a certain length range was extracted. PCR amplification was performed to obtain a cDNA library. The average insert size of the final cDNA library was 300 ± 50 bp. Twelve separate RNA-seq libraries were generated for the study group and control group.

RNA-seq

High-throughput sequencing of samples was executed on a NovaSeq 6000 using 2 × 150-bp paired-end sequences (PE150) according to a standard procedure. First, reads containing adapter contamination, low-quality bases, or unidentified bases were excluded by applying cutadapt v1.9 [74] (https://cutadapt.readthedocs.io/en/stable/#). Then, the quality of reads was authenticated using FastQC v0.10.1 [75] (www.bioinformatics.babraham.ac.uk/projects/fastqc/). The Q20, Q30, and GC contents of clean data were calculated, and subsequent analyses were based on these data. Sequencing and data acquisition were undertaken by LC-Bio Technology (Hangzhou, China).

Transcript assembly

HISAT (v2.0.4) [76] was utilized to preprocess raw data and map the processed valid data to the human reference genome (GRCh38). Annotation files in Gene Transfer Format were created to identify which genes these reads mapped to. The mapped reads from each sample were assembled using StringTie v1.3.0 [77] (https://ccb.jhu.edu/software/stringtie/index.shtml/) with default parameters. Then, all transcripts in the sample were joined using gffcompare (https://ccb.jhu.edu/software/stringtie/gffcompare.shtml) to build a whole transcriptome. Then, StringTie was used to calculate the FPKM to quantify the expression of the transcripts.

lncRNA identification

First, transcripts that overlapped with known mRNAs and transcripts shorter than 200 bp were removed. If the remaining transcripts duplicated the lncRNA, then these transcripts were deemed to be “known lncRNAs”. The rest of the non-annotated transcripts were assumed to recognize the underlying novel lncRNAs. Then, we employed CPC (0.9-r2) [78] and CNCI (v2.0) [79] to forecast the coding potential of the remaining transcripts of length ≥ 200 bp, mapped read coverage ≥3, and exon number ≥ 1. Transcripts with a CPC score < − 1 and CNCI score < 0 were assumed to be potential novel lncRNAs.

Differential expression of mRNAs and lncRNAs

StringTie was adopted to evaluate the expression of mRNAs and lncRNAs by calculation of FPKM. Multiple values were analyzed by the packages -edgeR [80] or DESeq2 [77] within R (R Institute for Statistical Computing, Vienna, Austria). P < 0.05 and | log2 (fold change) | ≥1 were used as the basis for screening DE mRNAs and lncRNAs between the DOR group and NOR group.

Target-gene prediction and functional analyses of lncRNAs

The interaction between lncRNAs and their adjacent genes is known as “cis-regulation”. We used Python scripts to predict the cis-regulatory relationships between mRNAs and lncRNAs in the 100-kbp range upstream and downstream of chromosomes. We undertook analyses of enrichment of the functions and signaling pathways of lncRNA-targeted mRNAs using GO and KEGG databases [81, 82], respectively.

Validation of lncRNA expression by RT–qPCR

Based on the analysis of sequencing results, literature review and research needs, the expression of seven lncRNAs (GNG12, MEG3, NEAT1, MST1L, RAB5A, GCN1 and ZEB2-AS1) and mRNAs (HAS3, MMP15, FN1, GLO1, RGS4, SUOX, SEMA3A) was measured by RT–qPCR in OGC samples from 12 patients to validate the RNA-seq results. Total RNA was extracted from samples using the RNeasy Micro Kit following the manufacturer’s instructions (Qiagen, Stanford, VA, USA). The content and purity of RNA were measured at an absorbance of 260/280 nm using Scandrop™ 100 (Analytik Jena, Jena, Germany). Total RNA was reverse-transcribed into cDNA using the TUREscript First Strand cDNA Synthesis Kit (Aidlab, Beijing, China) according to the manufacturer’s instructions. The specific primers designed for RT–qPCR using Primer-Blast (National Center for Biotechnology Information) and primer-primer6 together are shown in Supplemental Table 7. The analytikjena-qTOWER2.2 PCR System (Analytik Jena) was applied to conduct RT–qPCR. The reaction protocol was as follows: initial denaturation at 95 °C for 3 min, denaturing at 95 °C for 10 s, and annealing at 60 °C for 30 s for 40 cycles. Glyceraldehyde − 3-phosphate dehydrogenase, a housekeeping gene, was selected as a reference gene because of its relatively stable and high expression in ovarian granulosa cells [83]. The relative expression levels of lncRNAs and mRNAs were calculated by the 2−ΔΔCt method [84]. For each reaction, three independent biological replicates were employed.

Statistical analyses

SPSS 23.0 (IBM, Armonk, NY, USA) was used for statistical analyses. Prism 8.0.1 (GraphPad, San Diego, CA, USA) was employed for statistical analyses of RT–qPCR results and graphs. The Shapiro-Wilk test was used to assess the normality of the data distribution. Measurement data were expressed as the mean ± SD. Continuous variables were compared using the Student’s t-test or Welch’s t-test. As for categorical variables, the differences between groups were analyzed using the Chi-square test (χ2). P < 0.05 was deemed significant.

Availability of data and materials

All the raw data has been stored in GEO (Gene Expression Omnibus) with the login number GSE193136 (https://www.ncbi.nlm.nih.gov/ geo/query/acc.cgi?acc = GSE193136).

Abbreviations

AFC:

Antral follicle count

AMH:

Anti- Müllerian hormone

BMI:

Body Mass Index

BP:

Biological process

CC:

Cellular component

cDNA:

Complementary DNA

CNCI:

Coding-Non-Coding Index

CPC:

Coding Potential Calculator

DE:

Differentially expressed

DOR:

Diminished ovarian reserve

FPKM:

Fragments per kilobase of exon model per million mapped reads

FSH:

Follicle-stimulating hormone

GO:

Gene Ontology

ICSI:

Intracytoplasmic sperm injection

IVF:

In vitro fertilization

KEGG:

Kyoto Encyclopedia of Genes and Genomes

lncRNAs:

Long noncoding RNAs

MF:

Molecular function

NOR:

Normal ovarian reserve

OGCs:

Ovarian granulosa cells

ORF:

Open reading frame

PCOS:

Polycystic ovary syndrome

POF:

Premature ovarian failure

POI:

Premature ovarian insufficiency

PPAR:

Peroxisome proliferator-activated receptor

RT–qPCR:

Reverse transcription-quantitative polymerase chain reaction

RNA-seq:

RNA-sequencing

References

  1. He C, Wang K, Gao Y, Wang C, Li L, Liao Y, et al. Roles of Noncoding RNA in Reproduction. Front Genet. 2021;12:777510.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Navot D, Rosenwaks Z, Margalioth E. Prognostic assessment of female fecundity. Lancet (London, England). 1987;2(8560):645–7.

    Article  CAS  Google Scholar 

  3. Devine K, Mumford SL, Wu M, DeCherney AH, Hill MJ, Propst A. Diminished ovarian reserve in the United States assisted reproductive technology population: diagnostic trends among 181,536 cycles from the Society for Assisted Reproductive Technology Clinic Outcomes Reporting System. Fertil Steril. 2015;104(3):612–19.e3.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Balmaceda J, Schwarze J. Assessment of ovarian reserve - should we perform tests of ovarian reserve routinely? Human reproduction (Oxford, England). 2007;22(5):1492; author reply −3.

  5. Cohen J, Chabbert-Buffet N, Darai E. Diminished ovarian reserve, premature ovarian failure, poor ovarian responder--a plea for universal definitions. J Assist Reprod Genet. 2015;32(12):1709–12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Alborzi S, Madadi G, Samsami A, Soheil P, Azizi M, Alborzi M, et al. Decreased ovarian reserve: any new hope? Minerva Ginecol. 2015;67(2):149–67.

    CAS  PubMed  Google Scholar 

  7. Ponting C, Belgard T. Transcribed dark matter: meaning or myth? Hum Mol Genet. 2010;19:R162–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Kumar M, Goyal R. LncRNA as a Therapeutic Target for Angiogenesis. Curr Top Med Chem. 2017;17(15):1750–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Kong Y, Lu Z, Liu P, Liu Y, Wang F, Liang E, et al. Long Noncoding RNA: Genomics and Relevance to Physiology. Compr Physiol. 2019;9(3):933–46.

    Article  PubMed  Google Scholar 

  10. McDonel P, Guttman M. Approaches for Understanding the Mechanisms of Long Noncoding RNA Regulation of Gene Expression. Cold Spring Harb Perspect Biol. 2019;11(12).

  11. Cui D, Ma J, Liu Y, Lin K, Jiang X, Qu Y, et al. Analysis of long non-coding RNA expression profiles using RNA sequencing in ovarian endometriosis. Gene. 2018;673:140–8.

    Article  CAS  PubMed  Google Scholar 

  12. Pankiewicz K, Laudański P, Issat T. The Role of Noncoding RNA in the Pathophysiology and Treatment of Premature Ovarian Insufficiency. Int J Mol Sci. 2021;22(17).

  13. Huang X, Hao C, Bao H, Wang M, Dai H. Aberrant expression of long noncoding RNAs in cumulus cells isolated from PCOS patients. J Assist Reprod Genet. 2016;33(1):111–21.

    Article  PubMed  Google Scholar 

  14. Chen Y, Li Z, Chen X, Zhang S. Long non-coding RNAs: From disease code to drug role. Acta Pharm Sin B. 2021;11(2):340–54.

    Article  CAS  PubMed  Google Scholar 

  15. Jiang Z, Jin L, Shi W, Xi J, Hu Y, Liu X, et al. A combination of follicle stimulating hormone, estradiol and age is associated with the pregnancy outcome for women undergoing assisted reproduction: a retrospective cohort analysis. Sci China Life Sci. 2019;62(1):112–8.

    Article  CAS  PubMed  Google Scholar 

  16. Chen S, Wang Y, Liao L, Meng L, Li J, Shi C, et al. Similar Repair Effects of Human Placenta, Bone Marrow Mesenchymal Stem Cells, and Their Exosomes for Damaged SVOG Ovarian Granulosa Cells. Stem Cells Int. 2020;2020:8861557.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Xiong Y, Liu T, Wang S, Chi H, Chen C, Zheng J. Cyclophosphamide promotes the proliferation inhibition of mouse ovarian granulosa cells and premature ovarian failure by activating the lncRNA-Meg3-p53-p66Shc pathway. Gene. 2017;596:1–8.

    Article  CAS  PubMed  Google Scholar 

  18. Yao G, He J, Kong Y, Zhai J, Xu Y, Yang G, et al. Transcriptional profiling of long noncoding RNAs and their target transcripts in ovarian cortical tissues from women with normal menstrual cycles and primary ovarian insufficiency. Mol Reprod Dev. 2019;86(7):847–61.

    Article  CAS  PubMed  Google Scholar 

  19. Zhao W, Dong L. Long non-coding RNA HOTAIR overexpression improves premature ovarian failure by upregulating Notch-1 expression. Exp Ther Med. 2018;16(6):4791–5.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. Bouckenheimer J, Assou S, Riquier S, Hou C, Philippe N, Sansac C, et al. Long non-coding RNAs in human early embryonic development and their potential in ART. Hum Reprod Update. 2016;23(1):19–40.

    Article  PubMed  Google Scholar 

  21. Tu J, Chen Y, Li Z, Yang H, Chen H, Yu Z. Long non-coding RNAs in ovarian granulosa cells. J Ovarian Res. 2020;13(1):63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Zheng Y, Bian Y, Wu R, Chen W, Fu L, Li P, et al. High-Throughput Sequencing Profiles About lncRNAs and mRNAs of Ovarian Granulosa Cells in Polycystic Ovary Syndrome. Front Med. 2021;8:741803.

    Article  Google Scholar 

  23. Wang L, Fan H, Zou Y, Yuan Q, Hu X, Chen X, et al. Aberrant Expression of Long Non-coding RNAs in Exosomes in Follicle Fluid From PCOS Patients. Front Genet. 2020;11:608178.

    Article  CAS  PubMed  Google Scholar 

  24. Khalaj K, Miller JE, Lingegowda H, Fazleabas AT, Young SL, Lessey BA, et al. Extracellular vesicles from endometriosis patients are characterized by a unique miRNA-lncRNA signature. JCI. Insight. 2019;4(18).

  25. Wu J, Huang H, Huang W, Wang L, Xia X, Fang X. Analysis of exosomal lncRNA, miRNA and mRNA expression profiles and ceRNA network construction in endometriosis. Epigenomics. 2020;12(14):1193–213.

    Article  CAS  PubMed  Google Scholar 

  26. Guo S, Quan S, Zou S. Roles of the Notch Signaling Pathway in Ovarian Functioning. Reprod Sci. 2021;28(10):2770–8.

    Article  CAS  PubMed  Google Scholar 

  27. Zheng C, Liu S, Qin Z, Zhang X, Song Y. LncRNA DLEU1 is overexpressed in premature ovarian failure and sponges miR-146b-5p to increase granulosa cell apoptosis. J Ovarian Res. 2021;14(1):151.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Ruszkowska M, Nynca A, Paukszto L, Sadowska A, Swigonska S, Orlowska K, et al. Identification and characterization of long non-coding RNAs in porcine granulosa cells exposed to 2,3,7,8-tetrachlorodibenzo-p-dioxin. J Anim Sci Biotechnol. 2018;9:72.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Majewska M, Lipka A, Paukszto L, Jastrzebski JP, Gowkielewicz M, Jozwik M, et al. Preliminary RNA-Seq Analysis of Long Non-Coding RNAs Expressed in Human Term Placenta. Int J Mol Sci. 2018;19(7).

  30. Dang Y. The role and mechanism of non-coding RNA in primary ovarian insufficiency [ Ph.D. thesis]: Shandong University; 2016.

    Google Scholar 

  31. Li D, Wang X, Dang Y, Zhang X, Zhao S, Lu G, et al. lncRNA GCAT1 is involved in premature ovarian insufficiency by regulating p27 translation in GCs via competitive binding to PTBP1. Mol Ther Nucleic Acids. 2021;23:132–41.

    Article  CAS  PubMed  Google Scholar 

  32. Li D, Xu W, Wang X, Dang Y, Xu L, Lu G, et al. lncRNA DDGC participates in premature ovarian insufficiency through regulating RAD51 and WT1. Mol Ther Nucleic Acids. 2021;26:1092–106.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Wang F, Chen X, Sun B, Ma Y, Niu W, Zhai J, et al. Hypermethylation-mediated downregulation of lncRNA PVT1 promotes granulosa cell apoptosis in premature ovarian insufficiency via interacting with Foxo3a. J Cell Physiol. 2021;236(7):5162–75.

    Article  CAS  PubMed  Google Scholar 

  34. Harvey CN, Esmail M, Wang Q, Brooks AI, Zachow R, Uzumcu M. Effect of the methoxychlor metabolite HPTE on the rat ovarian granulosa cell transcriptome in vitro. Toxicol Sci. 2009;110(1):95–106.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Wang X, Meng K, Wang H, Wang Y, Zhao Y, Kang J, et al. Identification of small extracellular vesicle subtypes in follicular fluid: Insights into the function and miRNA profiles. J Cell Physiol. 2021;236(8):5633–45.

    Article  CAS  PubMed  Google Scholar 

  36. Chen Y, Chang H, Cheng J, Tsai H, Wu C, Leung P. Transforming growth factor-β1 up-regulates connexin43 expression in human granulosa cells. Hum Reprod (Oxford, England). 2015;30(9):2190–201.

    Article  CAS  Google Scholar 

  37. Craig J, Orisaka M, Wang H, Orisaka S, Thompson W, Zhu C, et al. Gonadotropin and intra-ovarian signals regulating follicle development and atresia: the delicate balance between life and death. Front Biosci. 2007;12:3628–39.

    Article  CAS  PubMed  Google Scholar 

  38. Wang Y, Hu SB, Wang MR, Yao RW, Wu D, Yang L, et al. Genome-wide screening of NEAT1 regulators reveals cross-regulation between paraspeckles and mitochondria. Nat Cell Biol. 2018;20(10):1145–58.

    Article  CAS  PubMed  Google Scholar 

  39. Bouckenheimer J, Fauque P, Lecellier CH, Bruno C, Commes T, Lemaître JM, et al. Differential long non-coding RNA expression profiles in human oocytes and cumulus cells. Sci Rep. 2018;8(1):2202.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Ernst EH, Nielsen J, Ipsen MB, Villesen P, Lykke-Hartmann K. Transcriptome Analysis of Long Non-coding RNAs and Genes Encoding Paraspeckle Proteins During Human Ovarian Follicle Development. Front Cell Dev Biol. 2018;6:78.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Zhen J, Li J, Li X, Wang X, Xiao Y, Sun Z, et al. Downregulating lncRNA NEAT1 induces proliferation and represses apoptosis of ovarian granulosa cells in polycystic ovary syndrome via microRNA-381/IGF1 axis. J Biomed Sci. 2021;28(1):53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Li M, Peng J, Zeng Z. Overexpression of long non-coding RNA nuclear enriched abundant transcript 1 inhibits the expression of p53 and improves premature ovarian failure. Exp Ther Med. 2020;20(5):69.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Mahboobeh Z, Pegah M, Fatemeh S, Elham K, Hanieh A, Milad R, et al. lncRNA ZEB2-AS1: A promising biomarker in human cancers. IUBMB Life. 2020;72(9):1891–9.

    Article  CAS  PubMed  Google Scholar 

  44. Zhang B, Fan DB, Liu L, Qin Y, Feng DQ. Knockdown of ZEB2-AS1 inhibits cell proliferation, invasion and induces apoptosis in osteosarcoma by combining with EZH2. Eur Rev Med Pharmacol Sci. 2020;24(12):6533–9.

    CAS  PubMed  Google Scholar 

  45. Xiang Z, Lv Q, Chen X, Zhu X, Liu S, Li D, et al. Lnc GNG12-AS1 knockdown suppresses glioma progression through the AKT/GSK-3β/β-catenin pathway. Biosci Rep. 2020;40(8).

  46. Yu N, Chen X, Du M, Li H, Wang Y, Jiang F, et al. Long Non-Coding RNA ZEB2-AS1 Augments Activity of Trophoblast Cells and Prevents the Development of Recurrent Spontaneous Abortion in Mice Through EZH2-Mediated CST3 Inhibition. Reprod Sci. 2022;29(3):963–74.

    Article  CAS  PubMed  Google Scholar 

  47. Patil K, Hinduja I, Mukherjee S. Alteration in angiogenic potential of granulosa-lutein cells and follicular fluid contributes to luteal defects in polycystic ovary syndrome. Hum Reprod. 2021;36(4):1052–64.

    Article  CAS  PubMed  Google Scholar 

  48. Castillo-Higuera T, Alarcón-Granados MC, Marin-Suarez J, Moreno-Ortiz H, Esteban-Pérez CI, Ferrebuz-Cardozo AJ, et al. A Comprehensive Overview of Common Polymorphic Variants in Genes Related to Polycystic Ovary Syndrome. Reprod Sci. 2021;28(9):2399–412.

    Article  CAS  PubMed  Google Scholar 

  49. Biase FH. Oocyte Developmental Competence: Insights from Cross-Species Differential Gene Expression and Human Oocyte-Specific Functional Gene Networks. Omics. 2017;21(3):156–68.

    Article  CAS  PubMed  Google Scholar 

  50. Santos Simões R, Carbonel AAF, Borges FT, Baracat MCP, da Silva Sasso GR, Simões MJ, et al. Analysis of hyaluronic acid in the endometrium of women with polycystic ovary syndrome. Gynecol Endocrinol. 2019;35(2):133–7.

    Article  PubMed  Google Scholar 

  51. Zhang GM, Gu CH, Zhang YL, Sun HY, Qian WP, Zhou ZR, et al. Age-associated changes in gene expression of goat oocytes. Theriogenology. 2013;80(4):328–36.

    Article  CAS  PubMed  Google Scholar 

  52. Guda MR, Velpula KK, Asuthkar S, Cain CP, Tsung AJ. Targeting RGS4 Ablates Glioblastoma Proliferation. Int J Mol Sci. 2020;21(9).

  53. Keren-Raifman T, Bera AK, Zveig D, Peleg S, Witherow DS, Slepak VZ, et al. Expression levels of RGS7 and RGS4 proteins determine the mode of regulation of the G protein-activated K(+) channel and control regulation of RGS7 by G beta 5. FEBS Lett. 2001;492(1–2):20–8.

    Article  CAS  PubMed  Google Scholar 

  54. Guerrero-Netro HM, Chorfi Y, Price CA. Effects of the mycotoxin deoxynivalenol on steroidogenesis and apoptosis in granulosa cells. Reproduction. 2015;149(6):555–61.

    Article  CAS  PubMed  Google Scholar 

  55. Wang H, He Y, Cheng D, Pu D, Tan R, Gao L, et al. Cypermethrin exposure reduces the ovarian reserve by causing mitochondrial dysfunction in granulosa cells. Toxicol Appl Pharmacol. 2019;379:114693.

    Article  CAS  PubMed  Google Scholar 

  56. Ożegowska K, Brązert M, Ciesiółka S, Nawrocki MJ, Kranc W, Celichowski P, et al. Genes Involved in the Processes of Cell Proliferation, Migration, Adhesion, and Tissue Development as New Potential Markers of Porcine Granulosa Cellular Processes In Vitro: A Microarray Approach. DNA Cell Biol. 2019;38(6):549–60.

    Article  PubMed  Google Scholar 

  57. Kranc W, Brązert M, Budna J, Celichowski P, Bryja A, Nawrocki MJ, et al. Genes responsible for proliferation, differentiation, and junction adhesion are significantly up-regulated in human ovarian granulosa cells during a long-term primary in vitro culture. Histochem Cell Biol. 2019;151(2):125–43.

    Article  CAS  PubMed  Google Scholar 

  58. Pashaiasl M, Ebrahimi M, Ebrahimie E. Identification of the key regulating genes of diminished ovarian reserve (DOR) by network and gene ontology analysis. Mol Biol Rep. 2016;43(9):923–37.

    Article  CAS  PubMed  Google Scholar 

  59. Esmaeili-Fard SM, Gholizadeh M, Hafezian SH, Abdollahi-Arpanahi R. Genome-wide association study and pathway analysis identify NTRK2 as a novel candidate gene for litter size in sheep. PLoS One. 2021;16(1):e0244408.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Qazi IH, Cao Y, Yang H, Angel C, Pan B, Zhou G, et al. Impact of Dietary Selenium on Modulation of Expression of Several Non-Selenoprotein Genes Related to Key Ovarian Functions, Female Fertility, and Proteostasis: a Transcriptome-Based Analysis of the Aging Mice Ovaries. Biol Trace Elem Res. 2021;199(2):633–48.

    Article  PubMed  Google Scholar 

  61. Yang X, Zhou Y, Peng S, Wu L, Lin HY, Wang S, et al. Differentially expressed plasma microRNAs in premature ovarian failure patients and the potential regulatory function of mir-23a in granulosa cell apoptosis. Reproduction. 2012;144(2):235–44.

    Article  CAS  PubMed  Google Scholar 

  62. Kirshenbaum M, Orvieto R. Premature ovarian insufficiency (POI) and autoimmunity-an update appraisal. J Assist Reprod Genet. 2019;36(11):2207–15.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Sharif K, Watad A, Bridgewood C, Kanduc D, Amital H, Shoenfeld Y. Insights into the autoimmune aspect of premature ovarian insufficiency. Best Pract Res Clin Endocrinol Metab. 2019;33(6):101323.

    Article  CAS  PubMed  Google Scholar 

  64. Lim W, Bae H, Bazer FW, Song G. Cell-specific expression and signal transduction of C-C motif chemokine ligand 2 and atypical chemokine receptors in the porcine endometrium during early pregnancy. Dev Comp Immunol. 2018;81:312–23.

    Article  CAS  PubMed  Google Scholar 

  65. Matsuda F, Inoue N, Manabe N, Ohkura S. Follicular growth and atresia in mammalian ovaries: regulation by survival and death of granulosa cells. J Reprod Dev. 2012;58(1):44–50.

    Article  CAS  PubMed  Google Scholar 

  66. Fan Y, Chang Y, Wei L, Chen J, Li J, Goldsmith S, et al. Apoptosis of mural granulosa cells is increased in women with diminished ovarian reserve. J Assist Reprod Genet. 2019;36(6):1225–35.

    Article  PubMed  PubMed Central  Google Scholar 

  67. Hall SE, Upton RMO, McLaughlin EA, Sutherland JM. Phosphoinositide 3-kinase/protein kinase B (PI3K/AKT) and Janus kinase/signal transducer and activator of transcription (JAK/STAT) follicular signalling is conserved in the mare ovary. Reprod Fertil Dev. 2018;30(4):624–33.

    Article  CAS  PubMed  Google Scholar 

  68. Frost ER, Ford EA, Peters AE, Reed NL, McLaughlin EA, Baker MA, et al. Signal transducer and activator of transcription (STAT) 1 and STAT3 are expressed in the human ovary and have Janus kinase 1-independent functions in the COV434 human granulosa cell line. Reprod Fertil Dev. 2020;32(12):1027–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Ding N, Harlow SD, Randolph JF Jr, Loch-Caruso R, Park SK. Perfluoroalkyl and polyfluoroalkyl substances (PFAS) and their effects on the ovary. Hum Reprod Update. 2020;26(5):724–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Tatone C, Benedetti E, Vitti M, Di Emidio G, Ciriminna R, Vento ME, et al. Modulating Intrafollicular Hormonal Milieu in Controlled Ovarian Stimulation: Insights From PPAR Expression in Human Granulosa Cells. J Cell Physiol. 2016;231(4):908–14.

    Article  CAS  PubMed  Google Scholar 

  71. Osuka S, Iwase A, Goto M, Takikawa S, Nakamura T, Murase T, et al. Thyroid Autoantibodies do not Impair the Ovarian Reserve in Euthyroid Infertile Women: A Cross-Sectional Study. Horm Metab Res. 2018;50(7):537–42.

    Article  CAS  PubMed  Google Scholar 

  72. Zhang Y, Ouyang X, You S, Zou H, Shao X, Zhang G, et al. Effect of human amniotic epithelial cells on ovarian function, fertility and ovarian reserve in primary ovarian insufficiency rats and analysis of underlying mechanisms by mRNA sequencing. Am J Transl Res. 2020;12(7):3234–54.

    CAS  PubMed  PubMed Central  Google Scholar 

  73. Ferraretti AP, La Marca A, Fauser BCJM, Tarlatzis B, Nargund G, Gianaroli L. ESHRE consensus on the definition of 'poor response' to ovarian stimulation for in vitro fertilization: the Bologna criteria. Human Reproduction. 2011;26(7):1616–24.

    Article  CAS  PubMed  Google Scholar 

  74. Kechin A, Boyarskikh U, Kel A, Filipenko M. cutPrimers: A New Tool for Accurate Cutting of Primers from Reads of Targeted Next Generation Sequencing. J Comput Biol. 2017;24(11):1138–43.

    Article  CAS  PubMed  Google Scholar 

  75. Beekman R, Chapaprieta V, Russiñol N, Vilarrasa-Blasi R, Verdaguer-Dot N, Martens JHA, et al. The reference epigenome and regulatory chromatin landscape of chronic lymphocytic leukemia. Nat Med. 2018;24(6):868–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35(Web Server issue):W345–9.

    Article  PubMed  PubMed Central  Google Scholar 

  79. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41(17):e166.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  CAS  PubMed  Google Scholar 

  81. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.

    Article  PubMed  PubMed Central  Google Scholar 

  82. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36(Database issue):D480–4.

    CAS  PubMed  Google Scholar 

  83. Wang YK, Li X, Song ZQ, Yang CX. Methods of RNA preparation affect mRNA abundance quantification of reference genes in pig maturing oocytes. Reprod Domest Anim. 2017;52(5):722–30.

    Article  CAS  PubMed  Google Scholar 

  84. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(−Delta Delta C(T)) Method. Methods. 2001;25(4):402–8.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgments

The authors thank the nurses and laboratory staff of the Department of Assisted Reproduction of the Affiliated Hospital of Shandong University of Traditional Chinese Medicine for their technical assistance.

Funding

This research was funded by Natural Science Foundation of Shandong Province (NO. ZR2021MH255), National Natural Science Foundation of China (NO. 81703958, 81774355, 81974577 and 82274573) and Key Technology Research and Development Program of Shandong (NO. 2019GSF108067).

Author information

Authors and Affiliations

Authors

Contributions

FL and HW conceived and designed this study. LD contributed to the acquisition and analyses of data. LD and XX contributed to drafting the manuscript. CY was responsible for sample collection. HMC and PCKL contributed to manuscript revision. All authors read and approved the final submitted version of this manuscript.

Corresponding authors

Correspondence to Fang Lian or Haicui Wu.

Ethics declarations

Ethics approval and consent to participate

The study was conducted in accordance with the Declaration of Helsinki, and approved by the Ethics Committee of the Reproductive Medicine Ethics Committee of the Affiliated Hospital of the Shandong University of Traditional Chinese Medicine (Ref. SZ2018090103). Written informed consent was obtained from all the participants.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no commercial or financial relationships that could be construed as a potential conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1 Supplemental Table 1

Alignment of the sequencing reads to the reference genome. Supplemental Table 2 Statistical analyses of lncRNAs in each sample. Supplemental Table 3 Statistical analyses of the DE lncRNAs. Supplemental Table 4 Statistical analyses of the DE lncRNA cis-regulated target genes. Supplemental Table 5 GO terms enriched for the DE mRNA genes. Supplemental Table 6 The lncRNA–mRNA coexpression pairs. Supplemental Table 7 Specific primer sequences of the validated mRNAs and lncRNAs.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Dong, L., Xin, X., Chang, HM. et al. Expression of long noncoding RNAs in the ovarian granulosa cells of women with diminished ovarian reserve using high-throughput sequencing. J Ovarian Res 15, 119 (2022). https://doi.org/10.1186/s13048-022-01053-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13048-022-01053-6

Keywords