Temporal transcriptomic landscape of postnatal mouse ovaries reveals dynamic gene signatures associated with ovarian aging.

The ovary is the most important organ for maintaining female reproductive health, but it fails before most other organs. Aging-associated alterations in gene expression patterns in mammalian ovaries remain largely unknown. In this study, the transcriptomic landscape of postnatal mouse ovaries over the reproductive lifespan was investigated using bulk RNA sequencing in C57BL/6 mice. Gene expression dynamics revealed that the lifespan of postnatal mouse ovaries comprised four sequential stages, during which 2517 genes were identified as differentially enriched. Notably, the DNA repair pathway was found to make a considerable and specific contribution to the process of ovarian aging. Temporal gene expression patterns were dissected to identify differences in gene expression trajectories over the lifespan. In addition to DNA repair, distinct biological functions (including hypoxia response, epigenetic modification, fertilization, mitochondrial function, etc.) were overrepresented in particular clusters. Association studies were further performed to explore the relationships between known genes responsible for ovarian function and differentially expressed genes identified in this work. We found that the causative genes of human premature ovarian insufficiency were specifically enriched in distinct gene clusters. Taken together, our findings reveal a comprehensive transcriptomic landscape of the mouse ovary over the lifespan, providing insights into the molecular mechanisms underlying mammalian ovarian aging and supporting future etiological studies of aging-associated ovarian disorders.

Introduction

The ovary is critical for reproductive success in female mammals and has two prominent functions: oogenesis and hormonogenesis (1). The ovary shelters the limited lifetime supply of oocytes, which can reach up to approximately two million in quantity by birth and reside in quiescent primordial follicles until puberty in humans (2). Additionally, ovaries cyclically secrete estrogen and progesterone in response to follicle-stimulating hormone and luteinizing hormone (LH) to maintain the menstrual cycle, which comprises follicular development and ovulation. There is a consensus that ovaries suffer from aging-associated dysfunctions chronologically before most of the other organs (3). Female fecundity declines from 25 to 30 years old in humans due to early-onset decreases in oocyte quantity and quality (4–6), whereas aging-related diseases derived from other organs, such as rheumatoid arthritis and chronic obstructive pulmonary disease, rarely develop before 60–65 years old (7). Moreover, ovaries exhibit higher sensitivity to environmental factors than most other organs. Prolonged dormancy could predispose oocytes to the accumulation of various types of damage, such as DNA double-strand breaks (DSBs), leading to the accelerated death of oocytes and subfertility of females during aging (8).

To guarantee offspring health, follicular development and ovulation are orchestrated by intricate molecular mechanisms and endocrine regulation. Once the physiological process is disturbed, aging-related dysfunctions might arise, such as early menopause with the final menstrual period (FMP) between the ages of 40 and 45 years and premature ovarian insufficiency (POI) with the FMP before 40 years old (9). Approximately 7.6% of females are affected by early menopause in populations of the UK, Scandinavia, Australia and Japan (10), while approximately 1% of females are diagnosed with POI globally (11,12). These females experience reduced fertility and menopausal symptoms with high susceptibility to chronic disorders such as cardiovascular disease, osteoporosis or type 2 diabetes, leading to an increased risk of early death (13–16). Thus, predictive or preventive approaches for ovarian aging are crucial to facilitate female fertility and maintain reproductive health.

Genetic factors are considered to account for up to 90% of phenotypic variation in menopausal age (17). However, the mechanisms underlying the regulation of reproductive lifespan in females remain largely unknown. Several genome-wide association studies have been carried out in humans to identify genetic loci associated with menopausal age in recent decades (18–20). Microarrays have also been utilized to pinpoint changed epigenetic or transcriptomic characteristics of oocytes and granulosa cells from subfertile or elderly females (21–23), but intact human ovaries are usually unavailable. From this perspective, Mus musculus has become a classic model. Comparisons between young and old mice have shown upregulation of oxidative stress responses and immune defense factors and downregulation of mitochondrial functions in old ovaries (24–26). However, studies on ovarian aging lack gene expression profiling with higher temporal resolution over the reproductive lifespan.

To this end, ovaries collected over the whole reproductive lifespan of mice (from neonatal, 3-week-old, 2-month-old, 3-month-old, 6-month-old, 9-month-old and 1-year-old mice) were utilized as a source of bulk transcriptomes from early development to aging. We found that the gene expression alterations of DNA repair genes make a considerable contribution to progressive mouse ovarian aging. Furthermore, distinct temporal gene expression patterns were identified, each of which possessed significant biological functions. These data provide potential biomarkers for different stages of ovarian aging. Additional association studies revealed that human POI causative genes were specifically enriched in clusters related to the meiotic cell cycle, DNA repair and the cellular response to hypoxia, implying the potential contribution of distinct gene clusters to the timing of ovarian aging.

Results

Gradually decreased follicle reserve in aging ovaries from C57BL/6 mice

