Dynamic Changes of DNA Methylation and Transcriptome Expression in Porcine Ovaries during Aging

Xi, Xiaoyu; Zou, Qin; Wei, Yingying; Chen, Yan; Wang, Xue; Lv, Daojun; Li, Peilin; Wen, Anxiang; Zhu, Li; Tang, Guoqing; Ma, Jideng; Li, Mingzhou; Li, Xuewei; Jiang, Yanzhi

doi:https://doi.org/10.1155/2019/8732023

BioMed Research International

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusions Data Availability Conflicts of Interest Acknowledgments Supplementary Materials References Copyright Related Articles

Research Article | Open Access

Volume 2019 | Article ID 8732023 | https://doi.org/10.1155/2019/8732023

Dynamic Changes of DNA Methylation and Transcriptome Expression in Porcine Ovaries during Aging

Xiaoyu Xi,¹Qin Zou,¹Yingying Wei,¹Yan Chen,¹Xue Wang,¹Daojun Lv,²Peilin Li,²Anxiang Wen,¹Li Zhu,³Guoqing Tang,³Jideng Ma,³and Mingzhou Li³ et al.

Academic Editor: Thomas Lufkin

Received15 Feb 2019

Revised25 Jun 2019

Accepted30 Sept 2019

Published30 Oct 2019

Abstract

The biological function of human ovaries declines along with aging. To identify the underlying molecular changes during ovarian aging, pigs were used as model animals. Genome-wide DNA methylation and transcriptome-wide RNA expression analyses were performed via high-throughput sequencing of ovaries from young pigs (180 days, puberty stage of first ovulation) and old pigs (eight years, reproductive exhaustion stage). The results identified 422 different methylation regions between old and young pigs; furthermore, a total of 2,243 mRNAs, 95 microRNAs, 248 long noncoding RNAs (lncRNAs), and 116 circular RNAs (circRNAs) were differentially expressed during both developmental stages. Gene ontology analysis showed that these genes related to different methylation and expression are involved in the ovarian aging cycle. Specifically, these are involved in cell apoptosis, death effector domain binding, embryonic development, reproduction and fertilization process, ovarian cumulus expansion, and the ovulation cycle. Multigroup cooperative control relationships were also assessed, and competing endogenous RNA (ceRNA) networks were constructed in the ovarian aging cycle. These data will help to clarify ovary age-associated potential molecular changes in DNA methylation and transcriptional patterns over time.

