An Integrated Analysis of circRNA and lncRNA Expression of Bovine Granulosa Cells Induced by Melatonin Reveals the Pathways Potentially Involved in Follicular Development

TL;DR

This study explores how melatonin affects bovine granulosa cells by analyzing changes in lncRNA and circRNA expression, revealing pathways involved in follicular development.

Contribution

The study identifies novel lncRNAs and circRNAs regulated by melatonin and their potential roles in granulosa cell function and follicular development.

Findings

Melatonin treatment altered the expression of 99 lncRNAs and 28 circRNAs in bovine granulosa cells.

Differentially expressed RNAs were enriched in pathways related to development, apoptosis, and reproductive functions.

An lncRNA/circRNA and miRNA regulatory network was constructed to explore their interactions in follicular development.

Abstract

Objective: Accumulating evidence demonstrates that melatonin is involved in modulating granulosa cell function and follicular development. lncRNAs (long non-coding RNAs) and circRNAs (circular RNAs) have been reported to participate in multiple biological processes. This study aimed to explore the candidate circRNAs and lncRNAs related to molecular mechanisms when exploring the role of melatonin in regulating ovarian function. Methods: Bovine ovary granulosa cells were collected 48 h after treatment with melatonin at 10−7 M. The lncRNA and circRNA profiles of bovine granulosa cells were further explored using high-throughput sequencing in the absence/presence of melatonin. The differentially expressed lncRNAs and circRNAs were analyzed through the annotation information of source transcripts for GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes). Results: We identified 99 differentially expressed lncRNAs and 28 differentially expressed circRNAs. Enrichment analysis of differentially expressed lncRNAs and circRNAs showed they were enriched in multiple pathways involved in development, apoptosis, and reproductive function, such as the mTOR (mammalian Target of Rapamycin) signaling pathway, FoxO (Forkhead box O) signaling pathway, MAPK (Mitogen-Activated Protein Kinase) signaling pathway, Hippo signaling pathway, TGF-beta (Transforming Growth Factor-β) signaling pathway, PI3K-Akt (Phosphatidylinositol 3-Kinase-Akt) signaling pathway, apoptosis, and Rap1 (Ras-related protein 1), most of which were mainly related to granulosa cell function and the crosstalk between granulosa cells and oocytes. The present analysis indicated the potential role of melatonin in granulosa cell function by regulating lncRNA and circRNA expression and, thus, mediating follicular development. An lncRNA/circRNA and miRNA regulatory network was also constructed to take their interactions into account. Conclusions: Our study offers details of lncRNA and circRNA expression in bovine granulosa cells and further provides insight into the potential role of melatonin in regulating reproduction by modulating lncRNA and circRNA expression.

Figures

Full text

1. Introduction

Follicular development is accompanied by a considerable amount of follicular atresia, around 99%, which ensures that the healthy follicles containing optimal-quality oocytes can reach the ovulatory stage [1]. Many factors are involved in follicular atresia, and granulosa cells are a major contributing factor affecting follicular development. The main reason for inducing antral follicular atresia is the apoptosis of granulosa cells [2]. Granulosa cell apoptosis can be initiated by inflammation and oxidative stress [3]. Numerous factors participate in the contribution of granulosa cells to the regulation of follicular development, such as steroid hormones and transforming growth factor beta superfamily derived from granulosa cells [4]. It is crucial to better understand the importance of interactions between granulosa cells and oocytes for oocyte maturation, which is beneficial for improving fertility. However, the regulatory mechanism of granulosa cells for promoting follicular development is still not fully characterized.