Ovaries were collected from C57BL/6 mice at 2 weeks, 1 month, 2 months and 6 months of age. Histological methods showed a dynamic panorama of the process of aging, including follicle reserve, activation, development, maturation, ovulation and atresia (Fig. 1A). In infants (2 weeks old), more than 300 primordial follicles constituted the ovarian reserve, while approximately 40 primary and secondary follicles and no antral follicles were observed. Before puberty (1 month old), several well-developed antral follicles were detected, preparing for the first ovulation. When undergoing puberty (2 months old), the number of antral follicles continued to increase, and the formation of the corpus lutea indicated that the mature cumulus–oocyte complex (COC) had been successfully ovulated. At 6 months old, when the mice were reaching the period of premenopause or perimenopause, although the numbers of antral follicles and corpus lutea remained high, the dwindling primordial follicles indicated the depletion of the ovarian reserve (Fig. 1B). In summary, the mouse ovaries exhibited well-organized aging processes, and the data collected at 6 months of age foreshadowed the end of the reproductive period of female mice.

Figure 1

Follicle development and atresia in mouse ovaries during aging. (A) Representative images of H&E-stained ovaries from 2-week-old, 1-month-old, 2-month-old and 6-month-old mice. Scale bars represent 100 μm. (B) The numbers of primordial follicles, primary follicles, secondary follicles, antral follicles, atretic follicles and corpus lutea in ovaries from 2-week-old, 1-month-old, 2-month-old and 6-month-old mice were counted. ‘×’ represents the average, horizontal lines represent the median, the upper and lower edges of the box represent the upper and lower quartiles, the upper and lower bars represent the maximum and minimum values and the dot represents the outlier. N = 6. All data were analyzed by one-way ANOVA followed by the least significant difference method with Benjamini–Hochberg adjustment. Different letters indicate values that are significantly different; P < 0.05.

Regularly changed transcriptomic dynamics of mouse ovaries during reproductive development and aging

To better explore the dynamic changes in global gene expression levels throughout the reproductive lifespan, seven groups of ovaries from newborn mice and mice at 3 weeks (before the onset of puberty), 2 months, 3 months (the end of puberty when the mice achieve full sexual maturity), 6 months, 9 months (the endocrine equivalent of human perimenopause) (27) and 1 year (reproductive aging) of age (28) were harvested and subjected to bulk RNA sequencing (Fig. 2A).

Figure 2

Global transcriptome patterns of mouse ovaries over the reproductive lifespan. (A) Schematic diagram of the mouse ovarian RNA-seq workflow. The four points above the timeline represent the beginning of each physiological period of mice. The seven points below the timeline represent the age of the mice when the ovaries were harvested and subjected to RNA-seq. (B) PCA of RNA-seq data. PC1 segregated three reproductive periods (neonatal, prepubertal and sexual maturity). PC2 separated aging from puberty. Similar colors indicate the same stage. (C) Expression patterns of typical genes contributing to PC1 or PC2 exhibited on PCA plots. The gradient of gray to red represents low to high gene expression levels.

RNA from the ovaries was sequenced to 100× depth using an Illumina HiSeq platform, with an average of 40.1 million high-quality paired-end reads per sample (Supplementary Material, Table S1). The filtered reads were aligned to the reference genome database of M. musculus (GRCm38/mm10) by TopHat2, with a mapping ratio of approximately 90%. Identification of differentially expressed genes (DEGs) revealed 812, 391, 14, 118, 30 and 41 upregulated DEGs and 1038, 356, 4, 49, 74 and 12 downregulated DEGs between sequential groups [false discovery rate (FDR) < 0.05; fold change (FC) of normalized read counts >2], indicating drastic transcriptomic changes before puberty, as expected (Supplementary Material, Fig. S1). Notably, these transcriptional alterations decreased drastically during sexual maturation (3 versus 2 months) but clearly arose again during aging, with an average of 108 DEGs between adjacent time points after the end of puberty. The number of alternative splicing (AS) events displayed a trajectory similar to that of the DEGs (Supplementary Material, Fig. S2A), which might be partially explained by the transcription levels of splicing factors (Supplementary Material, Fig. S2B). Furthermore, the Pearson correlation coefficient (PCC) revealed the consistency between AS events and the changing expression of multiple splicing factors between adjacent groups, with Pearson’s r = 0.975 and P-value <0.001 (Supplementary Material, Fig. S2B).

Transcriptomic atlas defined four stages during the mouse ovarian lifespan

