Transcriptome Analyses Identify Potential Key microRNAs and Their Target Genes Contributing to Ovarian Reserve

Female endocrinological symptoms, such as premature ovarian inefficiency (POI) are caused by diminished ovarian reserve and chemotherapy. The etiology of POI remains unknown, but this can lead to infertility. This has accelerated the search for master regulator genes or other molecules that contribute as enhancers or silencers. The impact of regulatory microRNAs (miRNAs) on POI has gained attention; however, their regulatory function in this condition is not well known. RNA sequencing was performed at four stages, 2-(2 W), 6-(6 W), 15-(15 W), and 20-(20 W) weeks, on ovarian tissue samples and 5058 differentially expressed genes (DEGs) were identified. Gene expression and enrichment were analyzed based on the gene ontology and KEGG databases, and their association with other proteins was assessed using the STRING database. Gene set enrichment analysis was performed to identify the key target genes. The DEGs were most highly enriched in 6 W and 15 W groups. Figla, GDF9, Nobox, and Pou51 were significantly in-creased at 2 W compared with levels at 6 W and 20 W, whereas the expression of Foxo1, Inha, and Taf4b was significantly de-creased at 20 W. Ccnd2 and Igf1 expression was maintained at similar levels in each stage. In total, 27 genes were upregulated and 26 genes interacted with miRNAs; moreover, stage-specific upregulated and downregulated interactions were demonstrated. Increased and decreased miRNAs were identified at each stage in the ovaries. The constitutively expressed genes, Ccnd2 and Igf1, were identified as the major targets of many miRNAs (p < 0.05), and Fshr and Foxo3 interacted with miRNAs, namely mmu-miR-670-3p and mmu-miR-153-3p. miR-26a-5p interacted with Piwil2, and its target genes were downregulated in the 20 W mouse ovary. In this study, we aimed to identify key miRNAs and their target genes encompassing the reproductive span of mouse ovaries using mRNA and miRNA sequencing. These results indicated that gene sets are regulated in the reproductive stage-specific manner via interaction with miRNAs. Furthermore, consistent expression of Ccnd2 and Igf1 is considered crucial for the ovarian reserve and is regulated by many interactive miRNAs.


Introduction
Diminished ovarian reserve, such as premature ovarian inefficiency (POI) is an increasing problem and leads to infertility in women under the age of 40 years. In previous studies, genes involved in POI have been identified (Table 1), but the critical regulators of this condition have not been fully elucidated. In addition to regulatory genes, non-coding RNAs are considered candidate regulators. miRNAs are well-known conserved regulators of gene expression [1], and their contribution to germ cell developmental control has been studied. Their involvement in ovarian follicles and oocyte development has also been demonstrated [2][3][4]. The correlation between miRNAs and POI has been studied, and several miRNAs have been suggested to be associated with POI [5]. Among these, miR-23a and miR-27a are involved in granulosa cell proliferation and apoptosis by controlling XIAP [6,7] and IGFBP2 [8].
The fundamental reasons for POI are not fully elucidated due to its complex regulation; however, advancements in array technology have enhanced the identification of major genes and other regulators, such as non-coding RNAs, including microRNAs (miR-NAs) [9,10] and long non-coding RNAs [11,12]. miRNAs are small non-coding RNAs that contribute to the regulation of gene expression at transcriptional and post-transcriptional levels [1,13]. Their regulatory functions in the ovary are also known [14,15], and their regulatory role [3,16] and expression in patients with POI have been studied [17,18].
The reproductive span of oocytes changes dramatically during the life cycle of women [19]. The number of oocytes in a female is pre-determined, and these are not reproducible and are diminished until menopause via menstrual cycles [20]. However, the reproductive cycle becomes irregular with endocrinological symptoms such as POI and other exogenous factors such as chemotherapy. With advancements in array technologies, some key genes regulating follicle development, ovarian reserve, and POI have been reported (Table 1). Table 1. List of genes related to premature ovarian inefficiency (POI) from previous studies.
Other studies have demonstrated the relationship between mRNA-miRNA interaction during ovarian follicle development and POI regulation [44,45]. Correlations between miRNAs and POI have been studied, and miR-23a, miR-27a, miR-22-3p, miR-146a, miR-196a, miR-290-295, miR-423, and miR-608 have been suggested to be associated with this condition [5]. Among these miRNAs, miR-23a and miR-27a are involved in granulosa cell proliferation and apoptosis by controlling XIAP [6,7] and IGFBP2 [8]. The apoptosis of granulosa cells is critical in POI, therefore, understanding the regulatory roles of miRNAs in granulosa cells enhances our understanding of the pathogenesis of POI [46].
The Piwil (Piwi) protein family is a small-RNA-bound effector complex [47,48], and the mouse genome encoded three Piwi proteins, PIWIL1/MIWI, PIWIL2/MILI, and PI-WIL4/MIWI2 [49]. Mammalian Piwil proteins associate with Piwi-interacting RNAs (piRNAs) and piRNA expression is largely restricted to the germline [50]. The deficiency of genes required for piRNA biogenesis leads to infertility in males [51]; however, females with this genotype are not infertile [52]. Whereas the role of regulatory miRNAs in POI has been studied, their regulatory function in this condition is not well known.
In this study, we aimed to identify the miRNAs that regulate the ovarian age-specific genes, which might contribute to the regulation of POI. We also analyzed the correlation between Piwil and specific miRNAs.