CircRNAs have been identified as non-coding RNA and come from mRNA processing before splicing. Accumulating evidence has revealed that circRNAs play a potentially considerable role in biological development [5]. Due to circRNAs’ function in cell processes, considerable interest has increasingly led to more research. The main reason is that circRNAs could reduce the miRNA-inhibiting effect on target genes by blocking the miRNA from binding to target genes [6]. Moreover, circRNAs can modulate the intracellular transport of RNA-binding proteins and the stability of mRNAs and then influence biological processes, for example, through circFoxo3 and circRNA cerebellar degeneration protein 1. circFoxo3 could arrest the cell cycle in the G1/S phase by binding to cyclin-dependent protein kinase 2 [7], and circRNA cerebellar degeneration protein 1 could affect the stability of mRNAs [8]. Furthermore, HuR, another widely studied RNA-binding protein, combined with circPABPN1(poly(A) binding protein nuclear 1) could inhibit PABPN1 translation [9], which means that circRNA and mRNA could also competitively combine and regulate the function of RNA-binding proteins. As regards the role of circRNA in reproductive function, increasing amounts of research indicate that circRNAs are involved in regulating ovarian follicle development. circRNA expression exhibits different characteristics in human granulosa cells from 30 and 38 years and in porcine granulosa cells derived from healthy and atretic antral follicles [10,11]. Recent research provides evidence for the circINHA/miR-10a-5p/CTGF (connective tissue growth factor) regulatory pathway in pig granulosa cell apoptosis and provides novel insights into the role of circRNAs in the modulation of ovarian physiological functions [12]. LncRNAs are also drawing increased attention as another important non-coding RNA, which also perform important roles in many biological processes by binding miRNA, such as embryonic development, oocyte meiotic resumption, and the progesterone-mediated oocyte maturation cell cycle [13,14,15]. Increasing evidence supports that lncRNAs are involved in the regulation of seasonal reproduction [16]. Therefore, lncRNAs and circRNAs may participate in their special roles in modulating the ovarian function and improving fertility through affecting granulosa cell function.

Communication between granulosa cells and oocytes is complex and crucial for oocyte quality and follicular development. An important product secreted by granulosa cells is melatonin, and accumulating researches suggest that melatonin is crucial for modulating reproductive function and follicular development, such as improving oocyte maturation, embryo development and inducing hormone secretion [2,17]. In addition, the higher concentration in follicular fluid is consistent with the role of melatonin in regulating ovarian function [18]. Previous studies showed that melatonin was involved in modulating the antioxidant activity, antiapoptosis and the steroid hormone secretion of granulosa cells, as well as promoting oocyte maturation and embryo development in bovines [19,20,21].

However, the lncRNA and circRNA expression profiles in the granulosa cells of bovines are still unclear. Specifically, melatonin regulates the lncRNA and circRNA expression characteristics in bovine granulosa cells. We investigated the features of lncRNAs and circRNAs in bovine granulosa cell treatment with melatonin compared to that without treatment using RNA sequencing. Furthermore, the potential functions of differentially expressed lncRNAs and circRNAs were analyzed using GO annotation and KEGG enrichment. We also constructed a regulatory network of lncRNAs/circRNAs and miRNAs to take into account their interactions. These results might detail new insights for bovine granulosa cell lncRNAs and circRNAs and provide a better understanding of the potential roles of lncRNAs and circRNAs in regulating bovine follicular development.

2. Materials and Methods

2.1. Bovine Granulosa Cell Culture and Cell Treatment

Granulosa cell isolation and collection were operated following our previously described work, with minor revisions [19,20]. Bovine ovaries were obtained from a local abattoir (Bengbu, China) and sent back to the laboratory within 3 h in a thermos cup. About 80 bovine ovaries were collected by washing three times using 70% alcohol, and then the ovaries were washed three times using sterile 0.9% NaCl to remove alcohol. The follicular fluids obtained from the follicle (diameter 5–8 mm) were centrifugated at 1500 rpm for 5 min. The collected cell pellets were digested by 0.25% trypsin with 0.025% EDTA (Gibco, Grand Island, NY, USA) for 5 min. After being digested, the cell pellets were centrifugated again and dispersed in Dulbecco’s Modified Eagle Medium (DMEM) (Gibco, Grand Island, NY, USA) supplemented with 10% fetal bovine serum (FBS; Hyclone, UT, USA) and antibiotics, including streptomycin (50 µg/mL), penicillin (50 IU/mL) (Pen-Strep, Invitrogen, Carlsbad, CA, USA) and plasmocin (25 µg/mL; invivogen, San Diego, CA, USA). Finally, the separated cells were seeded into the 60 mm cell culture dishes. Melatonin (Sigma, St. Louis, MO, USA) was diluted in DMSO (Dimethyl Sulfoxide, Sigma, St. Louis, MO, USA) to a concentration of 0.01 M. The final concentration of DMSO was adjusted to 0.2% in all dilutions, and 0.2% DMSO was set as the control. Then, granulosa cells were distributed into two groups as follows: the control group (the untreated cells incubated with culture medium containing 0.2% DMSO) and the melatonin group (the cells exposed to 10^−7^ M melatonin treatment). One day before treatment, 2–5 × 10^5^ cells were cultured in a 12-well plate to reach 70–80% confluence at the time of treatment, and the medium was replaced with fresh medium containing melatonin. The granulosa cells were cultured in an incubator at 37 °C with 5% CO_2_ and harvested 48 h after treatment. Three biological repeats were performed for all experiments. The granulosa cells were immediately frozen in liquid nitrogen after harvesting, followed by storage at −80 °C, for subsequent RNA-seq and qRT-PCR analysis. In this study, the experimental protocols were reviewed and approved by the Anhui Science and Technology University Institutional Committee on Animal Care and Use.