To explore the global gene expression profile of mouse ovaries during the whole lifespan, we conducted principal component analysis (PCA) on RNA sequencing data from 7 time points (Fig. 2B), and the results were consistent with the correlations between samples (Supplementary Material, Fig. S3). The proportion of the variance explained by the top three principal components (PC1–3) accounted for over 60% in total (Supplementary Material, Fig. S4). Specifically, PC1 and PC2 were the most important components for stage discrimination (Supplementary Material, Fig. S5). As shown in Fig. 2B, PC1 distinguished early developmental groups (0 day and 3 weeks) from adult groups, and PC2 further separated the adult groups by sexual maturation (2 and 3 months) and the aging process (6 months, 9 months and 1 year). All groups from the 7 time points finally fell into four stages, which were arranged according to age in a counterclockwise direction on the plot. According to the known developmental stages of the mouse reproductive system (27), we annotated the four stages as the neonatal stage (NN) containing the 0-day group, the prepubertal stage (PPB) containing the 3-week group, the pubertal stage (PB) containing the 2- and 3-month groups and the sexually mature or aged stage (MA) containing the 6-month, 9-month and 1-year groups.

Stage-specific enriched genes were identified against all the other stages (FDR < 0.05, FC of normalized read counts >2; Supplementary Material, Fig. S6A). There were 1538, 471, 56 and 452 stage-specific enriched genes recognized in NN, PPB, PB and MA, respectively. Among them, the gag-related retrotransposon gene Rtl3/Zcchc5 and genes involved in cell growth signaling pathways, such as Igf2 and Cdkn1c, were almost exclusively expressed in NN (Supplementary Material, Fig. S6B). The expression of Nlrp5, a maternal effector gene encoding a member of the subcortical maternal complex (SCMC); the ovarian primitive germ cell tumor marker Sall4 and several genes associated with microtubules, such as Kif17, were highly upregulated in PPB (Supplementary Material, Fig. S6B). Nnt, encoding an integral protein of the inner mitochondrial membrane that translocates protons, was typical of the few genes highly expressed in PB (Supplementary Material, Fig. S6B). The genes with elevated expression specifically in aging ovaries (MA) included Ly6a and Cd5l, which were associated with the inflammatory response, and Adcyap1/Pacap, encoding a secreted proprotein induced by LH and gonadotropins in granulosa cells of preovulatory follicles (Supplementary Material, Fig. S6B). Altogether, clustering of groups based on the transcriptomic atlas well defined four reproductive developmental stages, consistent with the physiological process. Furthermore, some of these stage-specific enriched genes could be potential biomarkers for ovarian development and aging, e.g. Nlrp5, Adcyap1 and Cd5l.

DNA repair pathways made considerable contributions to progressive mouse ovarian aging

Next, the biological significance of PC1 and PC2 was elucidated. A total of 1395 and 433 genes significantly contributing to PC1 and PC2, respectively (P-value <0.05, PCC > 0.85), were recognized. Representatives from the top 50 contributors are shown in Fig. 2C. The expression of Ptk7 (contribution = 0.0220%, the same below), which is essential for the Wnt signaling pathway, and Tbx2 (0.0217%), which is associated with the regulation of mesoderm differentiation, was elevated mainly in early developmental stages and downregulated upon PPB, while the expression of Serpina5 (0.0217%), which encodes a heparin-dependent serine protease inhibitor that acts in the reproductive tract, was induced across PC1 as maturation occurred. For PC2, Rif1 (0.0289%) and Paxip1 (0.0288%), which are responsible for the DNA damage response, were gradually transcriptionally repressed with aging (Fig. 2C). Reverse transcription and quantitative real-time PCR (qRT-PCR) assay were carried out to validate the expression patterns of Rif1 and Paxip1 (Supplementary Material, Fig. S7), and consistent results were obtained. RNA in situ hybridization of Rif1 further showed that Rif1 was strongly expressed in the cytoplasm of growing oocytes and their surrounding granulosa cells at the age of 3 weeks, with relatively weak expression in the stroma (Supplementary Material, Fig. S8). Consistently, the expression signal of Rif1 was nearly absent at the age of 9 months. Another contributor to PC2 is Nectin2 (0.0261%), encoding a modulator of T-cell signaling, displayed the opposite expression pattern (Fig. 2C). Interestingly, the well-known ovarian reserve marker Amh (0.0041%), also shown in the PC2 component, was expressed abundantly during PPB and maintained a moderate and steady expression level until MA (Supplementary Material, Fig. S6B). This observation preliminarily indicated that Amh might not be an ideal marker for the early onset of ovarian aging in mice.

Moreover, Gene Ontology (GO) analysis was performed on the contributors to PC1 and PC2 to achieve a global inspection of biological functions. Not surprisingly, many genes exerting effects on PC1 were overrepresented in pathways associated with early organ development (Fig. 3A), such as urogenital system development (1.25%) extracellular matrix organization (1.06%) and epithelial tube morphogenesis (1.18%). Additionally, GO terms such as synapse organization (1.31%), axon development (1.44%), positive regulation of neuron differentiation (1.24%) and neural crest cell development (0.29%; not shown on the plot) revealed the establishment of neural projections during early development, consistent with previous findings that invaded neural crest cells differentiated and formed a dense neural network within the ovarian medulla (29). Notably, among all the GO terms enriched in PC2 genes, DNA repair, the most significant one, possessed the largest contribution score of 0.67% (Fig. 3B). This finding strongly suggested that the gene expression alterations of DNA repair genes might make a considerable contribution to progressive mouse ovarian aging. Other enriched biological processes of PC2 also included the cytoskeleton and its regulation, DNA replication and cell cycle regulation.