Identification of Transcripts in Different Reproductive Stages of Mouse Ovaries
We identified a total of 5058 DEGs (Figures 1 and S1) including 2658 upregulated genes and 2400 downregulated genes in the reproductively aged mouse ovary based on comparisons among stages ( Table 2). Correlations of total gene count ( Figure 2A) from each stage and filtered count numbers were determined by respective FPKM and TPM analyses and displayed by heatmaps ( Figure S2A,B). The count data displayed by the PCA plot indicate the distance of expressed transcripts ( Figure 2B). The distribution of transcripts is displayed as a cluster in Figure S3. The DEGs among the reproductively aged mouse ovaries were compared and different comparisons were divided into six groups as follows: (1)

Gene Ontology (GO) Enrichment of DEGs
The overall biological processes associated with DEGs were analyzed, and the p value cutoff was set at <0.05. Gene ontology biological process, gene ontology molecular function, and gene ontology cellular component (GOCC) analyses were performed.
In the 6 W vs. 2 W group, the DEGs were mainly enriched in the steroid biosynthetic process, regulation of cholesterol metabolic process, organonitrogen compound biosynthetic process, and cholesterol biosynthetic process ( Figure S4A). In the 15 W vs. 2 W set, DEG sets comprising the regulation of cell migration, negative regulation of the cellular process, and extracellular matrix organization were enriched. Upregulated DEGs were mainly enriched in the sterol biosynthetic process, regulation of cholesterol metabolic process, and cholesterol biosynthetic process. Downregulated DEGs were enriched in the cellular protein metabolic process ( Figure S4B). Analysis revealed that lysosome, Golgi subcompartments, and Golgi membranes were enriched among cellular components ( Figure S6A).  In the 20 W vs. 2 W set, DEGs were highly enriched in cytokine-medicated signaling pathways, and downregulated DEGs were mostly enriched in cellular protein metabolic process ( Figure S4C), neuropilin binding ( Figure S5A), and the endoplasmic reticulum lumen ( Figure S6B). In the 15 W and 6 W sets, DEGs were enriched in the cellular response to cytokine stimulus ( Figure S4D) and cytosolic part ( Figure S6C). In the 20 W vs. 6 W set, DEGs were enriched in the regulation of transcription from RNA polymerase II promoter, regulation of cell proliferation, and positive regulation of transcription ( Figure S4E). Furthermore, they were enriched in cytokine activity ( Figure S5B) and platelet alpha granules ( Figure S6D). For 20 W vs. 15 W, DEGs were enriched in the regulation of transcription from RNA polymerase II promoter and positive regulation of transcription ( Figure S4F). They were also enriched in RNA polymerase II regulatory region sequence-specific DNA binding ( Figure S5C). These data indicated that the transcriptomes in the ovary are altered and stage-specific DEGs are increased in the reproductively aged 20 W ovary.

