Dynamic Changes of DNA Methylation and Transcriptome Expression in Porcine Ovaries during Aging.
The major challenges in women's reproductive health are the reduction of reproductive performance along with aging; moreover, oocyte quantity and oocyte quality are closely related to the reduction of follicular reserve in the ovary . The decrease of follicular reserve in the ovary is nonlinear and accelerated with age [2-4]. This leads to near-complete exhaustion by a mean age between 51 and 52 years, which is defined as menopause . Ovarian aging is affected by numerous factors and has been particularly linked to genetics [5, 6]. In conclusion, the heredity record of normal reproduction and numerous pathologies like premature ovarian insufficiency (POI) and polycystic ovary syndrome (PCOS) have been emphasized [6, 7].
The concept of epigenetics related to heritable changes in chromatin structure and gene expression does not involve changes of DNA sequence. The known classes of epigenetic modification are DNA methylation, histone modification, and the synthesis of noncoding RNA, including that of microRNAs (miRNAs) and long noncoding RNAs (lncRNAs) [8, 9]. Furthermore, these three epigenetic mechanisms are in fact forming a network . "Ovarian epigenetics" is a new field that has uncovered stimulating revelations. Recent research has documented that DNA methylation plays a key role in the regulation of ovarian cancer and ovarian diseases such as PCOS and POI [10, 11]. Moreover, several noncoding RNAs such as miRNAs and lncRNAs are crucial for the regulation of ovarian physiology [12-16]. However, these results mainly elucidate the epigenetic mechanism in ovarian health-associated phenotypes. For natural ovarian aging, epigenetic mechanisms in the ovarian context have been studied far less, due to the timing of menarche and menopause. Importantly, epigenetic modifications during the natural aging of ovaries are a perfect opportunity to understand the health-related phenotypes of ovaries. The reason is that these studies not only pointed out the mechanisms of ovarian aging but also elucidated the complex interaction networks of different ovarian phenotypes.
It is both difficult and unethical to investigate the ovarian mechanisms in women. While laboratory rodents (e.g., mice) are useful models for biomedical research, they offer comparatively limited use for the study of ovarian changes in mammals, due to their small body size and extremely short ovarian cycle . The big animal models such as equine, bovine, and ovine have been confirmed to be valuable and useful for investigation on the ovarian function in women . Therefore, pigs can be a valuable model to study human ovarian aging or disease due to its similar cycle length, luteinizing hormone (LH) receptor location and its function, length of ovulation, and LH surge  as well as their anatomical, physiological, and biochemical similarities to humans . In fact, a number of researchers have identified pig as an ideal model system to investigate the effects of metabolic syndrome and obesity in relation to the function and steroidogenesis of ovaries [21, 22].
Female natural ovarian aging is defined between the two time points (menarche and menopause) in a woman's life that open and close the reproductive system. The median ages for menarche and menopause are about 14 and 50 years, respectively . Both time points are central to ovarian function, and several recent genome-wide association studies (GWAS) explained the genetic background of traits, both for the timing of menarche [23-25] and menopause [26-28]. Correspondently, the puberty stage of the first ovulation and the reproductive exhaustion stage are two key points in the ovarian cycle. The median ages of the sow puberty stage of first ovulation and reproductive exhaustion stage are about 180 days and eight years, respectively [29, 30].
To identify the potential molecular changes during natural ovarian aging, the pig was used as a model animal. Genome-wide DNA methylation and transcriptome-wide RNA expression analyses were performed via highthroughput sequencing of ovaries from young (180 days, puberty stage of first ovulation) and old (eight years, reproductive exhaustion stage) sows. This enabled the determination of a number of differentially methylated and expressed genes or regulatory elements and previously unclear multigroup cooperative control networks in the ovarian aging cycle. These data will help to explain ovary age-associated changes in DNA methylation and transcriptional patterns over time.
2. Materials and Methods
2.1. Animal Material and Sample Preparation. Four healthy female Yanan pigs were utilized in this experiment. The Yanan pig is an indigenous Chinese pig breed, which emerged in the hilly areas of western Sichuan Province in the past. Due to the poor growth performance and carcass composition, it is endangered by extinction. These four pigs included two eight-year-old sows at the reproductive exhaustion stage and two 180-day young sows at the puberty stage of first ovulation. These pigs had no direct or collateral blood relationship during the last three generations. Piglets were weaned at the age of 28 [+ or -] 1 days. An initial diet started from the 30th to the 60th day after weaning and contained 3.40Mcahkg"1 of metabolizable energy with 20.0% crude protein (11.5 g/kg lysine). From the 61st to the 120th day, pigs were given a diet containing 14.0 MJ/kg of metabolizable energy comprising 18% of crude protein (9.0 g/kg lysine). From the 121st day, they received a dietary metabolizable energy and crude protein (8.0 g/kg lysine) of 13.5 MJ/kg and 16%, respectively. Pigs were allowed to access water and food ad libitum and were kept under similar conditions. The night before slaughtering, pigs were not allowed to feed and were given 2 h rest after transportation, then stunned electrically at 90 V and 50 Hz for 10 s, and exsanguinated to ameliorate pain. All animal experiments and procedures were conducted according to the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in June 2004) and were approved by the Institutional Animal Care and Use Committee in College of Animal Science and Technology, Sichuan Agricultural University, Sichuan, China, under permit no. SKY-S20150804. All research animals were obtained from Sichuan Weimu Modern Agricultural Science and Technology Co., Ltd., Chengdu, Sichuan 611536, P. R. China. Ovary tissues were rapidly collected from every carcass and directly frozen in liquid nitrogen after separation. All collected samples were stored at -80[degrees]C until the extraction of DNA and total RNA.
2.2. Whole Genome Bisulfite Sequencing (WGBS) and Data Analysis. Bisulfite treatment of 50-100 ng of purified genomic DNA was performed using the Zymo EZ DNA Methylation Lightning Kit. 50 to 100 ng of purified genomic DNA was treated with Zymo Lightning Conversion Reagent for 8 min at 98[degrees]C in a thermal cycler and then for 60 min at 54[degrees]C. Bisulfite-treated DNA was purified via spin column and was used to build a sequencing library with the help of the EpiGnome[TM] Kit (Epicentre). In this process, bisulfite-treated single-stranded DNA was arbitrary primed using a polymerase with the ability to read uracil nucleotides, to manufacture DNA with a particular sequence tag. 3' ends of the newly manufactured DNA strands were tagged with a 2nd specific sequence, resulting in ditagged DNA molecules at both their 5' and 3' ends with identified sequence tags. Then, these tags were used to add Illumina adapters P7 at the 5' and P5 at the 3' end of the original DNA strand by polymerase chain reaction (PCR). Only the complement to the original bisulfite-treated DNA was used as sequencing template; therefore, the resulting read always had a similar sequence like the original bisulfite-treated strands. Then, bisulfite genome DNA libraries were sequenced by the Illumina Hiseq 4000 platform with 150PE reads.
Reads were alimented to the Suscrofa 10.2 with Bismark tools after quality filtering . Bowtie2 was called by Bismark when mapped reads to reference. Parameters were multiseed length of 20 bp with 0 mismatches and minimum alignment score function L, 0, and -0.2 . For 150 bp reads, this would mean a minimum alignment score of -30 before an alignment would become invalid. This is approximately equal to four mismatches or -4 InDels of 1-4bp in the read . After reads were alimented, SAMtools were applied to deal with the bam file out form Bowtie2 . CpG island locus information was downloaded from the UCSC genome browser, and CpG island methylation level of each sample was calculated with bedtools . To identify differentially methylated regions between both stages, the R package DSS was used to call differentially methylated loci with p.threshold <0.001 first and then different methylation regions (DMR) with delta <0.1 with p.threshold < 0.001 [34, 35]. The function gene ontology (GO) of the DMR target gene was enriched by topGO . At last, other data statistics and visualizations were conducted in R and Python scripts.
2.3. RNA Sequencing and Data Analysis. A total amount of 5 y.g RNA of each sample was used as input material for RNA sample preparations. The rRNA-depleted RNA was used to generate sequencing libraries by NEBNext[R] Ultra[TM] Directional RNA Library Prep Kit for Illumina[R] (NEB, USA) according to the instructions of the manufacturer, and the Agilent Bioanalyzer 2100 system was used to assess library quality. After cluster generation, Illumina HiSeq 4000 was used for the sequencing of libraries and paired-end reads were generated at 150 bp. After removing the adapter, ploy-N, and low-quality reads from raw data, clean reads were obtained. These clean reads were aligned to the Ensemble (Susscrofa 10.2) using TopHat2 (v2.0.14) with default parameters .
StringTie software was used to assemble the mapped reads per sample , which worked in at least one of both replicates. The obtained transcripts were blasted (e value = 1e - 10) to Ensemble, and mapped transcripts were directly described as known lncRNA or mRNA. Salmon (v0.6.0) was used to calculate transcripts per million (TPMs) of both lncRNAs and mRNAs per sample . Then, coding potential calculator (CPC, 0.9) and Pfam Scan v1.5 were used to examine the transcript's coding potential [40, 41]. Transcripts predicted with coding potential were filtered out, and the transcripts without coding potential were regarded as candidate set of novel lncRNAs.
2.4. Small RNA Sequencing and Data Analysis. A total amount of 5 y.g RNA of each sample was used as input material for library preparation of small RNA. Sequencing libraries were created using NEBNext[R] Multiplex Small RNA Library Prep Set for Illumina[R] (NEB) according to the instructions of the manufacturer. In every sample index, codes were added to attribute sequences. The Agilent Bioanalyzer 2100 system was used to analyze library quality. The clustering of index-coded samples was performed on a cBot Cluster Generation system using TruSeq SR Cluster Kit v3-cBot-HS (Illumina), following the manufacturer's recommendations. After the generation of clusters, the preparations for library sequencing were conducted on an Illumina MiSeq platform and single-end reads were generated at 50 bp.
miRBase 21 was used as reference, and the software mirdeep2 was used to obtain the potential miRNA and to predict novel miRNA . The miRanda (v3.3a, default parameters) and cutoffs (score S [greater than or equal to] 140 and energy E [less than or equal to] -20.0) were used to predict the miRNA target . Then, miRNA expression was assessed by TPMs through the following criteria: normalization formula: normalized expres sion = mapped read count/total reads * 10,00,000.
2.5. Differential Expression Analysis and Gene Ontology Enrichment Analysis. Differentially expressed mRNAs, lncRNAs, miRNAs, and circRNAs were found by using the edgeR package . A differential expression P value <0.05 and fold change >2 were assigned as differentially expressed in different comparisons. GO enrichment was performed by TopGO . Other data statistics and visualizations were performed by self-written R scripts.
2.6. Network with MicroRNA Response Element (MRE) and Coexpression Network. The network of RNAs with MRE was established by target prediction of miRNA-mRNA, miRNA-lncRNA, and miRNA-circRNA with bioinformatics and visualized with the R package igraph. The coexpression network of mRNA, miRNA, lncRNA, and circRNA (circular RNAs) was created based on Pearson's correlation coefficient of expression. Pearson's correlation coefficient of 0.85 between two RNAs was considered relevant for network construction. A P value below 0.05 was considered statistically significant.
2.7. Bisulfite Sequencing PCR (BSP). The primers for BSP were designed by primer software V5.0 (Supplementary Table S1). The inspected DNA (bisulfite conversion) was treated using the EpiTect Fast DNA Bisulfite Kit (Qiagen) according to the manufacturer's protocols. PCR product was purified using the UNIQ-10 Spin Column DNA Gel Extraction Kit for PAGE (Sangon) and was then cloned with the pGM-T Fast Cloning Kit with competent cell (Tiangen). Ten effective clones were selected per gene, and then, an ABI 3730 DNA sequencer was used for sequencing. DNAMAN 7.0 (Lynnon Biosoft, USA) was used to analyze all sequences.
2.8. Quantitative Real-Time PCR (q-PCR). Total RNAs were extracted from ovaries using the HiPure Universal RNA Mini Kit (Magen, China) and reversely transcribed into cDNA using the oligo (dT) and random 6-mer primers, provided by the PrimeScript RT Master Mix kit (TaKaRa). The q-PCR was performed using a standard SYBR Premix Ex Taq kit (Takara, Dalian, China) on a Bio-RadCFX96 RealTime PCR Detection system (Bio-Rad, Hercules, CA, USA) following the manufacturer's directions. Three endogenous control genes, glyceraldehyde-3-phosphate dehydrogenase (GAPDH), [beta]-actin (ACTB), and small nuclear RNA (U6 snRNA), were used in this assay. The [2.sup.[DELTA][DELTA]Ct] method was used to determine the expression levels of objective mRNAs, miRNAs, lncRNAs, and circRNAs . These primers are shown in Supplementary Table S1.
3.1. Summary of Whole Genome Bisulfite Sequencing (WGBS) Data. To assess the changes in DNA epigenetic marks by aging that occurred during ovarian development, DNA methylation profiles were determined using WGBS between old pig (OP) ovarian tissues at reproductive exhaustion stage and young pig (YP) ovarian tissues at the puberty stage of first ovulation. Approximately 951,440,304 clean reads were obtained from the four ovarian samples, which provided about 30x sequencing depth. After removing unclear reads from clean reads, approximately 61-68% reads per sample were uniquely aligned to Ensemble Susscrofa 10.2 (Supplementary Table S2). The global methylation levels were 70.8% for YP and 72.9% for OP, respectively. There was no significant difference between the two ovarian development stages (Figure 1(a)), although a genomewide loss of DNA methylation was found in response to agerelated epigenetic modification . Analysis of the sequence context of cytosines showed methylation that occurred in the CpG and non-CpG (CHG and CHH) context in most chromosomal regions in each group. Differential methylation between both ovarian developmental stages mostly occurred at CpG context regions (Supplementary Table S3 and Supplementary Figure S1).
To understand the preferential location of the CG methylation on and surrounding the gene body region (GBR), the metagene profiles of CG methylation were investigated in the entire pig genome of the ovarian tissue. Both GBR and adjacent intergenic regions were heavily methylated; however, slightly higher levels of methylation were recorded in gene bodies than in neighboring intergenic regions (Figure 1(b)), which have similar distribution to those of human primordial germ and prenatal germline cells [47, 48] as well as pig skeletal muscle cells .
To examine the dynamics of methylation on a global scale, correlation analysis between genomic features and methylation levels was conducted. Negatively correlated methylation levels were observed across chromosomes with chromosomal length (Pearson's r = -0.636, P = 1.844 x [10.sup.-9]), and a positive correlation was observed with the GC content (the percentage of guanine and cytosine, r = 0.903, P < 2.2 x [10.sup-16]), gene number (being calculated for a 1 Mb window in chromosomes, r = 0.398, P < 2.2 x [10.sup.-16]), Cp [G.sub.o/e] (the ratio between the observed and expected numbers of CpG sites, r = 0.792, P < 2.2 x [10.sup.16]), and repeat region density (being calculated for a 1 Mb window in chromosomes, r = 0.126, P < 2.2 x [10.sup.-16]) (Supplementary Figure S2). These results are similar with previous reports [50, 51]. Among these genomic features, the GC content and CpGo/e showed the strongest correlations with the methylation level. The gene density also showed a moderate correlation with the methylation level, which may be due to the higher GC content examined in gene regions, which contributed to this (Figure 1(b)), and suggests a possible role of methylation dynamics in the gene transcription regulation .
3.2. Differential DNA Methylation Associated with Ovarian Aging. To further discover the differential CG methylation across ovarian aging, different methylation regions (DMRs) were identified between OP and YP. As a result, 422 DMRs were identified between both developmental stages and 303 of these DMRs were upmethylated while 119 DMRs were downmethylated in OP compared to YP (P < 0.001, Supplementary Database S1). The dynamical DMR level at both ovarian developmental stages suggested that DNA methylation may play a crucial role in ovarian aging. Furthermore, 146 of these 422 DMRs overlapped at gene body regions while only 12 of these DMRs were located in the gene promoter regions (Supplementary Database S1). This agreed with the preferential location of the CG methylation on the gene body regions (Figure 1(b)). Previous research also reported that more DMRs are enriched in gene bodies than in promoters of porcine skeletal muscle .
To characterize the role of genes associated with these DMRs, GO enrichment was performed. The results showed that hypermethylated genes related to DMRs were involved in numerous cellular functions like protein binding, death effector domain binding, and cysteinetype endopeptidase activity involved in the apoptotic signaling pathway (Figure 1(c) and Supplementary Database S1). Notably, hypomethylated genes that were related to DMRs depicted significant enrichment for different processes associated to embryonic skeletal/brain system development, embryonic digit morphogenesis, negative regulation of immune response, and apoptotic signaling pathway. For instance, the CASP8- and FADD-like apoptosis regulator gene (CFLAR) is an important molecule of the innate immune regulation network and a key suppressor of steatohepatitis  and was significantly upmethylated in old pigs. However, the Meckel syndrome Type 1 gene (MKS1) in Meckel-Gruber syndrome causes developmental malformations and cilia defects  and was significantly downmethylated in old pigs. Furthermore, the BSP results for methylation levels of the two genes were in accordance with the WGBS data between OP and YP (Supplementary Figure S3A).
3.3. DNA Methylation and Gene Expression in Gene Body. The effect of methylation on promoter regions was reported to be an important mechanism in regulating the gene transcription . However, the precise roles of DNA methylation in gene body are yet to be discovered. To search how gene expression affected intergenic methylation, the RNA sequence data were used (Supplementary Database S2) to correlate DMR-mRNA pairs. A significant negative correlation (r =-0.179, P = 5.72 x [10.sup-7]) was found between changes in the methylation levels in the gene bodies and gene expression levels. Similar to these findings, previous methylation analysis in gene bodies described a significantly negative correlation with the expression levels of mRNA . However, previous scientists reported a positive correlation with gene expression levels [56, 57] or no clear relationship patterns .
3.4. Transcriptome Profiles of Ovarian Aging. To assess transcriptional expression changes during ovarian aging, RNA and small RNA libraries were constructed for OP and YP samples, respectively, and transcriptome-wide profiling (mRNA, miRNA, lncRNA, and circRNA) was determined via high-throughput sequencing. For RNA sequencing libraries, an average of -75 million clean reads were obtained from each of the four samples and more than 69% of these reads could be uniquely aligned to the Ensemble Susscrofa 10.2 (Supplementary Table S2). Furthermore, for small RNA sequencing libraries, approximately 10.89-13.50 million clean reads were obtained from each of the four samples and 80.77-90.5% of these reads were uniquely aligned to the Ensemble Susscrofa 10.2 (Supplementary Table S2).
In total, 20,357 mRNAs were identified in these four samples, representing approximately 59.82% of the entire number of transcripts in pigs (Supplementary Database S2). Moreover, 1,196 miRNAs were identified in these four samples, and 869 potential novel miRNAs out of these miRNAs were detected that did not match the previously reported sequences (Supplementary Database S2). A total of 4,879 lncRNAs and 7,600 circRNAs were also identified in these four samples (Supplementary Database S2).
3.5. Differentially Expressed Transcriptomes Involved in Ovarian Aging. According to the chosen screening criteria (P [less than or equal to] 0.05), 2,243 differentially expressed genes (DEGs) were found in both stages of ovarian development, which included 1,660 upregulated genes and 583 downregulated genes (Figure 2(a) and Supplementary Database S3). GO analysis showed that these DEGs were mainly enriched in the extracellular region and involved in a variety of cellular functions including the defense response to virus, regulation of cell death, embryonic pattern specification, and reproduction process (Figure 3(a) and Supplementary Database S4). For instance, ISG15 (ISG15 ubiquitin-like modifier) repressed interferon-[alpha]/[beta] overamplification and autoinflammation , and the higher expression level of ISG15 in old pigs suggested that old females were more susceptible to viral infection. Furthermore, SRY-Box 9 (SOX9) battles of sexes with Forkhead Box L2 (FOXL2), which exert a decisive role during ovary maintenance process and somatic sex reprogramming of adult ovaries to tests . The higher expression level of SOX9 in OP indicated that the reproductive performance decline in old females.
Ninety-five differentially expressed miRNAs (DEMs) were identified during both stages of ovarian development, and 37 miRNAs of these DEMs were upregulated while 58 miRNAs were downregulated (Figure 2(b) and Supplementary Database S3). GO analysis showed that the target genes of DEMs are enriched in the G-protein coupled receptor signaling pathway, apoptotic signaling pathway, female pregnancy, fertilization process, embryonic development, and ovulation cycle (Figure 3(b) and Supplementary Database S4). For example, the upregulated miR-9 plays a key role in the determination of the neural fates in embryonic stem (ES) cell differentiation  and has prospective importance in recurrent ovarian cancer as biomarkers .
For differentially expressed lncRNAs (DELs) related to ovarian aging, 248 DELs were identified in both ages, which included 202 upregulated DELs and 46 downregulated DELs (Figure 2(c) and Supplementary Database S3). GO analysis showed that these target genes of DELs were significantly enriched for modulation by symbiont of host I-kappaB kinase/NF-kappaB cascade, meiotic cell cycle process involved in oocyte maturation, female pregnancy, and uterus development (Figure 3(c) and Supplementary Database S4). Furthermore, 116 differentially expressed circRNAs (DECs) were identified in both ages; 103 circRNAs of these DECs were upregulated while 13 circRNAs were downregulated (Figure 2(d) and Supplementary Database S3). These DEC source genes were enriched in transmembrane receptor protein serine/threonine kinase activity, in utero embryonic development, reproductive process, ovarian cumulus expansion, and ovulation cycle (Figure 3(d) and Supplementary Database S4).
3.6. Multiomics Coordinated Regulation in Ovarian Aging. To discover the coordinated regulation relationship among DNA methylation and several RNA species during ovarian aging, the number of genes and their percentages were first measured in every potential combination of regulation. Among DNA methylation, mRNA, and miRNA, 9,633 (OP) and 8,845 (YP) genes (representing 38.2% and 34.9% of the entire number of genes in swine genome, respectively) were methylated or simultaneously expressed. Except for each of these three combination groups, the number of genes that were methylated or expressed at the same time decreased and approximately 600 genes were not methylated or expressed simultaneously in OP or YP (Supplementary Table S4). Similarly, among DNA methylation, lncRNA, and miRNA, 1,465 (OP) and 1,354 (YP) genes (representing 18.1% and 16.7% of the entire gene number in swine genome, respectively) were methylated or simultaneously expressed. The lowest number of genes (175 for OP and 232 for YP) was not methylated or expressed simultaneously (Supplementary Table S4). The combination regulation relationships of gene transcriptional regulation and gene posttranscriptional regulation were also shown in each chromosome (Supplementary Figure S4). The results indicate that a large number of genes tended toward combination regulation between DNA methylation and transcriptome expression. The utilized combination regulation pattern might be more pronounced in OP. In addition, the DMRs and the differentially expressed transcriptome were screened to identify the overlap between DMRs and DEGs. The results showed that seven DEGs were differentially methylated during ovarian aging, most of which were found in the genebody regions, and only one gene was methylated in the gene promoter region (Supplementary Table S5). However, no consistent pattern of change was found between methylation and expression levels in these methylated genes in genebody regions. Interestingly, the gene that was methylated in the promoter region was hypomethylated and upexpressed simultaneously during aging in the ovary. These results indicate that changes of methylation occurred in the gene promoter region rather than the genebody region, which might be close to the regulation of the gene expression level during aging.
To further highlight the potentially coordinated control roles of DNA methylation and several RNA species involved in ovarian aging, GO terms of genes related to DMRs and four differentially expressed RNA species were overlapped. Among DMRs and these differentially expressed RNAs, several intersections of molecular functions and signaling pathways were found that are involved in the ovarian aging cycle such as embryonic development, apoptotic process, reproduction and fertilization, ovarian cumulus expansion, and female pregnancy (Figure 4 and Supplementary Database S4).
Based on the competing endogenous RNA (ceRNA) hypothesis , ceRNA networks were constructed among these differentially expressed RNAs, which shared a common binding site of the MRE (Figure 5(a) and Supplementary Database S5). The network consists of a large number of interrelated RNAs that include 20 miRNAs, 170 mRNAs, 27 lncRNAs, and eight circRNAs, and miRNAs played a primary role by targeting other RNA species with MRE in the network. At the upper part of the network, all novel miRNAs were only targeted by circRNAs; however, at the lower part of the network, three miRNAs (miR-9, miR-91, and miR-9-2) could simultaneously target nine mRNAs and one circRNA, respectively. In the center of the network, several miRNAs were found including miR-125b, miR-504, and miR-92b-5P, suggesting important roles during the ovarian aging cycle. In fact, previous studies showed that these miRNAs exert important effects on ovarian function. For instance, miR-125b was a negative regulator for p53induced apoptosis , promoted both the proliferation and migration of type II endometrial carcinoma cells , and suppressed ovarian cancer cell proliferation . MiR-504 stimulated cancer cell apoptosis and hindered cancer cell proliferation [67, 68]. MiR-92b promoted the growth of non-small-cell lung cancer cells  and controlled the G1/S cell cycle phase in human embryonic stem cells . In addition, miR-9 can act as a candidate tumor suppressor gene in recurrent ovarian cancer  and regulate neural development . Due to the important effects of the four miRNAs on ovarian function, their corresponding target genes were further screened and identified. Furthermore, previous studies showed that these target genes of the four miRNAs, which included Fas Associated Via Death Domain (FADD), Period Circadian Clock 2 (PER2), Cadherin EGF LAG Seven-Pass G-Type Receptor 1 (CELSR1), and EGFLike Domain Multiple 7 (EGFL7) were related to reproductive deficits , embryogenesis , and differentiation of embryonic stem cells .
The ceRNA network based on MRE prompted the further analysis of the coexpression relationship among these RNAs. Therefore, found miRNAs (miR-125b, miR504, miR-92b-5P, and miR-9) were selected that are related to ovarian function. Their target genes (FADD, EGFL7, PER2, and CELSR1), target lncRNAs (MSTRG62621, MSTRG114143, and MSTRG167556), and target circRNAs (circ000675 and circ13607) were used to construct a coexpression network using Pearson's correlation coefficient according to the differentially expressed RNA data (Figure 5(b) and Supplementary Database S5). Compared to the network based on MRE, the coexpression network showed more complex relationships among these RNAs.
Firstly, circ13607 was correlated to the miR-504 expression level although they did not share the same MRE. The coexpression relationship could be mediated by FADD, MSTRG.114143, and miR-92b-5p because miR-504 and miR-92b-5p can target MSTRG.114143 and FADD with MRE. Secondly, miR-92b-5p correlated to FADD, PER2, and CELSR1 expression levels; however, miR-92b-5p and PER2 did not share the same MRE, suggesting that other mediators might mediate the coexpression relationship between both. The coexpression relationship between PER2 and miR-92b5p were mediated in two ways: CELSR1, PER2, and circ13607 competed with miR-125b and CELSR1 and circ13607 competed with miR-92b-5p.
To verify the reliability of the obtained RNA-sequence data, the expression levels of all RNAs were verified by q-PCR and the results are shown in Figure 5(b). This comparison confirmed that the q-PCR results for these RNAs expression were in accordance with the RNA-sequence data (Supplementary Figure 3B-E).
4.1. Dynamic Changes of DNA Epigenetic Marks in Ovarian Aging. This study described a combined DNA methylation and transcriptome analysis of OP ovaries and YP ovaries and found progressive changes of age-differential methylation and gene expression landscapes during ovarian aging. DNA methylation modification is a key epigenetic mechanism involved in the vital processes of mammalian development, such as imprinting, transcriptional silencing, and X-chromosome inactivation [47, 75], as well as in the regulation of ovarian cancer and diseases such as PCOS and POI  and oocyte aging . This study first reports the DNA methylome profiles during natural ovarian aging and identified several differentially methylated genes involved in embryonic development, death effector domain binding, and cell apoptotic signaling pathways, which indicated that DNA methylation played a vital role during ovarian aging. Similarly, DNA methylation among other tissues was also found to be responsible for tissue aging. For example, in human aging brain tissue, DMR genes were linked to neurodevelopment-related terms, comprising the regulation of neurogenesis and cell projection organization ; however, DNA methylation in porcine aging skeletal muscle was reported to be involved in the protein degradation process and also responsible for muscular atrophy . Several cell death pathway genes such as cysteinyl aspartate proteases 10 (CASP10) and CFLAR, which played important roles in apoptosis , were found in the current DMRs data. This further indicated that oocyte apoptosis led to the depletion of ovarian reserves during ovarian aging .
4.2. Differentially Expressed Transcriptome Related to the Ovarian Aging Cycle. In other genetic systems, it has been reported that several RNA species play critical roles in the regulation of ovarian physiology, and these RNA species are involved in protein-coding RNAs [79-81] and noncoding RNAs such as miRNAs [12-14] and lncRNAs [15, 16]. The current study presents the global transcriptome expression profile of ovarian aging and identified a large number of differentially expressed RNAs. Genes associated with these differentially expressed RNAs were mainly involved in the ovarian aging cycle such as cell apoptosis, reproduction and fertilization process, embryonic development, ovarian cumulus expansion, and ovulation cycle.
Furthermore, several species of RNA could operate as ceRNA, which communicated with and regulated each other by using MREs as language and competed for binding to common miRNAs, which affected the stability of target genes or the translational activity . In this experiment, ceRNA networks were developed based on MRE and coexpression levels and the ceRNA networks were mainly involved in essential ovary developmental biological processes. For example, as a ceRNA, FADD and PER2 competed for binding to miR-92b-5p, thus affecting the expression level of target genes. MiR-92b-5p might promote lung cancer and non-small-cell lung cancer cell growth  and controlled the G1/S cell cycle phase in human ES cells ; furthermore, FADD and PER2 genes were reported to be related to reproductive deficits  and embryogenesis . Understanding the crosstalk between these novel RNAs will provide insights into the regulatory networks of genes with implications in ovarian aging.
4.3. Model Organism for Human Ovarian Aging or Disease Research. These conclusions promote the further research on pigs as model animal for ovarian aging or research of human diseases. Previous research identified genes that were differentially expressed in an age-dependent manner using microarray analysis of mouse oocytes, ovary, and ovarian surface epithelial cells as well as human oocytes [79, 82-86]. Although microarray analysis has been performed in porcine to identify the mechanisms involved in the muscle and brain aging [51,87], the aging process of ovaries has not been specifically investigated. Here, YP (180 days old) and OP (8 years old) were studied to search the DNA methylation and transcriptional expression changes in ovaries during the aging process. Although this study has its limitations, since whole ovaries were sampled, which contain multiple cell types [86, 88], hundreds of differentially methylated and expressed genes were identified to be related to ovarian aging, and these global DNA methylation and transcriptome expression profiles of pig ovaries of different ages constitute a useful resource.
In future, investigations of age-related differential DNA methylation and gene expression in individual cell types are required. Furthermore, this study only selected two age groups and limited samples. Studying pigs of additional successive ages and incorporating more samples are compulsory to further understand the differences in epigenetic modifications related to age, in addition to complicated mechanisms underlying the aging process. Moreover, in addition to providing novel evidence for biomedical studies, genomic/epigenomic searches on pigs can be helpful to discover the underlying molecular basis of the economic traits of pigs. Such knowledge can be used to progress the efficiency of artificial selection to produce more piglets.
Microarray analyses were performed to identify the differentially methylated and expressed genes in ovary tissues of OP and YP. Through GO enrichment analyses and the construction of ceRNA networks, the functions of differentially methylated and expressed genes, correlated pathways, and mutual regulatory relationships were analyzed between coding and noncoding genes. The obtained results contribute to the understanding of the aging process in ovaries and provide a basis for the future search for the molecular mechanisms of ovarian aging.
All data generated as part of this study are included in this published article and its supplementary information files. Requests for the raw data should be sent to the corresponding author.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
The authors would like to thank all students and teachers of the Zoology Lab, Swine Genetics and Breeding of Sichuan Agricultural University, China, for experimental conditions and advice. The authors also thank the Sichuan Weimu Modern Agricultural Science and Technology Co., Ltd. for managing and slaughtering the research flock. This study was supported by the National Key R&D Program of China (2018YFD0501204), Chengdu Local Good Varieties of Livestock and Poultry Resources Protection and Exploitation and Utilization of Construction Projects, and Sichuan Agricultural Science and Technology Transformation Fund Project-Breeding and Demonstration of Grain-Saving and Resistance Lean Pig Line (2018NZZJ017).
Table S1: primer sequences for q-PCR. Table S2: summary of sequence data and read-alignment statistics. Table S3: percentage of cytosine methylation after extraction. Table S4: association analysis of regulated genes between both ovarian development stages. Table S5: overlap between different methylation regions and differentially expressed genes. Figure S1: DNA methylation levels of CpG ChG and CHH on each chromosome. Figure S2: Pearson's correlation between DNA methylation level and chromosomal features. (A) Scatter plot and trend line (Pearson's correlation), indicating the correlation between the chromosome length and methylation level. (B) Scatter plot and trend line (Pearson's correlation), indicating the correlation between the chromosome GC content and methylation level. (C) Scatter plot and trend line (Pearson's correlation), indicating the correlation between the gene number in chromosome and methylation. (D) Scatter plot and trend line (Pearson's correlation), indicating the correlation between the chromosome CGI ratio and methylation level. (E) Scatter plot and trend line (Pearson's correlation), indicating the correlation between the chromosome repeat number and methylation level. Figure S3: verification of whole genome bisulfite and RNA sequencing data. (A) DNA methylation was verified via bisulfite sequencing PCR. The expression levels of mRNAs (B), miRNAs (C), lncRNAs (D), and circRNAs (E) were verified via q-PCR and adjusted by three endogenous control genes (porcine GAPDH, ACTB, and U6 snRNA). Figure S4: DNA methylation levels and RNA expression levels on each chromosome. Supplementary Database S1: summary and function enrichment of different methylation regions. Supplementary Database S2: identified mRNAs, miRNAs, lncRNAs, and circRNAs in all samples. Supplementary Database S3: differentially expressed mRNAs, miRNAs, lncRNAs, and circRNAs between both stages. Supplementary Database S4: GO analysis of differentially expressed RNAs and function overlap. Supplementary Database S5: targeted by MRE and coexpression of mRNAs, miRNAs, lncRNAs, and circRNAs. (Supplementary Materials)
 S. Titus, F. Li, R. Stobezki et al., "Impairment of BRCA1related DNA double strand break repair leads to ovarian aging in mice and humans," Science Translational Medicine, vol. 5, no. 172, p. 172ra21, 2013.
 M. J. Faddy, R. G. Gosden, A. Gougeon, S. J. Richardson, and J. F. Nelson, "Accelerated disappearance of ovarian follicles in mid-life: implications for forecasting menopause," Human Reproduction, vol. 7, no. 10, pp. 1342-1346, 1992.
 K. R. Hansen, N. S. Knowlton, A. C. Thyer, J. S. Charleston, M. R. Soules, and N. A. Klein, "A new model of reproductive aging: the decline in ovarian non-growing follicle number from birth to menopause," Human Reproduction, vol. 23, no. 3, pp. 699-708, 2008.
 J. S. Younis, "Ovarian aging: latest thoughts on assessment and management," Current Opinion in Obstetrics and Gynecology, vol. 23, no. 6, pp. 427-434, 2011.
 E. Pelosi, E. Simonsick, A. Forabosco, J. E. Garcia-Ortiz, and D. Schlessinger, "Dynamics of the ovarian reserve and impact of genetic and epidemiological factors on age of menopause," Biology of Reproduction, vol. 92, no. 5, p. 22, 2015.
 T. Laisk-Podar, C. M. Lindgren, M. Peters et al., "Ovarian physiology and GWAS: biobanks, biology, and beyond," Trends in Endocrinology & Metabolism, vol. 27, no. 7, pp. 516-528, 2016.
 H. Teede, A. Deeks, and L. Moran, "Polycystic ovary syndrome: a complex condition with psychological, reproductive and metabolic manifestations that impacts on health across the lifespan," BMC Medicine, vol. 8, no. 1, pp. 1-10, 2010.
 S. Guil and M. Esteller, "DNA methylomes, histone codes and miRNAs: tying it all together," The International Journal of Biochemistry & Cell Biology, vol. 41, no. 1, pp. 87-95, 2009.
 A. D. Goldberg, C. D. Allis, and E. Bernstein, "Epigenetics: a landscape takes shape," Cell, vol. 128, no. 4, pp. 635-638,2007.
 R. Calicchio, L. Doridot, F. Miralles, C. Mehats, and D. Vaiman, "DNA methylation, an epigenetic mode of gene expression regulation in reproductive science," Current Pharmaceutical Design, vol. 20, no. 11, pp. 1726-1750, 2014.
 Y. Qian, J. Tu, N. L. S. Tang et al., "Dynamic changes of DNA epigenetic marks in mouse oocytes during natural and accelerated aging," The International Journal of Biochemistry & Cell Biology, vol. 67, pp. 121-127, 2015.
 A. Diez-Fraile, T. Lammens, K. Tilleman et al., "Age-associated differential microRNA levels in human follicular fluid reveal pathways potentially determining fertility and success ofin vitrofertilization," Human Fertility, vol. 17, no. 2, pp. 90-98, 2014.
 J. M. Moreno, M. J. Nunez, A. Quinonero et al., "Follicular fluid and mural granulosa cells microRNA profiles vary in in vitro fertilization patients depending on their age and oocyte maturation stage," Fertility and Sterility, vol. 104, no. 4, pp. 1037.e1-1046.e1, 2015.
 A. Schneider, S. J. Matkovich, B. Victoria et al., "Changes of ovarian microRNA profile in long-living ames dwarf mice during aging," PLoS One, vol. 12, no. 1, Article ID e0169213, 2017.
 Q. Gao, H. Ren, M. Chen et al., "Long non-coding RNAs regulate effects of 1[beta]-crystallin B2 on mouse ovary development," Molecular Medicine Reports, vol. 14, no. 5, pp. 4223-4231, 2016.
 H. Wang, Z. Fu, C. Dai et al., "LncRNAs expression profiling in normal ovary, benign ovarian cyst and malignant epithelial ovarian cancer," Scientific Reports, vol. 6, no. 1, 2016.
 S. R. Milligan, A. A. Khan, and B. C. Thorne, "Delayed pseudopregnancy in the mouse," Reproduction, vol. 60, no. 1, pp. 49-51, 1980.
 G. P. Adams, J. Singh, and A. R. Baerwald, "Large animal models for the study of ovarian follicular dynamics in women," Theriogenology, vol. 78, no. 8, pp. 1733-1748, 2012.
 A. J. Ziecik, K. Derecka, B. Gawronska, A. Stepien, and G. Bodek, "Nongonadal LH/hCG receptors in pig: functional importance and parallels to human," Seminars in Reproductive Medicine, vol. 19, no. 1, pp. 19-30, 2001.
 M. J. Welsh, C. S. Rogers, D. A. Stoltz et al., "Development of a porcine model of cystic fibrosis," Transactions of the American Clinical & Climatological Association, vol. 120, pp. 149-162, 2009.
 A. E. Newell-Fugate, J. N. Taibl, S. G. Clark et al., "Effects of diet-induced obesity on metabolic parameters and reproductive function in female Ossabaw minipigs," Comparative Medicine, vol. 64, no. 1, pp. 44-49, 2014.
 A. E. Newell-Fugate, J. N. Taibl, M. Alloosh et al., "Effects of obesity and metabolic syndrome on steroidogenesis and folliculogenesis in the female ossabaw mini-pig," PLoS One, vol. 10, no. 6, Article ID e0128749, 2015.
 E. W. Demerath, C. T. Liu, N. Franceschini et al., "Genomewide association study of age at menarche in African-American women," Human Molecular Genetics, vol. 22, no. 16, pp. 3329-3346, 2013.
 J. R. Perry, F. Day, C. E. Elks et al., "Parent-of-origin-specific allelic associations among 106 genomic loci for age at menarche," Nature, vol. 514, no. 520, pp. 92-97, 2014.
 K. L. Lunetta, F. R. Day, P. Sulem et al., "Corrigendum: rare coding variants and X-linked loci associated with age at menarche," Nature Communications, vol. 6, no. 1, p. 10257, 2015.
 J. R. B. Perry, Y.-H. Hsu, D. I. Chasman et al., "DNA mismatch repair gene MSH6 implicated in determining age at natural menopause," Human Molecular Genetics, vol. 23, no. 9, pp. 2490-2497, 2014.
 J.-A. Pyun, S. Kim, N. H. Cho et al., "Genome-wide association studies and epistasis analyses of candidate genes related to age at menarche and age at natural menopause in a Korean population," Menopause, vol. 21, no. 5, pp. 522-529, 2014.
 F. R. Day, K. S. Ruth, D. J. Thompson et al., "Large-scale genomic analyses link reproductive ageing to hypothalamic signaling, breast cancer susceptibility and BRCA1-mediated DNA repair," Nature Genetics, vol. 47, no. 11, pp. 1294-1303, 2015.
 M. F. Rothschild and J. P. Bidanel, "Genetic aspects of domestication, common breeds and their origin,," in The Genetics of the Pig, M. F. Rothschild and A. Ruvinsky, Eds., pp. 313-344, CAB International, Wallingford, UK, 1998.
 Q. Li, X. Yuan, Z. Chen et al., "Heritability estimates and effect on lifetime reproductive performance of age at puberty in sows," Animal Reproduction Science, vol. 195, pp. 207-215, 2018.
 F. Krueger and S. R. Andrews, "Bismark: a flexible aligner and methylation caller for bisulfite-seq applications," Bioinformatics, vol. 27, no. 11, pp. 1571-1572, 2011.
 H. Li, B. Handsaker, A. Wysoker et al., "The sequence alignment/map format and SAMtools," Bioinformatics, vol. 25, no. 16, pp. 2078-2079, 2009.
 A. R. Quinlan and I. M. Hall, "BEDTools: a flexible suite of utilities for comparing genomic features," Bioinformatics, vol. 26, no. 6, pp. 841-842, 2010.
 H. Wu, T. Xu, H. Feng et al., "Detection of differentially methylated regions from whole-genome bisulfite sequencing data without replicates," Nucleic Acids Research, vol. 43, no. 21, p. e141, 2015.
 Y. Park and H. Wu, "Differential methylation analysis for BSseq data under general experimental design," Bioinformatics, vol. 32, no. 10, pp. 1446-1453, 2016.
 A. Alexa and R. Jorg, "topGO: Enrichment Analysis for Gene Ontology," R Package Version 2.0, 2010, http://bioconductor. uib.no/2.7/bioc/html/topGO.html.
 D. Kim, G. Pertea, C. Trapnell, H. Pimentel, R. Kelley, and S. L. Salzberg, "TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions," Genome Biology, vol. 14, no. 4, p. R36, 2013.
 M. Pertea, G. M. Pertea, C. M. Antonescu, T.-C. Chang, J. T. Mendell, and S. L. Salzberg, "StringTie enables improved reconstruction of a transcriptome from RNA-seq reads," Nature Biotechnology, vol. 33, no. 3, pp. 290-295, 2015.
 R. Patro, G. Duggal, M. I. Love, R. A. Irizarry, and C. Kingsford, "Salmon provides fast and bias-aware quantification of transcript expression," Nature Methods, vol. 14, no. 4, pp. 417-419, 2017.
 L. Kong, Y. Zhang, Z. Q. Ye et al., "CPC: assess the proteincoding potential of transcripts using sequence features and support vector machine," Nucleic Acids Research, vol. 35, no. 2, pp. 345-349, 2007.
 M. Punta, P. C. Coggill, R. Y. Eberhardt et al., "The Pfam protein families database," Nucleic Acids Research, vol. 40, no. D1, pp. D290-D301, 2012.
 M. R. Friedlander, S. D. Mackowiak, N. Li, W. Chen, and N. Rajewsky, "miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades," Nucleic Acids Research, vol. 40, no. 1, pp. 37-52, 2011.
 D. Betel, M. Wilson, A. Gabow, D. S. Marks, and C. Sander, "The microRNA.org resource: targets and expression," Nucleic Acids Research, vol. 36, pp. 149-153, 2008.
 M. D. Robinson, D. J. McCarthy, and G. K. Smyth, "edgeR: a bioconductor package for differential expression analysis of digital gene expression data," Bioinformatics, vol. 26, no. 1, pp. 139-140, 2010.
 K. J. Livak and T. D. Schmittgen, "Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method," Methods, vol. 25, no. 4, pp. 402-408, 2012.
 V. L. Wilson, R. A. Smith, S. Ma, and R. G. Cutler, "Genomic 5-methyldeoxycytidine decreases with age," Journal of Biological Chemistry, vol. 262, no. 2, pp. 9948-9951, 1987.
 S. Gkountela, K. X. Zhang, T. A. Shafiq et al., Cell, vol. 161, no. 6, pp. 1425-1436, 2015.
 F. Guo, L. Yan, H. Guo et al., "The transcriptome and DNA methylome landscapes of human primordial germ cells," Cell, vol. 161, no. 6, pp. 1437-1452, 2015.
 L. Shen, J. Du, Y. Xia et al., "Genome-wide landscape of DNA methylomes and their relationship with mRNA and miRNA transcriptomes in oxidative and glycolytic skeletal muscles," Scientific Reports, vol. 6, no. 1, 32186 pages, 2016.
 M. Li, H. Wu, Z. Luo et al., "An atlas of DNA methylomes in porcine adipose and muscle tissues," Nature Communications, vol. 3, no. 1, 2012.
 L. Jin, Z. Jiang, Y. Xia et al., "Genome-wide DNA methylation changes in skeletal muscle between young and middle-aged pigs," BMC Genomics, vol. 15, no. 1, p. 653, 2014.
 J. K. Kim, M. Samaranayake, and S. Pradhan, "Epigenetic mechanisms in mammals," Cellular and Molecular Life Sciences, vol. 66, no. 4, pp. 596-612, 2009.
 P.-X. Wang, Y.-X. Ji, X.-J. Zhang et al., "Targeting CASP8 and FADD-like apoptosis regulator ameliorates nonalcoholic steatohepatitis in mice and nonhuman primates," Nature Medicine, vol. 23, no. 4, pp. 439-449, 2017.
 C. Cui, B. Chatterjee, D. Francis et al., "Disruption of Mks1 localization to the mother centriole causes cilia defects and developmental malformations in Meckel-Gruber syndrome," Disease Models & Mechanisms, vol. 4, no. 1, pp. 43-56, 2011.
 P. A. Jones, "Functions of DNA methylation: islands, start sites, gene bodies and beyond," Nature Reviews Genetics, vol. 13, no. 7, pp. 484-492, 2012.
 M. P. Ball, J. B. Li, Y. Gao et al., "Targeted and genome-scale strategies reveal gene-body methylation signatures in human cells," Nature Biotechnology, vol. 27, no. 4, pp. 361-368, 2009.
 V. K. Rakyan, T. A. Down, N. P. Thorne et al., "An integrated resource for genome-wide identification and analysis of human tissue-specific differentially methylated regions (tDMRs)," Genome Research, vol. 18, no. 9, pp. 1518-1529, 2008.
 K. Lokk, V. Modhukur, B. Rajashekar et al., "DNA methylome profiling of human tissues identifies global and tissue-specific methylation patterns," Genome Biology, vol. 15, no. 4, p. r54, 2014.
 X. Zhang, D. Bogunovic, B. Payelle-Brogard et al., "Human intracellular ISG15 prevents interferon-a/S over-amplification and auto-inflammation," Nature, vol. 517, no. 7532, pp. 89-93, 2015.
 R. A. Veitia, "FOXL2 versus SOX9: a lifelong "battle of the sexes"," BioEssays, vol. 32, no. 5, pp. 375-380, 2010.
 A. M. Krichevsky, K.-C. Sonntag, O. Isacson, and K. S. Kosik, "Specific microRNAs modulate embryonic stem cell-derived neurogenesis," Stem Cells, vol. 24, no. 4, pp. 857-864, 2006.
 A. Laios, S. O'Toole, R. Flavin et al., "Potential role of miR-9 and miR-223 in recurrent ovarian cancer," Molecular Cancer, vol. 7, no. 1, p. 35, 2008.
 L. Salmena, L. Poliseno, Y. Tay, L. Kats, and P. P. Pandolfi, "A ceRNA hypothesis: the rosetta stone of a hidden RNA language?," Cell, vol. 146, no. 3, pp. 353-358, 2011.
 M. T. N. Le, C. Teh, N. Shyh-Chang et al., "MicroRNA-125b is a novel negative regulator of p53," Genes & Development, vol. 23, no. 7, pp. 862-876, 2009.
 F. Jiang, T. Liu, Y. He et al., "MiR-125b promotes proliferation and migration of type II endometrial carcinoma cells through targeting TP53INP1 tumor suppressor in vitro and in vivo," BMC Cancer, vol. 11, no. 1, 425 pages, 2011.
 Y. Guan, H. Yao, Z. Zheng, G. Qiu, and K. Sun, "MiR-125b targets BCL3 and suppresses ovarian cancer proliferation," International Journal of Cancer, vol. 128, no. 10, pp. 2274-2283, 2011.
 N. Kikkawa, T. Kinoshita, N. Nohata et al., "microRNA-504 inhibits cancer cell proliferation via targeting CDK6 in hypopharyngeal squamous cell carcinoma," International Journal of Oncology, vol. 44, no. 6, pp. 2085-2092, 2014.
 R. Cui, Y. Guan, C. Sun et al., "A tumor-suppressive microRNA, miR-504, inhibits cell proliferation and promotes apoptosis by targeting FOXP1 in human glioma," Cancer Letters, vol. 374, no. 1, pp. 1-11, 2016.
 L. Lei, Y. Huang, and W. Gong, "Inhibition of miR-92b suppresses nonsmall cell lung cancer cells growthand motility by targeting RECK," Molecular & Cellular Biochemistry, vol. 387, no. 1-2, pp. 171-176, 2014.
 S. Sengupta, J. Nie, R. J. Wagner, C. Yang, R. Stewart, and J. A. Thomson, "MicroRNA 92b controls the G1/S checkpoint gene p57 in human embryonic stem cells," Stem Cells, vol. 27, no. 7, pp. 1524-1528, 2009.
 M. Coolen, S. Katz, and L. Bally-Cuif, "miR-9: a versatile regulator of neurogenesis," Frontiers in Cellular Neuroscience, vol. 7220 pages, 2013.
 V. Pilorz and S. Steinlechner, "Low reproductive success in Per1 and Per2 mutant mouse females due to accelerated ageing?," Reproduction, vol. 135, no. 4, pp. 559-568, 2008.
 H. Z. Imtiyaz, X. Zhou, H. Zhang, D. Chen, T. Hu, and J. Zhang, "The death domain of FADD is essential for embryogenesis, lymphocyte development, and proliferation," Journal of Biological Chemistry, vol. 284, no. 15, pp. 9917-9926, 2009.
 A. Durrans and H. Stuhlmann, "A role for Egfl7 during endothelial organization in the embryoid body model system," Journal of Angiogenesis Research, vol. 2, no. 1, p. 4, 2010.
 H. Sasaki and Y. Matsui, "Epigenetic events in mammalian germ-cell development: reprogramming and beyond," Nature Reviews Genetics, vol. 9, no. 2, pp. 129-140, 2008.
 G. Oh, S. Ebrahimi, S. C. Wang et al., "Epigenetic assimilation in the aging human brain," Genome Biology, vol. 17, no. 1, p. 76, 2016.
 F. Liu, F. Li, L. Luo et al., "Genetic variants in cell death pathway genes and HBV-related hepatocellular carcinoma among a Chinese Han population," Apoptosis, vol. 22, no. 8, pp. 1035-1047, 2017.
 M. Tiwari, S. Prasad, A. Tripathi et al., "Apoptosis in mammalian oocytes: a review," Apoptosis, vol. 20, no. 8, pp. 1019-1025, 2015.
 M. L. Gr0ndahl, C. Yding Andersen, J. Bogstad, F. C. Nielsen, H. Meinertz, and R. Borup, "Gene expression profiles of single human mature oocytes in relation to age," Human Reproduction, vol. 25, no. 4, pp. 957-968, 2010.
 H. Wei, X. Liu, J. Yuan et al., "Age-specific gene expression profiles of rhesus monkey ovaries detected by microarray analysis," BioMed Research International, vol. 2015, Article ID 625192, 15 pages, 2015.
 V. Govindaraj, H. Krishnagiri, P. Chakraborty, M. Vasudevan, and A. J. Rao, "Age-related changes in gene expression patterns of immature and aged rat primordial follicles," Systems Biology in Reproductive Medicine, vol. 63, no. 1, pp. 37-48, 2017.
 T. Hamatani, G. Falco, M. G. Carter et al., "Age-associated alteration of gene expression patterns in mouse oocytes," Human Molecular Genetics, vol. 13, no. 19, pp. 2263-2278, 2004.
 A. Zimon, A. Erat, T. Von Wald et al., "Genes invoked in the ovarian transition to menopause," Nucleic Acids Research, vol. 34, no. 11, pp. 3279-3287, 2006.
 N. M. Steuerwald, M. G. Bermudez, D. Wells, S. Munne, and J. Cohen, "Maternal age-related differential global expression profiles observed in human oocytes," Reproductive BioMedicine Online, vol. 14, no. 6, pp. 700-708, 2007.
 A. A. Sharov, G. Falco, Y. Piao et al., "Effects of aging and calorie restriction on the global gene expression profiles of mouse testis and ovary," BMC Biology, vol. 6, no. 1, p. 24, 2008.
 N. Gava, C. L. Clarke, C. Bye, K. Byth, and A. deFazio, "Global gene expression profiles of ovarian surface epithelial cells in vivo," Journal of Molecular Endocrinology, vol. 40, no. 6, pp. 281-296, 2008.
 M. T. Ven0, T. B. Hansen, S. T. Ven0 et al., "Spatio-temporal regulation of circular RNA expression during porcine embryonic brain development," Genome Biology, vol. 16, no. 1, p. 245, 2015.
 A. Uyar, S. Torrealday, and E. Seli, "Cumulus and granulosa cell markers of oocyte and embryo quality," Fertility and Sterility, vol. 99, no. 4, pp. 979-997, 2013.