2.2. Total RNA Extraction

Total RNA from granulosa cells treated with 10^−7^ M melatonin for 48 h was extracted using Trizol regent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol and digested with RNase-Free DNase to remove residual DNAs. The purity and concentration of total RNA were assessed using a Nano Drop 2000 (Thermo, Waltham, CA, USA). The RNA integrity was evaluated using a QubitRNA Assay Kit (Life Technologies, Carlsbad, CA, USA) with the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The samples with RNA Integrity Number (RIN) ≥ 7 were subjected to subsequent analysis.

2.3. cDNA Library Construction and RNA Sequencing

The libraries were constructed using TruSeq Stranded Total RNA with Ribo-Zero Gold according to the manufacturer’s instructions (Illumina, San Diego, CA, USA). Briefly, ribosomal RNA was depleted from the total RNA using Ribo-Zero Gold. The rRNA-depleted RNA was broken at random by adding fragmentation buffer. The randomly fragmented RNA was adapted to synthesize the first-strand cDNA using the SuperScript II Reverse Transcriptase (Invitrogen, San Diego, CA, USA). Subsequently, the first-strand cDNA was used as the template strand, and we manufactured the second cDNA strand using the DNA polymerase I system. After purification, A-tailing and adaptors were tagged into the double-stranded cDNA. After adenylating 3′ Ends and ligating adapters, the cDNA was purified using AMPure XP beads (Beckman Coulter, Indianapolis, IN, USA). The suitable fragments (150–200 bp in length) were chosen as templates for PCR amplification. Finally, DNA fragments were enriched by PCR amplification using primer cocktail, and then the library was constructed. At last, these libraries were sequenced on the Illumina HiSeqTM 2500 sequencing platform (Illumina, San Diego, CA, USA) with paired-end 150 bp reads. High-throughput sequencing was carried out by OE Biotech Co., Ltd. (Shanghai, China).

2.4. RNA-Seq Data Preprocessing

Briefly, raw reads were generated from high-throughput sequencing platform Illumina HiSeqTM 2500 sequencer. In order to obtain high-quality reads that could be used for later analysis, raw reads needed to be further quality filtered. Adapter sequences, reads containing poly-N (N > 10%), and N-bases or low-quality reads were filtered out using Trimmomatic software (v0.39). Finally, the remaining high-quality clean reads were obtained to analyze circRNAs and lncRNA. Clean reads were mapped to the bovine reference genome (Bos taurus ARS-UCD1.2) by hisat2 (v2.1.0). BWA software (v0.7.17) was used to align the sequencing reads of each sample with reference genome and generate SAM file. Then, PCC signals (paired chiastic clipping signals) were scanned using CIRI2 software (v2.0.6) from the SAM file, and circRNA sequences were predicted based on junction reads and GT-AG cleavage signals.