1. Introduction

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 [1]. 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 [2]. 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 [9]. “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 [17]. 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 [18]. 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 [19] as well as their anatomical, physiological, and biochemical similarities to humans [20]. 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 [3]. 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 high-throughput 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 ± 1 days. An initial diet started from the 30^th to the 60^th day after weaning and contained 3.40 Mcal·kg⁻¹ of metabolizable energy with 20.0% crude protein (11.5 g/kg lysine). From the 61^st to the 120^th 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 121^st 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°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°C in a thermal cycler and then for 60 min at 54°C. Bisulfite-treated DNA was purified via spin column and was used to build a sequencing library with the help of the EpiGnome™ 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 2^nd 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 [31]. 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 [32]. 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–4 bp in the read [31]. After reads were alimented, SAMtools were applied to deal with the bam file out form Bowtie2 [32]. CpG island locus information was downloaded from the UCSC genome browser, and CpG island methylation level of each sample was calculated with bedtools [33]. 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 [36]. 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 μ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® Ultra™ Directional RNA Library Prep Kit for Illumina® (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 [37].

StringTie software was used to assemble the mapped reads per sample [38], 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 [39]. 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 μg RNA of each sample was used as input material for library preparation of small RNA. Sequencing libraries were created using NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (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 [42]. The miRanda (v3.3a, default parameters) and cutoffs (score S ≥ 140 and energy E ≤ −20.0) were used to predict the miRNA target [43]. Then, miRNA expression was assessed by TPMs through the following criteria: normalization formula: normalized expression = 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 [44]. A differential expression value <0.05 and fold change >2 were assigned as differentially expressed in different comparisons. GO enrichment was performed by TopGO [36]. 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 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 Real-Time PCR Detection system (Bio-Rad, Hercules, CA, USA) following the manufacturer’s directions. Three endogenous control genes, glyceraldehyde-3-phosphate dehydrogenase (GAPDH), β-actin (ACTB), and small nuclear RNA (U6 snRNA), were used in this assay. The 2^−ΔΔCt method was used to determine the expression levels of objective mRNAs, miRNAs, lncRNAs, and circRNAs [45]. These primers are shown in Supplementary Table S1.

3. Results

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 genome-wide loss of DNA methylation was found in response to age-related epigenetic modification [46]. 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).

(a)

(b)

(c)

(d)

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).

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 [49].

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, = 1.844 × 10⁻⁹), and a positive correlation was observed with the GC content (the percentage of guanine and cytosine, r = 0.903, < 2.2 × 10⁻¹⁶), gene number (being calculated for a 1 Mb window in chromosomes, r = 0.398, < 2.2 × 10⁻¹⁶), CpG_o/e (the ratio between the observed and expected numbers of CpG sites, r = 0.792, < 2.2 × 10⁻¹⁶), and repeat region density (being calculated for a 1 Mb window in chromosomes, r = 0.126, < 2.2 × 10⁻¹⁶) (Supplementary Figure S2). These results are similar with previous reports [50, 51]. Among these genomic features, the GC content and CpG_o/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 [52].

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 (, 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 [51].

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 cysteine-type 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 [53] 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 [54] 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 [55]. 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, = 5.72 × 10⁻⁷) 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 [51]. However, previous scientists reported a positive correlation with gene expression levels [56, 57] or no clear relationship patterns [58].

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 (), 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-α/β overamplification and autoinflammation [59], 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 [60]. The higher expression level of SOX9 in OP indicated that the reproductive performance decline in old females.

(a)

(b)

(c)

(d)

Figure 2

Differentially expressed transcriptome between both ovarian development stages. The x-axis indicates log₂FC, and the y-axis indicates the−log₁₀ value. Criteria of |FC| < 1 and 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.

(a)

(b)

(c)

(d)

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 [61] and has prospective importance in recurrent ovarian cancer as biomarkers [62].

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 [63], 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-9-1, 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 p53-induced apoptosis [64], promoted both the proliferation and migration of type II endometrial carcinoma cells [65], and suppressed ovarian cancer cell proliferation [66]. MiR-504 stimulated cancer cell apoptosis and hindered cancer cell proliferation [67, 68]. MiR-92 b promoted the growth of non-small-cell lung cancer cells [69] and controlled the G1/S cell cycle phase in human embryonic stem cells [70]. In addition, miR-9 can act as a candidate tumor suppressor gene in recurrent ovarian cancer [62] and regulate neural development [71]. 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 EGF-Like Domain Multiple 7 (EGFL7) were related to reproductive deficits [72], embryogenesis [73], and differentiation of embryonic stem cells [74].

(a)

(b)

The ceRNA network based on MRE prompted the further analysis of the coexpression relationship among these RNAs. Therefore, found miRNAs (miR-125b, miR-504, 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-92b-5p 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. Discussion

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 [10] and oocyte aging [11]. 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 [76]; however, DNA methylation in porcine aging skeletal muscle was reported to be involved in the protein degradation process and also responsible for muscular atrophy [51]. Several cell death pathway genes such as cysteinyl aspartate proteases 10 (CASP10) and CFLAR, which played important roles in apoptosis [77], were found in the current DMRs data. This further indicated that oocyte apoptosis led to the depletion of ovarian reserves during ovarian aging [78].

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 [63]. 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 [69] and controlled the G1/S cell cycle phase in human ES cells [70]; furthermore, FADD and PER2 genes were reported to be related to reproductive deficits [72] and embryogenesis [73]. 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.

5. Conclusions

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.

Data Availability

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.

Acknowledgments

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).