Figure 3

Representative GO terms in PC1 and PC2. The top 20 biological processes identified in PC1 (A) and PC2 (B) are shown. The y-axis shows the calculated contributions (%). The sizes of the circles represent the numbers of genes belonging to the terms. The color intensity reflects the negative log10-transformed adjusted P-value, where 0.05 was converted to 1.30.

Clustering identified distinct gene expression patterns of mouse ovaries over the lifespan

To dissect the transcriptional landscape during the ovarian lifespan of mice, after filtering genes with undetectable expression at any time point, we performed a pairwise comparison across all time points to produce global DEGs. Upon the recognition and removal of outliers, 4592 global DEGs (FDR < 0.05; FC of normalized read counts >2) were segregated into 12 distinct expression patterns (Fig. 4A and B; Supplementary Material, Table S2) by unsupervised quality threshold (QT) clustering. This division was further confirmed and visualized by t-distributed stochastic neighbor embedding (t-SNE) dimensionality analysis (Fig. 4C).

Figure 4

Gene expression dynamics and clustering in mouse ovaries during aging. (A) Heatmap showing the expression signatures of 4592 global DEGs during aging in 12 groups identified by QT clustering. The expression values of each gene are represented by raw Z-scores. (B) Line plots show the time trajectories of each cluster. Mean expression values were highlighted in red. (C) t-SNE plot characterizing genes with different expression patterns. Genes with ambiguous patterns were assigned the label C13. The x-axis and y-axis on the t-SNE plot were termed as ‘tSNE1’ and ‘tSNE2’, respectively.

We noted that the top genes contributing to PC1 and PC2 analyzed above could also be classified into specific clusters; for example, Ptk7 belonged to Cluster 6, Tbx2 belonged to Cluster 4, Serpina5 belonged to Cluster 3, Rif1 belonged to Cluster 9, Paxip1 belonged to Cluster 8 and Nectin2 belonged to Cluster 7. Furthermore, GO and KEGG enrichment analysis revealed that Clusters 1–10 had distinct biological functions (Fig. 5; Supplementary Material, Fig. S9), strongly suggesting a structured transcriptional response during ovarian development and aging.

Figure 5

Biological functions of DEG clusters revealed by GO term enrichment analysis. The terms specifically belonging to one DEG cluster are highlighted. P-values were adjusted by Benjamini–Hochberg correction.

Among the 2291 significant GO terms generated from the 10 clusters, five classes of GO terms were specifically attributed to particular clusters, i.e. hypoxia response overrepresented in Cluster 2, epigenetic modification emphasized in Cluster 4, fertilization highlighted in Cluster 8, DNA repair exclusive to Cluster 9 and mitochondrial functions converging on Cluster 10 (Fig. 5). In detail, the genes in Cluster 1 were enriched in immune response-related processes, such as Gpnmb or Zbtb7b, and presented stepwise elevation of expression from birth to senescence (Supplementary Material, Fig. S10A), indicating chronic inflammation in mouse ovaries with age, which is consistent with a previous study (30). Genes responding to hypoxia, such as Hyou1 and Foxo1, were harbored by Cluster 2, with a pattern that peaked sharply at PPB and slightly declined through the remaining stages (Supplementary Material, Fig. S10B). Cluster 2 also overrepresented genes associated with steroid metabolic processes (Fig. 5), such as Amh and Inha (Supplementary Material, Table S2). Genes encoding proteins involved in hormone biosynthesis, including Cyp11a1 and Hsd3b1, were remarkably concentrated in Cluster 3, experiencing transcriptional activation from NN to PB and henceforth fixing their expression at a stable level (Supplementary Material, Fig. S10C), suggesting a close correlation with sexual maturation.

Likewise, the timing of downregulated expression appeared firmly associated with distinct functions, in which Cluster 4 differed from Cluster 5. Genes required for germ cell development, the meiotic cell cycle and epigenetic modification, especially nucleosome organization and DNA methylation involved in gamete generation, were enriched in Cluster 4 (Fig. 5); these genes included Dnmt3a and Mael, which are essential for de novo methylation and piRNA-mediated silencing of transposable elements, respectively (Supplementary Material, Fig. S10D). In contrast, Cluster 5 showed overrepresentation of embryonic organ morphogenesis, including mesenchymal cell proliferation, extracellular matrix organization and urogenital system development, with representative genes including Taf7l and Has2 (Supplementary Material, Fig. S10E). In addition, Cluster 5 also contained several genes such as Zfp36l2 (Supplementary Material, Table S2), encoding RNA-binding proteins involved in maternal mRNA decay process, which is essential for oocyte meiotic maturation (31). Dual expression peaks at NN and MA (Clusters 6 and 7) were observed for genes related to epithelial cell proliferation and regulation of angiogenesis (Sfrp2 and Fgf2; Supplementary Material, Fig. S10F), as well as activation of the immune response (C1qa; Supplementary Material, Fig. S10G) and tumor necrosis factor production (Cyba; Supplementary Material, Fig. S10G). These two clusters differed only in relative expression levels between NN and MA.