Xiaoyu Xi, (1) Qin Zou, (1) Yingying Wei, (1) Yan Chen, (1) Xue Wang, (1) Daojun Lv, (2) Peilin Li, (2) Anxiang Wen, (1) Li Zhu [ID], (3) Guoqing Tang, (3) Jideng Ma, (3) Mingzhou Li [ID], (3) Xuewei Li [ID], (3) and Yanzhi Jiang [ID] (1)
(1) Department of Zoology, College of Life Science, Sichuan Agricultural University, Ya'an, Sichuan 625014, China
(2) Sichuan Weimu Modern Agricultural Science and Technology Co., Ltd., Chengdu, Sichuan 611536, China
(3) Institute of Animal Genetics and Breeding, College of Animal Science and Technology, Sichuan Agricultural University, Chengdu, Sichuan 611130, China
Correspondence should be addressed to Yanzhi Jiang; firstname.lastname@example.org
Xiaoyu Xi, Qin Zou, and Yingying Wei contributed equally to this work.
Received 15 February 2019; Revised 25 June 2019; Accepted 30 September 2019; Published 30 October 2019
Academic Editor: Thomas Lufkin
Caption: Figure 1: DNA methylation associated with ovarian aging. (a) DNA methylation levels between both ovarian development stages. YP-1and YP-2 represent sample 1 and sample 2 from young sows at the puberty stage, respectively; OP-1 and OP-2 represent sample 1 and sample 2 from old sows at the reproductive exhaustion stage, respectively. (b) Distribution of CG methylation reads on and around the gene body region. Abbreviations: TSS, transcription start site; TTS, transcription termination site. Gene ontology (GO) function enrichment of the (c) hypermethylated genes and (d) hypomethylated genes related to the different methylation regions (DMRs).
Caption: Figure 2: Differentially expressed transcriptome between both ovarian development stages. The x-axis indicates log2FC, and the y-axis indicates the-[log.sub.10] P value. Criteria of [absolute value of FC] < 1 and P value <0.05 were used to screen differently expressed RNAs. (a) Differentially expressed mRNAs between both ovarian development stages. (b) Differentially expressed miRNAs between both ovarian developmental stages. (c) Differentially expressed lncRNAs between both ovarian developmental stages. (d) Differentially expressed circRNAs between both ovarian developmental stages.
Caption: Figure 3: Gene ontology (GO) function enrichment of differentially expressed RNAs between both ovarian development stages. GO analysis of differentially expressed (a) mRNAs, (b) miRNA target genes, (c) lncRNA target genes, and (d) circRNA source genes.
Caption: Figure 4: Overlaps of gene ontology (GO) function enrichment among different DNA methylation regions and different expression RNAs. The top different methylation regions (DMRs) and expressed RNAs (miRNAs, mRNAs, lncRNAs, and circRNAs) were both GO enriched. The x-axis indicates the RNA class and DMR enriched at GO term, the upper half of y-axis indicates the P value of enrichment, and the lower half of y-axis indicates the ratio of differentially expressed and background of GO term.
Caption: Figure 5: Competing endogenous RNA (ceRNA) network of ovarian aging. (a) The ceRNA network was based on miRNA- mRNA, miRNA-lncRNA, and miRNA-circRNA interactions with microRNA response elements (MREs). (b) Coexpression network between four RNAs classes. Coexpressed RNAs pairs were identified using strict screening criteria (Pearson's correlation coefficients >0.85 or <-0.85, P < 0.01).
|Printer friendly Cite/link Email Feedback|
|Title Annotation:||Research Article|
|Author:||Xi, Xiaoyu; Zou, Qin; Wei, Yingying; Chen, Yan; Wang, Xue; Lv, Daojun; Li, Peilin; Wen, Anxiang; Zhu|
|Publication:||BioMed Research International|
|Date:||Nov 30, 2019|
|Previous Article:||The Potential of Hibiscus sabdariffa Linn in Inducing Glucagon-Like Peptide-1 via SGLT-1 and GLPR in DM Rats.|
|Next Article:||Curcumin Inhibits ERK/c-Jun Expressions and Phosphorylation against Endometrial Carcinoma.|