Supplementary Materials

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)

References

S. Titus, F. Li, R. Stobezki et al., “Impairment of BRCA1-related DNA double strand break repair leads to ovarian aging in mice and humans,” Science Translational Medicine, vol. 5, no. 172, p. 172ra21, 2013.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
A. D. Goldberg, C. D. Allis, and E. Bernstein, “Epigenetics: a landscape takes shape,” Cell, vol. 128, no. 4, pp. 635–638, 2007.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
J. M. Moreno, M. J. Núñez, A. Quiñonero 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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
Q. Gao, H. Ren, M. Chen et al., “Long non-coding RNAs regulate effects of β-crystallin B2 on mouse ovary development,” Molecular Medicine Reports, vol. 14, no. 5, pp. 4223–4231, 2016.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
S. R. Milligan, A. A. Khan, and B. C. Thorne, “Delayed pseudopregnancy in the mouse,” Reproduction, vol. 60, no. 1, pp. 49–51, 1980.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Google Scholar
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.
View at: Google Scholar
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.
View at: Publisher Site | Google Scholar
E. W. Demerath, C. T. Liu, N. Franceschini et al., “Genome-wide association study of age at menarche in African-American women,” Human Molecular Genetics, vol. 22, no. 16, pp. 3329–3346, 2013.
View at: Publisher Site | Google Scholar
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.
View at: Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
H. Li, B. Handsaker, A. Wysoker et al., “The sequence alignment/map format and SAMtools,” Bioinformatics, vol. 25, no. 16, pp. 2078-2079, 2009.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
Y. Park and H. Wu, “Differential methylation analysis for BS-seq data under general experimental design,” Bioinformatics, vol. 32, no. 10, pp. 1446–1453, 2016.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
L. Kong, Y. Zhang, Z. Q. Ye et al., “CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine,” Nucleic Acids Research, vol. 35, no. 2, pp. 345–349, 2007.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
M. R. Friedländer, 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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Google Scholar
S. Gkountela, K. X. Zhang, T. A. Shafiq et al., Cell, vol. 161, no. 6, pp. 1425–1436, 2015.
View at: Publisher Site
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
X. Zhang, D. Bogunovic, B. Payelle-Brogard et al., “Human intracellular ISG15 prevents interferon-α/β over-amplification and auto-inflammation,” Nature, vol. 517, no. 7532, pp. 89–93, 2015.
View at: Publisher Site | Google Scholar
R. A. Veitia, “FOXL2 versus SOX9: a lifelong “battle of the sexes”,” BioEssays, vol. 32, no. 5, pp. 375–380, 2010.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
L. Lei, Y. Huang, and W. Gong, “Inhibition of miR-92b suppresses nonsmall cell lung cancer cells growth and motility by targeting RECK,” Molecular & Cellular Biochemistry, vol. 387, no. 1-2, pp. 171–176, 2014.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
M. Coolen, S. Katz, and L. Bally-Cuif, “miR-9: a versatile regulator of neurogenesis,” Frontiers in Cellular Neuroscience, vol. 7, 220 pages, 2013.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
M. Tiwari, S. Prasad, A. Tripathi et al., “Apoptosis in mammalian oocytes: a review,” Apoptosis, vol. 20, no. 8, pp. 1019–1025, 2015.
View at: Publisher Site | Google Scholar
M. L. Grøndahl, 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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
N. M. Steuerwald, M. G. Bermúdez, D. Wells, S. Munné, 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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
M. T. Venø, T. B. Hansen, S. T. Venø et al., “Spatio-temporal regulation of circular RNA expression during porcine embryonic brain development,” Genome Biology, vol. 16, no. 1, p. 245, 2015.
View at: Google Scholar
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.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2019 Xiaoyu Xi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1542

Downloads

1518

Citations