KEGG Analysis of DEGs
The upregulated DEGs in each group were enriched in cytokine-cytokine receptors for 6 W vs. 2 W ( Figure S7A) and cell adhesion molecules for 15 W vs. 2 W ( Figure S7B). In the 20 W vs. 2 W set, DEGs were enriched in cytokine-cytokine receptors ( Figure S7 C). Hormones, steroids, ovarian signaling, and cortisol synthesis were upregulated in the 15 W vs. 6 W group ( Figure S7D). DEGs of 20 W vs. 6 W were enriched in TNF signaling, IL-17 signaling, MAPK signaling, and cytokine-cytokine receptor ( Figure S7E). TNF signaling and IL-17 signaling were highly enriched for DEGs that were upregulated ( Figure S7F).

Protein-Protein Interaction (PPI) Network Construction and Clusters Analyses
The PPI networks were built using STRING based on DEGs logFC values and are shown in Table 3 and Figure S8. The PPI networks of upregulated and downregulated DEGs were generated, and their biological processes, molecular functions, and KEGG pathways were analyzed. In the 6 W vs. 2 W set, Klk1, which has a role in the positive regulation of steroid hormone biosynthetic process, and the well-known Figla were downregulated (Table 3A and Figure S8A,B). Steroid biosynthesis-related Hsd17b7 and cholesterol metabolism-related Star were upregulated for 15 W vs. 2 W (Table 3B and Figure S8C). The RNA polymerase regulator, Atf3 was downregulated (Table 3B and Figure S8D), and this phenomenon was also observed in the 15 W vs. 6 W set. Steroid hormone biosynthesisrelated Akr1c18 and Hsd17b7 were also enriched in 20 W vs. 2 W (Table 3C and Figure S8E). Eif3j1, which has translation initiation factor activity was upregulated for 15 W vs. 6 W (Table 3D and Figure S8F). In the 20 W vs. 6 W set, Oog1 was downregulated, and it is involved female gamete generation (Table 3E and Figure S8G). For 20 W vs. 15 W, the IL-17 signaling pathway-related Cxcl1 and Atf3, and the regulation of transcription from RNA polymerase II promoter were upregulated in response to endoplasmic reticulum stress (Table 3F and Figure S8H). Table 3. Predicted protein interaction for up-and downregulated transcripts in the reproductively aged mouse ovary. Table 2. Fold changes among the groups were selected and their predictive interaction was analyzed using the STRING database. Redundant genes were indicated when they first appeared in the   Functions are indicated in previous page

BC100530
Biological process Negative regulation of endopeptidase activity Molecular Function Cysteine-type endopeptidase inhibitor activity

S100a8
Functions are indicated in previous page

Identification of Key Target Gene and Their Validation
Among the DEGs identified in each reproductive stage, we selected 90 highly expressed genes after filtering and selected 13, Bmp15, Ccnd2, Figla, Foxo1, Foxo3, Fshr, Gdf9, Igf1, Inha, Nobox, Smad3, Taf4b, and Pou5f1, as key target genes based on their correlation with the ovaries (Table 4). Their expression levels and ranked correlations were calculated and filtered using the FPKM and TPM values (p < 0.05). The expression levels of Ccnd2 and Igf1 were similar throughout the reproductivelyaged mouse ovaries. The levels of Figla, GDF9, Nobox, and Pou5f1 were decreased as reproductive aging progressed. Foxo1, Inha, and Taf4b demonstrated a more than two-fold increase in expression. Finally, in reproductively aged individuals, Fshr increased at 20 W.

Correlation with miRNAs
Interacting miRNAs with the key target genes were identified based on the ranked correlation (p < 0.05) and are listed in Table 5. Their interactions with miRNAs were confirmed using databases (Table S3). Continuously expressed genes, such as Ccnd2 and Igf1, are target genes of various miRNAs (Table 5). Genes regulating oocyte and follicle development, such as Figla, Foxo3, and Fshr, were determined to interact with only one miRNA, mmu-miR-669d-2-3p, mmu-miR-153-3p, and mmu-miR-670-3p, respectively.
In total, 27 genes were increased based on the regulatory effects of miRNAs, and 26 genes were decreased between the 2 W and 20 W sets (p < 0.05, Table S4A). During the reproductive period, some genes were downregulated following transient up-regulation. The genes expressed at high levels at 20 W and decreased at 2 W (Table S4B, left panel), as well as the genes expressed at high levels at 2 W and decreased at 20 W, set are listed in Table S4B, right panel. Genes upregulated following downregulation in the 20 W set (Table S4C, left panel) and 2 W set (Table S4C, right panel) are also listed with their interacting miRNAs.