For lncRNAs, the result of alignment with the reference genome was stored in a binary file, called a bam file. Using the Stringtie software (v2.2.1) to assemble the reads, the new transcript was spliced. Then, candidate lncRNA transcripts were selected by comparing the gene annotation information of the reference sequence produced by Cuffcompare software (v2.2.1). Finally, new transcripts with coding potential were screened out by CPC (Coding Potential Calculator), CNCI (Coding-Non-Coding-Index), Pfam and PLEK to obtain lncRNA predicted sequences. Aligning the sequencing reads of each sample with the sequence of known lncRNA sequences and lncRNA prediction sequences by bowtie2 (v2.5.1), and using eXpress (v1.5.1) to make gene quantitative analysis, the FPKM (fragments per kb per million reads) value and counts value (the number of reads for each gene in each sample) were obtained. The characteristics of obtained lncRNAs were then analyzed, including the type, length distribution and number of exons (lncRNA exons do not encode functional proteins).

The expression levels of circRNAs and lncRNAs were measured from the read counts and were normalized by the RPM (reads per million mapped reads) and FPKM, respectively. The raw read counts were normalized using the estimate Size Factors function of the DESeq R package (v1.48.0), and the nbinomTest function was subsequently employed to calculate the p-values and fold-change values for differential expression analysis. Foldchange values were carried out for melatonin (melatonin1, melatonin2, melatonin3) and control (control1, control2, control3) groups using DESeq. CircRNAs and lncRNAs were considered to be differentially expressed between the samples with a fold change > 2 and p < 0.05 (not FDR-corrected). There were three biological and technical replicates used to achieve the claimed statistical power.

2.5. Bioinformatics Analysis