Most known oocyte-specific genes accumulated in Clusters 8 and 9, which both displayed a single peak at PPB, with different original expression levels at NN, with overrepresentation of genes with distinct functions. Genes enriched in Cluster 8, such as Zp1 and Zp2, were annotated as single fertilization (Fig. 5 and Supplementary Material, Fig. S10H), while Cluster 9 mainly included genes associated with cell cycle checkpoints and DNA repair, such as Ticrr and Brca2 (Fig. 5 and Supplementary Material, Fig. S10I). In particular, the expression pattern of Cluster 9 was constantly reduced during aging, leading the transcription level at MA to be lower than that at NN. Cluster 10 showed abundant expression at 6 months and encompassed genes regulating mitochondrial organization and oxidative phosphorylation, which have been recently demonstrated to be important pathways in primate ovarian aging (32) (Fig. 5). Typical genes of this cluster included Cox7a2, which encodes a component of the mitochondrial respiratory chain, and Chchd2, which is associated with the mitochondrial intermembrane (Supplementary Material, Fig. S10J).

All these results were further recapitulated by KEGG enrichment analysis (Supplementary Material, Fig. S9). Of note, Cluster 9 was more specifically linked with base excision repair and Fanconi anemia pathway-mediated repair of interstrand DNA cross-links, indicating the different contributions of various DNA repair pathways to ovarian aging and implying that the ovarian capacity for DNA repair might decrease after PPB (Supplementary Material, Figs S9 and S10I).

Causative genes of human POI were significantly associated with Cluster 9

As expected, global DEGs with distinct biological functions were distributed in clusters (Fig. 6A). To examine whether there was any relationship between known genes essential for normal physiological functions and the genes in our identified clusters, we utilized mouse and human genes annotated as being involved in embryonic organ development and aging, and carried out Fisher’s exact test. As shown in Fig. 6B and C, genes responsible for embryonic development in mice were enriched in Cluster 5 (adjusted P-value = 1.29 × 10−5), and their orthologs in humans showed the same enrichment, with even greater significance (adjusted P-value = 5.99 × 10−6). This observation was also consistent with the annotated biological functions of Cluster 5, as described above. Mouse aging-related genes were distributed in all 10 clusters without preference (Fig. 6D). However, human aging-related genes were highlighted in Cluster 6 (adjusted P-value = 0.005), which was associated with tissue remodeling, ossification and regulation of angiogenesis (Figs 5 and 6E).

Figure 6

Distinct distribution patterns of genes with different functions in the clusters. The pie charts and the bar charts show the distribution patterns in the clusters (%) and the estimated odds ratios of (A) global DEGs, (B–E) mouse and human genes annotated with embryonic organ development or aging from the Amigo database, (F) genes related to human DSD and (G) pathogenic genes of human POI generated by a literature search. Odds ratios are presented as the mean ± CI. The dashed lines indicate odds ratio = 1. A two-tailed Fisher’s exact test with FDR adjustment was used for multiple testing. *P < 0.05; **P < 0.01; ***P < 0.001.

Two ovary-related disorders were further investigated. Causative genes of POI (Supplementary Material, Table S3), a typical and common disease with early onset of reduced oocyte quantity or quality, were generated by a literature search. In addition, the list of genes associated with disorders of sexual development (DSD), a pathological state defined as atypical development of chromosomal, gonadal or anatomical sex (33), was mainly derived from a previous report (34). Genes associated with DSD had a slight tendency to accumulate in Cluster 2 (P-value = 0.049, adjusted P-value > 0.05) (Fig. 6F). Interestingly, significant enrichments of causative genes of POI were observed in Cluster 2 (adjusted P-value = 0.015) and Cluster 9 (adjusted P-value = 0.015) (Fig. 6G). The POI-causative genes in Cluster 2 included Amh, Bmpr1b, Inha, etc. (Supplementary Material, Table S3), many of which encoded members or ligands of the transforming growth factor-β superfamily (35). Representative POI-causative genes in Cluster 9 included Brca2, Mcm8 and Msh4 (Supplementary Material, Table S3), which were involved in the DNA DSB repair and mismatch repair pathways, respectively. Loss-of-function variants of these genes have been demonstrated to cause early onset of ovarian dysfunction in humans and mice (36).