Interaction with Piwil Gene Family
The Piwil gene family is essential for developmental regulation. Therefore, we additionally analyzed the correlation between DEGs from reproductively aged mouse ovaries and this family (Table S6A). Target genes and each Piwil gene were analyzed and are listed in Table S6B. Target genes of Piwil included well-known early oocyte development and follicle assembly genes (Table S6B). Piwil4 was determined to interact with mmu-miR-210-5p and mmu-miR-3470b, whereas Piwil1 and Piwil2 were found to interact with more miRNAs (Table 6). Table 6. Correlation of Piwil gene family expression and frequently expressed target genes and microRNAs in the reproductively aged mouse ovary. Target genes of the Piwil gene family and interactive miRNAs of each Piwil gene were analyzed based on the FPKM correlation (p < 0.05, score 90+).

Discussion
In this study, we demonstrated the stage-specific upregulation and downregulation of DEGs and their correlation with miRNAs in reproductively aged mouse ovaries. We identified 24,958 transcripts from each reproductive stage of ovaries, and 5058 genes were identified. In total, 2658 genes were upregulated and 2400 genes were downregulated. Furthermore, 53 upregulated and downregulated miRNAs were identified. GO analyses demonstrated that upregulated DEGs were enriched in different cellular processes in a stage-specific manner. DEGs from the early phase set of 2 W, ovaries were mainly enriched in the steroid biosynthetic process, regulation of cholesterol metabolic process, organonitrogen compound biosynthetic process, and cholesterol biosynthetic process. In ovaries from middle-phase sets, 6 W and 15 W, DEG sets associated with the regulation of cell migration, negative regulation of the cellular process, and extracellular matrix organization were enriched. KEGG analyses suggested that the upregulated genes were mostly enriched in cell adhesion molecules.
We also focused on the regulatory roles of miRNAs in reproductively aged mouse ovaries. The roles of miRNAs in folliculogenesis [4] and oocyte development have emerged [3]. We focused on the altered gene expression in a specific manner, and interactive regulatory miRNAs of target genes were identified based on a score >90 (p < 0.05). mmu-miR-136-5p, mmu-miR-335-5p, mmu-miR-665-3p, and mmu-miR-18a-5p targeted several upregulated genes in 20 W mouse ovaries (Table S5A). miRNAs targeted genes were involved in spermatogenesis, follicular development process. ESR1, SMAD2, and Spata1 are the target of mmu-miR-18a-5p. ESR1 has been involved in the genetic variation of female infertility, and Sparta is an acrosomal protein in sperm and related to genetic mutation in sperm. Its role on female infertility is still unknown. mmu-miR-136-5p targets Prdm16 and GDF6, genes involved in follicular development. mmu-miR-335-5p targets Calu, regulator of chromosome condensation. The major target of mmu-miR-665-3p is Pou2f3, which belongs to same family of Oct4, master regulator of pluripotency.
Among the highly expressed DEGs, well-known ovary-associated genes were identified. The master regulator of oocyte development, Figla was downregulated in later-stage ovaries. Figla (factor in the germline alpha), encodes a germ cell-specific basic helix-loophelix transcription factor first identified as an activator of oocyte genes [53,54]. Figla regulates primordial follicle formation in the fetal ovary [55] and antagonizes spermatogenesis during embryo development [56]. This gene is related to POI, and its expression was decreased according to reproductive age. Its interacting miRNA was mmu-miR-669d-2-3p, and the main function of this is to modify sperm-related gene expression [57]. This result supports the known function of Figla and supports its decreased expression throughout the reproductive span.
Foxo3 and Fshr are also mainly regulated by miRNAs. Foxo3 is involved in the apoptosis of granulosa cells [58] and involved in the early development of follicles on its own [59] or with other genes [60]. Foxo3 interacts with mmu-miR-153-3p, which is known to suppress tumor growth in ovarian carcinoma [61]. Fshr (follicle-stimulating hormone receptor) is one of the most important receptors in the ovary and other female reproductive organs. It is associated with the granulosa cells of the follicle and is located in the granulosa cells [62]. Interactive mmu-miR-670-3p is expressed in the newborn ovary [63,64], however, its roles are still not well known.
Klk1 has been involved in positive regulation of the steroid hormone biosynthetic process, is related to the renin-angiotensin system, and was enriched in the early and middle phase sets. The renin-angiotensin system, with angiotensin, is known for its relation to ovarian follicle development [65] and acquisition of dominancy [66]. Components are expressed in an estrogen-dependent manner in the uterus and are involved in cell proliferation [67]. Furthermore, they also involved the relationship with umbilical veins [68].
The Piwil family has been identified as comprising regulatory proteins responsible for stem cell and germ cell differentiation [69], and it is co-expressed with early developmental genes [70]. To date, the Piwil family is known to be more involved in male fertility than female fertility [50]. POI does not result in infertility from the beginning; therefore, we hypothesized that fine regulatory components might be involved in the onset of POI. Furthermore, they would be responsive to E2, which is another important molecule in the ovary [71]. The main targets of this family are enriched in genes such as FOXL2, NANOG, POU5F1, LIN28A, SOX2, PRDM1, NANOS3, DAZL, and DDX4 (Table 6). These results indicate that Piwil might contribute to the regulation of certain DEGs in a stagespecific manner.
Among this family, Piwil4 interacted predominantly with two miRNAs, mmu-miR-210-5p, and mmu-miR-3470b. Further, mmu-miR-26a-5p, miRNA that was determined to interact with Piwil2, mostly interacted with downregulated genes, such as APCDD1, CTDSP2, EZH2, MGA4A, PATZ1, and RHOQ, based on the 2 W vs. 20 W comparison (Table S5 B). These results indicate that Piwil2 and mmu-miR-26a-5p might be a candidate that regulates POI. Piwi-interacting RNAs (piRNAs) are involved in infertility in males, but not in females. We hypothesized that Piwil2 does not lead to terminal infertility; however, it affects the ovarian reserve via downregulated genes. APCDD1 is involved in adipogenic differentiation [72], CTDSP2 is a target gene of FOXO and regulates cell cycle progression [72], EZH2 is an important regulator in the female reproductive tract [73], MGA4A is involved in embryo lethality and female infertility [74], and RHOQ regulates mitochondrial function [75]. All genes regulated by miR-26a-5p were related to female infertility and cell cycle progression. Therefore, their correlation with POI needs to be further studied. These findings also correlated with the target genes regulated by mmu-miR-26a-5p were male fertility-related genes.
In conclusion, we identified DEGs and interactive miRNAs in reproductively aged mouse ovaries. Stage-specific upregulation and downregulation of DEGs indicated the regulatory roles of miRNAs in the ovary and their correlation with ovarian reserve. We used three reproductively aged stages, which correspond to the reproductive span of a human. However, differences between species limit the interpretation of experimental results. POI is an endocrinological symptom only in women; therefore, identified miRNAs might have distinct roles in humans. Further studies are needed to elucidate the functions of these miRNAs as regulators of the ovarian reserve. We will further study the roles of these miRNAs during in vitro follicular development based on non-human primate models, which are physiologically closer to humans [76][77][78].