Enrichment of differentially expressed circRNAs and lncRNAs was performed through the annotation information of source transcripts for GO and KEGG analysis. The GO term analysis could provide the differentially expressed circRNAs and lncRNAs involved in biological functions, including biological processes, molecular functions and cellular components, using the database (http://geneontology.org/, accessed on 15 March 2025). The differentially expressed circRNAs and lncRNAs taking part in signaling pathways were identified using the KEGG database (https://www.kegg.jp/kegg/, accessed on 15 March 2025). The cluster profiles R package (v4.6.2) and KOBAS software (v2.0) were used to perform the functional annotation of circRNAs and lncRNAs. p < 0.05 (not FDR-corrected) was considered as significant for functional enrichment. Finally, circRNA/lncRNA and miRNA interactions were predicted with miRanda software (v3.3a).

2.6. Gene Expression Validation via Quantitative PCR (qPCR)

Four differentially expressed circRNAs and lncRNAs were selected and identified using qPCR. Approximately 1 μg of total RNA was reverse transcribed into cDNA using a RevertAid First Strand cDNA Synthesis Kit k1622 (Thermo Scientific, Waltham, MA, USA). Reverse transcription was carried out in 20 μL reaction volume. Briefly, 1 μg RNA templete, 1 μL Oligo(dT)18 primer, 4 μL 5×Reaction Buffer, 1 μL RiboLock RNase Inhibitor, 2 μL 10 mM dNTP Mix, 1 μL RevertAid M-MuLV RT, and nuclease-free water were mixed well, and the reaction mixtures used the following conditions: 42 °C for 60 min, 70 °C for 5 min. The cDNA was stored at −20 °C.

qPCR was performed using LightCycler 480 SYBR Green I Master Mix on a LightCycler 480II Real-Time PCR System (Roche, Penzberg, Germany). A total of 10 μL reaction volume was prepared containing 5 μL LightCycler 480 SYBR Green I Master Mix, 1 μL reverse-transcribed cDNA, 0.5 μM specific primer, and 3 μL RNase and DNase-free water. The qPCR was performed as follows: 95 °C for 5 min, followed by 45 cycles of denaturation at 95 °C for 20 s, annealing at particular temperatures for 20 s, extension at 72 °C for 20 s. Melting curve was established from 65 °C to 95 °C to identify specific PCR product. All primers designed for circRNAs and lncRNAs are shown in Table S1. All samples were amplified in triplicate, and the expression of identified circRNAs and lncRNAs was normalized to β-actin used as a control. The relative expression values of circRNAs and lncRNAs were analyzed using the 2^−ΔΔCt^ method.

2.7. Statistical Analysis

The qPCR data were displayed as the mean ± standard deviations (SD) of three replicates. Significance differences between melatonin and control groups were determined using one-way ANOVA with SPSS 17. The significant difference is presented as p < 0.05.

3. Results

3.1. Summary of Raw Sequence Reads

After low-quality bases and N-bases or low-quality reads were removed, a total of 309,386,934 and 307,886,210 clean reads were obtained in control and melatonin with more than 90.25% of Q30, respectively (Table 1). Greater than 95% of the reads were successfully coordinated to the Bos taurus reference genome (ARS-UCD1.2).

3.2. Characterization of circRNAs and lncRNAs Expressed in Bovine Granulosa Cells Treated with Melatonin

A total of 3443 circRNAs were identified in the melatonin and control groups (Supplementary Table S1). After annotation, circRNAs were divided into five subclasses, denoted as sense-overlapping (3032), exonic (166), intergenic (110), intronic (57) and antisense (78). Most circRNAs were produced from sense overlapping according to the genomic location and features (Figure 1A). The majority of the circRNAs contained more than one exon, and the peaks of exon numbers were 1–6 (Figure 1D). Most host genes (1159/1873) only generated a single circRNA, while some (714/1873) produced multiple circRNAs. Further analysis of the identified circRNAs indicated that they were widely distributed among the chromosomes, and chromosome 1 had the highest number of circRNAs, followed by chromosome 2 and 3 (Figure 1F). The length of circRNAs was mostly distributed on 200–700 and more than 2000 nucleotides (Figure 1H).

Additionally, the length of transcripts > 200 bp and exon count ≥ 1 were used to screen out the lncRNAs through CPC, CNCI, Pfam and PLEK (Figure 1B, Supplementary Table S2). As a result, 6767 lncRNA transcripts (6503 known and 264 novel) were found (Supplementary Table S3). Among these, lncRNAs with lengths greater than 2000 bp accounted for the highest proportion (Figure 1I), and most lncRNAs contained two to six exons (Figure 1E). The identified lncRNAs in granulosa cells were distributed across 29 autosomes and X-chromosomes, with the most on chromosome 7 (322 lncRNAs) and the fewest on chromosome 28 (88 lncRNAs). Moreover, three lncRNAs were located on chromosome MT (Figure 1G). According to the genomic loci, lncRNAs were classified as genic exonic, genic intronic, intergenic downstream and intergenic upstream. Most (1136) lncRNAs originated from antisense genic intronic regions (Figure 1C).

3.3. Identification of Differentially Expressed circRNAs and lncRNAs

The differentially expressed circRNAs and lncRNAs were analyzed according to the RNA sequencing, and the significant differences in these circRNAs and lncRNAs were taken into account based on the parameter of log2 fold change > 1 and p < 0.05 (not FDR-corrected). The average RPM and FPKM values were calculated for each circRNA and lncRNA in the granulosa cells from the melatonin and control groups, respectively (Figure 2A,B). The differential expression analysis exhibited that there were 28 circRNAs and 99 lncRNAs in the bovine granulosa cells before and after treatment with melatonin, among which 13 circRNAs and 58 lncRNAs were significantly upregulated and 15 circRNAs and 41 lncRNAs were significantly downregulated (Figure 2C,D, Supplementary Tables S4 and S5; the differential expression results did not remain significant after FDR correction). These differentially expressed circRNAs and lncRNAs, some downregulated and others upregulated, were shown in volcanic plots (Figure 2E,F). A hierarchical clustering heatmap was used to visualize and cluster the differentially expressed and similar expression pattern of circRNAs and lncRNAs between the melatonin and control groups (Figure 2G,H).

3.4. GO and KEGG Analysis of Differentially Expressed circRNAs and lncRNAs

CircRNAs are produced from the host genes by back-splicing. Therefore, analyzing GO classification terms of circRNA host genes can offer new insights. To explore these insights, circRNAs and lncRNAs associated with biological processes, molecular functions, and cellular components were analyzed for GO term enrichment (Figure 3). After melatonin treatment, we found that circRNAs and lncRNAs participated in the biological processes related to ovarian function, including biological regulation, the developmental process, metabolic process, reproduction, reproductive process and signaling. For the cellular components, the important categories were cell, cell part and cell junction, which may be beneficial to the communication between granulosa cells and oocytes. As regards the molecular functions, the dominant terms were binding, catalytic activity and enzyme regulator activity. Interestingly, receptor regulator activity was also involved in the cellular components, meaning melatonin regulates the bovine granulosa cell function. This finding further confirms the regulatory network between melatonin and granulosa cell function.

The KEGG database was used to explore the potential roles of differentially expressed circRNAs and lncRNAs in the bovine granulosa cells. After melatonin treatment, several pathways were significantly enriched (Table 2), including the mTOR signaling pathway, tight junction, MAPK signaling pathway, Hippo signaling pathway, apoptosis and Rap1. Most of these pathways were mainly related to the granulosa cell function and the crosstalk between granulosa cells and oocytes. Additionally, the differentially expressed lncRNAs involved in the pathways were also analyzed (Figure 4, Supplementary Table S6). Notably, pathways, such as apoptosis-multiple species, mTOR signaling pathway, insulin signaling pathway, Gap junction, pI3K-Akt signaling pathway, TGF-beta signaling pathway, HIF-1 signaling pathway and FoxO signaling pathway, were enriched in the melatonin groups, which are related to follicular development, interactions, reproduction and the endocrine system. Therefore, the differentially expressed circRNAs and lncRNAs generated by the melatonin treatment probably play an important role in granulosa cells, thereby influencing oocyte maturation and follicular development.

3.5. Predictive Analysis of Target miRNAs of Differentially Expressed circRNAs and lncRNAs

The miRNA targets of the circRNAs and lncRNAs were predicted using bioinformatics tools. The majority of the circRNAs and lncRNAs had miRNA binding sites (Supplementary Tables S7 and S8). Moreover, the predictive analysis results showed that (109 and 165) miRNAs bind to the differentially expressed circRNAs and lncRNAs, respectively. Subsequently, we performed predictive analysis of the interactions between differentially expressed circRNAs/lncRNAs and target miRNAs based on the RNAseq data. The miRNA-circRNA regulatory network of granulosa cells treated with melatonin revealed 8 differentially expressed circRNAs interacting with 109 miRNAs, with circRNA-2289 and circRNA-1862 matched to 63 miRNAs (Supplementary Tables S7 and S8). In addition, the top 300 circRNAs/lncRNAs-miRNA interactions were identified based on TargetScan and miRanda analysis (Figure 5 and Figure 6).

3.6. Validation of circRNA and lncRNA Expression by qPCR

Four circRNAs (circRNA_1862, circRNA_3095, circRNA_2867 and circRNA_2319) and lncRNAs (XR_003031433.1, XR_003030822.1, XR_003035979.1, XR_001495452.2) were selected to validate the sequencing data via qPCR after melatonin treatment (Figure 7, Supplementary Table S9). The qPCR results of selected circRNAs and lncRNAs were consistent with the sequencing results of circRNAs and lncRNAs, verifying the reliability of the sequencing data.

4. Discussion

Follicular development is a complex process from small antral follicles to growing antral follicles, finally reaching the dominant follicles ready for ovulation. During this process, oocytes are surrounded by granulosa cells, and their interaction is crucial for follicular and oocyte developmental capacity [22]. However, considerable follicles undergo atresia by up to 99% during follicular development [23]. How to improve follicular development, oocyte maturation and oocyte quality is still a worthy research topic, which is crucial for female reproductive performance. It is well known that granulosa cells play an important role in the process of follicular development and oocyte maturation. FSH and LH target granulosa cells and then regulate cell proliferation and differentiation, which can induce granulosa cells, producing a large amount of fluid to promote follicular development [22]. Accompanied by the proliferation and differentiation of granulosa cells, a number of hormones and factors can be synthesized from granulosa cells, such as BMP6, estradiol, progesterone, inhibin, etc., which are closely related to follicular development and oocyte maturation [24,25]. On the other hand, follicular atresia can also be induced by granulosa cell apoptosis [26]. In addition to the factors mentioned above, it is widely accepted that melatonin is another important regulator involved in modulating the proliferation, apoptosis and hormone secretion of granulosa cells [19,27]. However, the mechanisms by which melatonin influences circRNA and lncRNA expression in bovine granulosa cells remain unclear. RNA sequencing can provide insights into the candidate circRNAs and lncRNAs related to the molecular mechanism.

Accumulating studies have demonstrated that circRNAs play a potential role in ovarian function. Meng et al. [11] showed that 62 circRNAs were differentially expressed in the porcine granulosa cells from healthy antral and atretic antral follicles. Among these, circ_KIF16B may regulate its targeted gene of the apoptotic gene TP53 and its downstream target PHLDA3, contributing to the regulation of porcine granulosa cell apoptosis. Additionally, the differentially expressed circRNAs, circRNA_103827 and circRNA_104816, were possibly involved in ovarian aging and might be potential indicators of a compromised follicular microenvironment, which could be used to improve female infertility management [10]. Recent studies have explained that circRNAs may play an important role in oocyte maturation, embryo development, follicular development and granulosa cell function. Cao et al. [28] demonstrated that circRNAs were abundantly and dynamically expressed in a developmental-stage-specific manner in cumulus cells and oocytes, and maternally expressed circARMC4 was essential for porcine oocyte meiotic maturation and early embryo development. CircDDX10 is related to the number of oocytes and good-quality embryo rate and may participate in the regulation of ovarian function and senescence by affecting the proliferation, apoptosis and steroid hormone synthesis of granulosa cells, which is expected to become a novel biomarker for predicting the outcomes of assisted reproduction techniques [29]. Moreover, circKif2, circVcan, circMast4, circMIIt10 and their target miRNAs and protein networks provide new insights into the complex interactions associated with periovulatory follicular development [30]. The expression profiles and characterization of circRNAs at two phases of follicular development, thereby, provide an improved understanding of the roles of circRNAs in the sheep hypothalamus and their involvement in follicular development and ovulation [31]. It is generally accepted that circRNAs could modulate gene expression by regulating microRNAs. Guo et al. [11] provided evidence that circINHA promoted granulosa cell proliferation and inhibited pig granulosa cell apoptosis via CTGF as a competing endogenous RNA that directly bound to miR-10a-5p to modulate ovarian physiological functions. Similarly, ssc-circINHA-001 serves as a ceRNA to guarantee the expression of INHBA by sponging miR-214-5p, miR-7144-3p, and miR-9830-5p simultaneously, thereby increasing activin production and inhibiting granulosa cell apoptosis and follicular atresia [32]. In the present study, we investigated the characteristics of circRNA expression in granulosa cells treated with melatonin via RNA sequencing. Thus, 28 differentially expressed circRNAs were found in the bovine granulosa cells treated with or without melatonin. The differentially expressed circRNAs were further analyzed to explore their potential functions via bioinformatics analysis. The GO classification terms indicated that circRNAs participated in biological processes, including biological regulation, signaling, and molecular functions, such as cell part and cell junction. KEGG analysis showed that some signal pathways, including the mTOR signaling pathway, tight junction, MAPK signaling pathway, Hippo signaling pathway and apoptosis, were enriched in granulosa cells after treatment with melatonin, most of which are crucial for regulating the bovine granulosa cell function, communication between granulosa cells and oocytes, oocyte maturation and follicular development. Therefore, the present research may provide some new perspectives to explain the mechanisms of melatonin in regulating ovarian function.

In addition, lncRNAs have exhibited important roles in development, reproduction, sex hormone responses and oocyte meiosis by regulating the expression of neighboring coding genes [33]. lncRNAs are differentially expressed in human mature MII oocytes, and these lncRNAs could be involved in many processes that regulate follicular genesis [34]. The presence and dynamics of genes that encode lncRNAs and the comprehensive transcriptome analysis of lncRNAs in primordial, primary and small antral follicles suggest that these lncRNAs may mediate functions in the cyclic recruitment and differentiation of human follicles [35]. Dong et al. [36] demonstrated that differential lncRNAs might be involved in biological processes by modulating gene expression and further providing insight into the role of lncRNA involvement in regulating reproduction in the ovary of tongue sole. The transcriptome data confirm that lncRNA MSTRG.28645 is more highly expressed in small tailed than in sheep but is lowly expressed in Xinji fine wool sheep. Furthermore, MSTRG.28645 is one of the key lncRNAs related to fecundity, which is a critical regulator in the secretion of hormones in the granulosa cells, thus affecting the fecundity of the sheep [37]. lncRNA NEAT1 regulates granulosa cell proliferation, apoptosis, and steroidogenesis via the miR-204-5p/ESR1/MAPK axis [38]. In line with these studies, our results showed 99 differentially expressed lncRNAs in the bovine granulosa cells treated with or without melatonin. The pathway analysis of the differentially expressed lncRNAs from bovine granulosa cells exhibited that apoptosis and follicular development-related pathways, including the “MAPK signaling pathway”, “FoxO signaling pathway”, “HIF-1 signaling pathway”, “NF-kappa B signaling pathway”, “TGF-beta signaling pathway”, “PI3K-Akt signaling pathway” and “Apoptosis”, were enriched. Granulosa cells are crucial for follicular development, and the main reason is that apoptosis of granulosa cells could induce antral follicular atresia [2]. Moreover, the gap junction signaling pathway is enriched in granulosa cells treated with melatonin, which is critical for follicular development and oocyte maturity through communication between the oocyte and the granulosa cells [39]. Furthermore, insulin and the mTOR signaling pathway also play important roles in follicular development [40,41]. The present analysis indicated that lncRNAs may be involved in follicular development by mediating the signaling pathway of apoptosis and development in granulosa cells in response to melatonin treatment, suggesting a potential role of melatonin in follicular development.

It is well established that melatonin plays an important role in modulating ovarian function by promoting granulosa cell proliferation, hormone secretion, as well as follicular development and oocyte maturation, which is critical for improving animals’ reproductive performance. In addition, epidemiological studies have indicated a possible oncostatic property of melatonin on different types of tumors. Melatonin could be an excellent candidate for the prevention and treatment of several cancers, such as breast cancer, prostate cancer, gastric cancer and colorectal cancer, through antioxidant activity, modulation of melatonin receptors MT1 and MT2, stimulation of apoptosis, regulation of pro-survival signaling and tumor metabolism, inhibition of angiogenesis, metastasis, and induction of epigenetic alteration [42,43]. There is a significant connection between lncRNAs/CircRNAs and the formation, growth and spread of cancer cells, indicating the potential for these lncRNAs and circRNAs to contain coding sequences for small peptides and to control their ability to multiply, move and invade other tissues, respectively [44,45]. For example, the lncRNA HOXB–AS3 can produce a 53-amino-acid peptide, inhibiting the growth of cancer cells [46]. The lncRNA TPT1–AS1 has been shown to promote the growth of tumors in ovarian cancer [47]. CircRNA contributes to modulating cell metabolism through acting as a sponge for miRNAs, with the potential role to impact the ferroptotic process [48]. CircSnx12 is involved in controlling the processes of ferroptosis in ovarian cancer through inhibiting SLC7 A11 and binding with miR-194-5p [49]. Melatonin reverses Warburg-type metabolism and possibly reduces glutaminolysis, thereby attenuating various oncogenic molecules associated with ovarian cancer progression and invasion [50]. Moreover, melatonin inhibits triple-negative breast cancer progression through the lnc049808-FUNDC1 pathway. Therefore, melatonin could be used as a potential therapeutic agent for triple-negative breast cancer [51]. Melatonin also inhibits the biological functions of osteosarcoma cells by repressing the expression of lncRNA JPX through regulating the Wnt/β-catenin signaling pathway. This suggests that melatonin might be applied as a potentially useful and effective natural agent in the treatment of osteosarcoma [52]. In addition, melatonin decreases the expression of circ_0017109 and suppresses non-small-cell lung cancer cell migration, invasion, and proliferation through decreasing TOX3 expression via direct activation of miR-135b-3p [53]. Therefore, it can be expected that melatonin influences cancer progression by modulating lncRNA and circRNA expression, which could be used as a potential therapeutic approach.

5. Conclusions

This study concludes that melatonin alters bovine granulosa cell lncRNA and circRNA expression. The present findings further indicate that melatonin may regulate granulosa cell function by mediating the expression of lncRNAs and circRNAs, thereby exerting a regulatory effect on follicular development. Collectively, these results provide insights into the candidate circRNAs and lncRNAs associated with the molecular mechanisms underlying the regulatory role of melatonin in ovarian function.