The qRT-PCR was further performed to confirm the expression patterns of several known POI-causative genes, and the results were consistent with those obtained from RNA-seq data (Supplementary Material, Fig. S7). For example, Amh and Bmpr1b, which function in hormone signaling, maintained moderate transcription levels in the ovaries of mice during MA and were harbored in Cluster 2 (Fig. 4B and Supplementary Material, Fig. S7). Nobox, which is essential for early folliculogenesis, exhibited a profile corresponding to Cluster 4, which was annotated as germ cell development (Figs 4B and5 and Supplementary Material, Fig. S7). Furthermore, qRT-PCR provided additional evidence that POI-causative genes related to DNA repair in Cluster 9, such as Mcm8, exhibited stepwise reduction in ovarian expression upon PPB (Fig. 4B and Supplementary Material, Fig. S7). All these results strongly suggested the involvement of genes enriched in Clusters 2 and 9 in aging-associated human ovarian disorders.

Discussion

Reproductive aging in mammals has attracted considerable attention for a long time. Although many studies have investigated senescence-associated changes in ovarian transcriptomes to explore the underlying mechanism, most data were produced through isolated comparisons between young and old groups, lacking dynamic expression profiles throughout the whole reproductive lifespan (24,26,32,37). Recently, Cardoso-Moreira et al. performed a bulk RNA sequencing study to provide transcriptomic profiles of seven organs across early developmental time points for humans and six other animal models, highlighting the importance of revealing the developmental trajectories of gene expression during mammalian organ development (38). In this study, by utilizing a mouse model and bulk RNA sequencing, we specifically focused on the transcriptomic landscape of postnatal mouse ovarian development and aging, and identified 12 clusters with distinct gene expression profiles.

To evaluate data quality, we compared our early postnatal (0 day–2 months) data with those from Cardoso-Moreira’s group. Not surprisingly, approximately 6000 out of 18 140 genes (with detected expression in both groups) showed high PCCs (>0.9) (Supplementary Material, Fig. S11A and B), demonstrating apparent consistency between the two studies. Moreover, some genes with similar expression profiles in early development showed divergent trajectories during the subsequent aging process in our results, suggesting their different effects on senescence. For example, Mael and Sfrp2 shared identical expression patterns before adolescence in both studies (PCC = 0.82; P-value = 0.042), but we found that Sfrp2 was upregulated soon after the juvenile stage (Cluster 6), while Mael maintained a unidirectional declining pattern (Cluster 4) (Fig. 4B and Supplementary Material, Fig. S11C). Likewise, the transcription of genes such as Foxo1 (Cluster 2), Hsd3b1 (Cluster 3) and Taf7l (Cluster 5) experienced drastic alterations during early development but remained stable throughout the remaining lifespan, while a large number of genes such as Cyba (Cluster 7), Zp1 (Cluster 8), Ticrr (Cluster 9) and Chchd2 (Cluster 10) exhibited profound expression shifts after adulthood (Supplementary Material, Figs S10 and S11C).

Therefore, our ovarian transcriptomic landscape over the mouse reproductive lifespan provided more details on ovarian gene expression dynamics and highlighted the important transition from the adult stage to aging. These temporal signatures were further annotated with distinct biological functions. For example, Cluster 3 showed peak gene expression at PB and was involved in hormone metabolic processes, indicating that the pattern had a close correlation with sexual maturation. Cluster 6, exhibiting high transcription levels at NN and MA, showed overrepresentation of GO terms including epithelial to mesenchymal transition, which is known to play dual roles in both early-onset organ development and later-onset cancer progression. Moreover, genes regulating fertilization, which are essential for reproductive success, were uniquely enriched in Cluster 8, with a single peak at PPB.

Intriguingly, DNA repair genes critical for oocyte quality were exclusively concentrated in Cluster 9, the expression pattern of which showed a sustained and rapid decline after puberty, leading to lower transcriptional levels at senescence than at birth. As shown in Fig. 3B, the DNA repair pathway was also the most important component of PC2, which segregated different stages during aging.

DNA repair has been demonstrated to be of great importance for maintaining oocyte quality (8), but it remains unknown why and how DNA repair capacity declines with increasing age. Herein, we first revealed the temporal expression profile of DNA repair genes during postnatal ovarian development and aging. Two representative genes are Rif1 and Paxip1 (Fig. 2C). Paxip1 is involved in regulating H3K4 methylation and promoting DNA damage responses (39). Homozygous deficiency in Paxip1 led to lethality at embryonic day 9.5 due to accumulating unrepaired DSBs (40). Specifically, conditional knockout of Paxip1 in adult mouse testes halted the spermatogenesis, indicating the important role of Paxip1 in germ cell development (41). Rif1 is a multifaceted genome maintenance regulator involved in telomere homeostasis, DSB repair pathway choice, replication timing and replication fork protection (42–45). Rif1 has been previously reported to be abundantly expressed in embryonic germ cells, with a low amount of expression in adult ovaries (46). In the present work, we showed that both genes exhibited fast transcriptional activations from NN to PPB and continuous expression decline from PPB to MA, strongly suggesting their roles in maintaining mouse ovarian function. Accelerated decreases in the expression of these genes might indicate a high risk of early onset of ovarian aging. Consistently, loss-of-function mutations in DNA repair genes have frequently been reported in aging-related human ovarian disorders, including POI (36,47).