Ethics
All the experiments were conducted under the control of AAALAC guidelines, and the animal experimental plan was reviewed and approved by IACUC of Seoul National University Hospital (No.18-0029-C1A1).

Preparation of Ovaries and Extraction of Total RNA
The experimental scheme of this study is presented in Figure 1. Ovaries of 2, 8, 15, and 20 weeks old C57BL/6 strain female mice were collected immediately after cervical dislocation. The lipid pad near the ovaries was removed and the collected ovaries were washed with HBSS (Invitrogen, Waltham, MA, USA). And then, ovaries were minced into pieces by a surgical blade (Feather Safety Razor, Osaka, Japan). The pieces were digested with Collagenase Type IV (1 mg/mL, Invitrogen) at 37 • C for 45 min and centrifugated. Pelleted samples washed with PBS (Invitrogen) and total RNAs were extracted using Trizol kit (Invitrogen) according to manufacturer's instruction. Extracted RNAs were quantified and qualified for further Total Omics Transcriptome analyses.

mRNA Sequencing
The libraries were prepared for 150 bp paired-end sequencing using TruSeq Stranded mRNA Sample Preparation Kit (Illumina, CA, USA). Namely, mRNA molecules were purified and fragmented from 1µg of total RNA using oligo (dT) magnetic beads. The fragmented mRNAs were synthesized as single-stranded cDNAs through random hexamer priming. By applying this as a template for second-strand synthesis, double-stranded cDNA was prepared. After the sequential process of end repair, A-tailing, and adapter ligation, cDNA libraries were amplified with PCR (Polymerase Chain Reaction). The quality of these cDNA libraries was evaluated with the Agilent 2100 BioAnalyzer (Agilent, CA, USA). They were quantified with the KAPA library quantification kit (Kapa Biosystems, MA, USA) according to the manufacturer's library quantification protocol. Following cluster amplification of denatured templates, sequencing was progressed as paired-end (2 × 150 bp) using Illumina NovaSeq 6000 sequencer (Illumina, CA, USA).

microRNA Sequencing
The libraries were prepared for 50 bp single-end sequencing using the NEXTflex Small RNA-Seq Kit (Bioo Scientific Corp). Namely, small RNA molecules were isolated from 1 µg of total RNA via the adapter ligation. The isolated small RNAs were synthesized as single-stranded cDNAs through the RT (Reverse transcription) priming. By applying this as a template for second-strand synthesis, double-stranded cDNA was prepared through PCR (Polymerase Chain Reaction). And, the fragments around 150 bp were extracted for sequencing through size selection by gel electrophoresis. The quality of these cDNA libraries was evaluated with the Agilent 2100 BioAnalyzer (Agilent, CA, USA) followed by quantification with the KAPA library quantification kit (Kapa Biosystems, MA, USA) according to the manufacturer's protocol. Following cluster amplification of denatured templates, sequencing was progressed as single-end (50 bp) using Illumina sequencing platform (Illumina, CA, USA).

Read Quality of the Experiment
The identification was performed by the Sanger method, using Illumina 1.9. Total sequences were 32219029 before fast qc and 31058105 after fast qc. Sequence length ranged from 36 to 101 ( Figure S1A). The percentage of the GC was 49 and the GC content has not emerged from theoretical distribution ( Figure S1B). No rRNA contamination was observed. Sequences filtered as the poor quality was none ( Figure S1C) and the duplicated sequence was 36.32% ( Figure S1D).

Data Processing
We evaluated the data quality by sample clustering based on Pearson's correlation matrices between different samples. And a heatmap was drawn corresponding to the different expressions of probes.

Differentially Expressed Genes (DEGs)
We employed the "limma" R language package to screen the DEGs between uterine leiomyoma and normal myometrium. The adjusted p-value < 0.05 and |log2fold change (FC)| > 1 were considered statistically significant.

KEGG and GO Enrichment Analyses of DEGs
Gene Ontology (GO, http://geneontology.org/, accessed on 1 February 2021) provides an ontology of defined terms to represent gene functions (molecular function, cellular component and biological process). Besides, Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/ accessed on 1 March 2021) is a database resource for understanding high-level functions and utilities of the biological system. A package in R language called "clusterprofiler" was used to determine the biological significance of DEGs. The package is capable of providing GO and KEGG enrichment analyses and visualization for users to obtain more valuable biological information [36]. p value < 0.05 was considered a significant enrichment.

Gene Set Enrichment Analysis (GSEA)
To explore the potential functions of selected key genes and microRNAs in mouse ovaries, and samples of datasets divided into different groups following the expression levels of the key genes, respectively. GSEA was utilized to explore whether priority determined biological processes datasets were enriched in these groups derived from DEGs between the two groups [30]. The criteria were set as p-value < 0.05 and FDR < 0.25.

Protein Interaction
Prediction of interaction with other functional genes was predicted by STRING database (https://string-db.org, accessed on 1 April 2021).