Due to the high level of genetic heterogeneity (48), it remains a major challenge to identify potential deleterious variants in sporadic POI patients. Multiple studies have demonstrated that mouse models carrying homologous pathogenic variants identified in POI patients manifest human POI phenotypes. However, comprehensive gene expression profiles associated with ovarian aging are unavailable in mice. Herein, we generated a detailed ovarian transcriptomic landscape over the mouse reproductive lifespan and revealed that orthologs of human POI-causative genes were significantly enriched in Clusters 2 and 9. We hypothesized that genes enriched in these two clusters might serve as a good reference for etiology studies of POI and perhaps other human aging-related ovarian disorders.

Interestingly, mutations of genes in Cluster 9, except for known POI-causative genes, also caused female reproductive system defects (Supplementary Material, Fig. S12). For example, mice with a knockout of Dlgap5/Hurp, which acts as a potential cell cycle regulator, exhibited female subfertility without effects on any other functions (49). Not coincidently, Ezhip/AU022751 encodes a cofactor of Polycomb Repressive Complex 2 and is expressed predominantly in the gonads (50). Ezhip-null mice exhibited elevated H3K27me2/3 deposition during late oocyte maturation and accelerated exhaustion of primordial follicle reserve during aging, leading to a considerable reduction in the number of pups each month (50). Fmn2 is a maternal effector gene required for progression through metaphase of meiosis I. Fmn2-deficient mice suffered recurrent silent abortion when mating due to meiotic incompetence in oocytes (51). Kash5/Ccdc155 encodes a member of linker of the nucleoskeleton and cytoskeleton, and is required for homolog pairing during meiotic prophase. Similarly, Kash5−/− mice appeared overtly normal without any obvious hematological pathologies but displayed severe infertile phenotypes with no follicles detected in the ovaries. All these genes emerged as candidates for aging-related human ovarian dysfunctions, and pathogenic variants of these genes deserve further intensive investigation via genetic analysis.

To zdate, state-of-the-art single-cell sequencing has generated new opportunities to dissect complicated cell populations in the ovaries, elucidate interaction networks between cells and model the trajectories of cell lineages during ovarian development (32,52,53). Recently, Fan et al. provided a single-cell atlas of the human adult ovary and revealed the participation of the complement system in follicular remodeling; a representative of this group was C1qa (Cluster 7), which also exhibited gradual transcriptional activation with aging in our results (Supplementary Material, Fig. S10G) (52). Although single-cell sequencing has compelling advantages in establishing high-resolution cell type-specific transcriptional profiles, it is intractable to monitor gene expression trajectories over a long time course because of the tremendous expenditure and enormous computational resources. From this perspective, our results provide a foundation from which to track the molecular mechanisms underlying the alterations of ovarian function during reproductive aging.

To our knowledge, this is the first report of a comprehensive chronological transcriptomic atlas of mouse ovaries throughout the murine reproductive lifespan. By illustrating the profound dynamic changes in differential expression, we identified DNA repair as the most representative pathway for aging, and characterized 10 groups of genes with distinct temporal signatures and coherent biological functions. These findings provide new prospective research directions for mechanistic studies on ovarian aging and etiological studies of aging-associated human ovarian disorders.

Materials and Methods

Hematoxylin and eosin staining of mouse ovaries

Ovaries from C57BL/6 mice (Shanghai Laboratory Animals Center, China) were harvested at the indicated time points, fixed overnight at room temperature in 4% paraformaldehyde, dehydrated, embedded in paraffin and then cut into 5 μm thick sections. Six females were used in each group, and only one ovary was collected from each mouse. One in every five sections was stained with hematoxylin and eosin (H&E) for histological morphology observation and quantification as previously described (54). Representative pictures are shown.

Follicle and corpus luteum counting

Only follicles containing oocytes with a visible nucleus and normal morphology were counted and classified as primordial, primary, secondary or antral follicles. Follicle classification was determined by Pedersen’s system (55). Primordial follicles have small oocytes surrounded by a ring of flattened pregranulosa cells. Primary follicles have growing (enlarged) oocytes surrounded by a ring of cuboidal granulosa cells, while oocytes of secondary follicles are surrounded by two or several layers of cuboidal granulosa cells. Antral follicles have grown oocytes with more than five layers of granulosa cells and a follicle cavity, and they gradually expand to form COCs. Follicles with apoptotic oocytes and scattered granulosa cells were counted as atretic follicles. The corpus luteum is composed of granulosa cells after ovulation. The counting results were calculated as the sum of all H&E-stained sections through the entire ovary.

RNA extraction and sequencing

Ovaries of C57BL/6 mice were harvested at the indicated time points (postnatal day 0, 3 weeks, 2 months, 3 months, 6 months, 9 months and 1 year) and immediately subjected to RNA extraction. An average of 3–8 ovaries from the same time point were pooled in a sample depending on the ovarian size, and 4–5 samples were utilized at each time point. Total RNA was extracted from ovaries using the Allprep DNA/RNA/Protein mini Kit (Qiagen, USA), and RNA samples with RIN values greater than 7 were selected for RNA sequencing. The subsequent procedures, including rRNA depletion, strand-specific library construction and whole-transcriptome sequencing, were performed by Shanghai Biotechnology Corporation (China) on an Illumina HiSeq 2500 platform with a 150 bp paired-end read length. Each sample had approximately 10-G base-calling raw data with Q20 ratios all greater than 90%.

RNA sequencing data processing and availability

The raw sequence data reported in this work are available under accession number CRA003645 that is publicly accessible at https://bigd.big.ac.cn/gsa. Raw data were trimmed by Trimmomatic (version: 0.38) to remove the reads with adapters and low-quality bases. The cleaned reads were aligned to the UCSC Mus musculus reference genome (version: GRCm38/mm10) using Tophat2 (56) (Supplementary Material, Table S1). The quality of alignment was evaluated with RSeQC (57), including sequencing saturation, coverage of the reference genome and distribution of reads on reference genes. Only the uniquely mapped reads were counted by HT-Seq (58). Genes with fragments per kilobase of exon model per million reads mapped (FPKM) > 1 in replicate samples from at least one time point were regarded as expressed genes and retained in the ovary expression matrix for subsequent analysis.

Pairwise DEGs were identified using the R package DESeq2 (FDR < 0.05; log2-transformed FC of normalized read counts >1 or < −1), and global DEGs were defined as the union of all the pairwise DEGs across any two time points, i.e. normalized read counts from each time point were compared to those from every other time point, resulting in n*(n − 1)/2 comparisons altogether, where n was the number of time points.

Identification of genes and pathways with significant contributions to principal components

Read counts yielded by HT-seq were subsequently normalized using DESeq2 with the recommended parameters, as previously reported (59). PCA was performed based on log2-transformed normalized counts, and hierarchical clustering was conducted with default parameters using the R package FactoMineR according to previous work (60). Correlations among different samples were visualized by heatmap (Supplementary Material, Fig. S3).

The correlation coefficients between each gene and the dimension of PCA (P-value <0.05, two-tailed) were calculated by the function ‘dimdesc’ in the R package FactoMineR. Genes with PCCs > 0.85 or < −0.85 were considered contributors to PCs. Representative genes indicated in Fig. 2C were selected from the top 50 contributors. GO enrichment analysis was performed using the R package clusterProfiler (61) in the mouse annotation database org.Mm.eg.db with P-value cut-off < 0.01 and Q-value cut-off < 0.05, adjusted by the Benjamini–Hochberg method. All the significant GO or KEGG terms acquired in this work were generated using the same parameters. The contributions of genes to each PC were obtained by FactoMineR using the formula ‘(var.cos2* 100)/(total cos2 of the component)’. The contributions of GO terms were obtained by simply adding the contributions of all involved genes together. The top 20 GO terms were selected after removing the redundant terms.

Gene expression clustering in the ovaries of mice during aging

QT clustering, an unsupervised partitioning algorithm proposed for gene clustering, was utilized to perform time series analysis (62). The Z-score of each global DEG was calculated based on the log2-transformed average FPKM values per time point and subsequently entered into the function qtclust in the R package flexclust (63) with the maximum radius of clusters = 1.0 and the minimal size of clusters = 25. There were 725 outliers recognized by QT clustering from the overall 5763 global DEGs, and the retained 5038 genes were further clustered into 33 candidate groups. Based on the silhouette coefficients and the PCCs, several raw clusters were discarded due to low quality or were combined together due to similarities (Supplementary Material, Fig. S13). The final patterns were visualized by heatmap based on the Z-score. t-SNE dimensionality analysis was performed with the R package Rtsne with perplexity = 100. GO and KEGG enrichment analyses were conducted using the same parameters for PCA.

Enrichment analysis of mouse and human gene sets in clusters

Mouse and human genes annotated as being involved in embryonic organ development or aging were acquired from GO:0048568 and GO:0007568 (AmiGO 2: http://amigo.geneontology.org/), respectively. The list of POI-causative genes was generated by searching online databases, mainly OMIM (http://www.omim.org) and PubMed (http://www.ncbi.nlm.nih.gov/pubmed). The search keywords used included ‘primary ovarian insufficiency’, ‘premature ovarian failure’, ‘POI’ and ‘POF’. Only genes with multiple lines of evidence (such as clinical reports, mouse models with reproductive phenotypes, ovary-associated functional study, etc.) were classified as causative genes of human POI. Among the genes of target sets in Clusters 1–10, the two-tailed Fisher’s exact test was performed with FDR adjustment for multiple testing.

Statistical analysis

For multiple comparisons of follicle counts and gene expression patterns, one-way ANOVA was performed first to determine whether the observations were significantly changed over the time course (P-value < 0.05). If so, the least significant difference method with Benjamini–Hochberg adjustment in R was utilized to identify the differentiated groups. The error bars of the results are specified in the figure legends.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.