Animal models, oocyte retrieval, and producing offspring production

Four-week-old mice were randomly divided into three groups, and animal models were produced as shown in Fig. 1A. Briefly, each mouse was intraperitoneally injected with 7.5 IU PMSG (Pregnant Mare Serum Gonadotropin, Ningbo Hormone Product Co. Ltd., China), and 7.5 IU hCG (Human Chorionic Gonadotropin, Ningbo Hormone Product Co. Ltd., China) 46–48 h later after PMSG. The interval between two superovulation cycles was seven days. According to the super-ovulated cycles, the groups were, respectively, marked as R1 (one super-ovulated cycle), R3 (three super-ovulated cycles), and R5 (five super-ovulated cycles). At the last injection of hCG, mice were mated with normal male mice immediately to produce offspring, or used to retrieve metaphase II (MII) oocytes 13–14 h later. Ovulated cumulus-oocyte complexes (COCs) were collected from the oviduct ampulla. The surrounding cumulus cells were removed by incubating COCs with 0.1% hyaluronidase. The cleaned oocytes were collected and stored at − 80 °C immediately.

Open in a separate window

Fig. 1

Multi-superovulation has adverse effects on the transcriptome of oocytes. A schedule of repeated superovulation and producing offspring; B and C differential expression genes (DEGs) are identified in R3 and R5; red, up-regulated genes; blue, down-regulated genes; D numbers of DEGs in R3 and R5; E and F GO analysis of DEGs in R3 and R5; (G and H) KEGG analysis of DEGs in R3 and R5

For germinal vesicle oocytes (GV) oocytes, COCs were collected from ovaries at 46–48 h after the last injection of PMSG. The cumulus cells were removed via capillary pipetting.

Transcriptome analysis of oocytes

Smart-seq was used to investigate the transcriptome of the oocytes and was performed by Gene Denovo Biotechnology Co. (Guangzhou, China). Fifty MII oocytes retrieved from≥five female mice were pooled together for each sample, and each group contained three replicates. Total RNA was purified using a TRIzol reagent kit (Invitrogen, USA) according to the handling book. During the purification of total RNA, genomic DNA was removed using DNase. RNA quality was assessed using Agilent 2100 Bioanalyzer (Agilent Technologies, USA) and gel electrophoresis. The purified RNA was reversely transcribed into cDNA using the NEBNext Ultra RNA library Prep Kit for Illumina (NEB, #7530, USA). The cDNA library was sequenced using the Illumina Novaseq 6000 platform. Raw reads were filtered by fastp (version 0.18.0). Clean data were mapped to the reference genome using HISAT2.2.4 with “-rna-standness RF” and other parameters set to defaults. The mapped reads were assembled by StringTie v1.3.1 in a reference-based approach. The expression abundance and variation in transcription regions were quantified by the fragment per kilobase of transcript per million mapped reads (FPKM) value, which was calculated using RSEM software. Differential expression analysis was performed with DESeq2 software [22] via reads count, and differentially expressed genes (DEGs) with false discovery rates (FDR)≤0.05 and absolute fold changes≥two were identified. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using Omicsmart platform (https://www.omicsmart.com/). Trend analysis was performed by ShortTime-series Expression Miner software with the parameters (-pro 20 -ratio 1.0000 [log2(2)=1, log2(1.5)=0.5849625, log2(1.2)=0.2630344]).

Poly (A) tail assay

A poly (A) tail assay (PAT) was performed according to a previous study [19, 21] with minor modifications. Briefly, total RNA was purified from 100 oocytes using the RNAprep Pure Micro Kit (Tiangen, #DP420, Beijing, China). Adaptor P1 (5’-P-GGTCACCTTGATCTGAAGC-NH2–3’) was anchored at the 3’ end with T4 RNA ligase (NEB). Then, it was used for cDNA synthesis using Hifair III 1st Strand cDNA Synthesis Kit (Yeasen, #11139ES60, Shanghai, China) with the P1-antisense primer P2 (5’-GCTTCAGATCAAGGTGACCTTTTT-3’). The product was used as template to amplify specific genes with gene-specific primer and P2 primer at the following conditions: 30 s at 94 °C, 30 s at 58 °C, and 40 s at 72 °C for 35 cycles. Then, the products were subsequently examined using 2.5% agarose gel electrophoresis. The primers used are shown in Table S1.

Whole-genome bisulfite sequencing (WGBS-seq) for small samples

To investigate the whole-genomic DNA methylation in oocytes, WGBS-seq for small samples was used [23]. After hCG injection, the mice were killed by cervical dislocation and the metaphase II (MII) oocytes were collected from the ampulla of fallopian tube. 100 MII oocytes (mouse number: R1:R3:R5=10:20:40) were pooled together for each sample. The samples were subsequently lysed with lysis buffer for 1 h at 37 °C. Bisulfite treatment was performed on the lysed samples using EZ DNA Methylation-Gold kit (Zymo Research) according to the manufacturer’s instructions. The treated DNA was ligated with adapters, and then it was used as a template to construct a sequencing library. The quality of the library was evaluated using a Qubit 3.0 fluorometer, Agilent 2100 Bioanalyzer, and StepOnePlus™ Real-Time PCR System. The library was sequenced using a Hiseq 2500, which was performed by ANOROAD.

To obtain clean data, the raw data were filtered via Trimmomatic(v0.36) software to remove the low quality and adapter sequences. The clean data were mapped to the mouse reference genome (GRCm38) using Bismark(v0.16.3) in paired mode with Bowtie2. The unique mapped reads were used for further analysis. CGIs (CpG islands) and repeat elements were extracted from the UCSC Genome browser. The CGI shore was defined as a 2k region around the CGI, and promoter was defined as the upstream 2k region from transcriptional start site. Differentially methylated regions (DMRs) were identified using DSS software based on β-binomial distribution model and tested by Wald method. DMRs were recognized as fragment length≥200 bp, absolute difference of methylation between groups≥0.25, CpG sites≥5, and absolute Areastat value≥100. KEGG analysis was performed using the online tools (https://www.omicshare.com/tools/Home/Soft/pathwaygsea).

Bisulfite treatment and sequencing

Bisulfite treatment of the oocytes was performed as described in our previous study [24]. Briefly, MII oocytes were retrieved from the ampulla of fallopian tube. Each sample included five oocytes. The samples were incubated with lysis buffer for 37 min, and then treated with 3 M NaOH for 15min. After that, the mixture was embedded in 2% low melting agarose beads. Beads treated with bisulfite for 4 h were washed with TE (Tris–HCl and EDTA), 0.3 M NaOH, and water, respectively. Then, beads with DNA were used as template for nest PCR. PCR products from at least 20 samples were pooled together, and then cloned to a T-vector for sequencing. Sequencing was repeated more than four times for each target sequence. At least ten available clones were used for each target sequence. The methylation level was analyzed using BiQ Analyzer which could filtered the possible repeated sequences.

Body weight, and glucose (GTT) and insulin tolerance test (ITT)

Body weight (g) was tested at the age of birth, three weeks, six weeks, eight weeks, and 12 weeks, respectively. GTT and ITT were performed at eight weeks of age as described in previous studies [24]. Briefly, glucose at a dose of 2 g/kg body weight was given by intraperitoneal injection after 16 h fasting. Blood glucose levels were tested at 0, 30, 60, 90, and 120 min after injection, respectively. ITT was carried out after 4 h fasting. Insulin (Actrapid®, Novo Nordisk) was injected at 10 IU/kg body weight, and blood glucose levels were tested at 0, 30, 60, 90, and 120 min after injection, respectively.

Expression trend analysis

To further understand the effects of the number of superovulation cycles on gene expression, we analyzed the expression trends of genes using ShortTime-series Expression Miner [25]. A total of 455 of 937 genes were filtered (fold change≥2 and correlation between repeats>0). A total of 482 genes were included in the trend analysis. The genes were clustered into 8 profiles (Fig. 2A–H), and the expression of genes in profiles 0, 1 and 3 was significantly reduced with the increase of superovulation times (Fig. 2A, B and D). We further analyzed the enrichment of genes in profiles 0, 1, and 3. GO analysis showed that genes in profile 0 were enriched in metabolic and reproductive processes including prevention of polyspermy, negative regulation of fertilization, and isocitrate metabolism (Fig. S4A). The genes in profile 1 were enriched in protein maturation and glycosylation (Fig. S4B), and the genes in profile 3 were enriched in osteoblast differentiation, and catalytic processes (Fig. S4C). KEGG analysis showed that genes in profile 0 were enriched in protein processing, and acid metabolism (Fig. S4D), genes in profile 1 were enriched in glycan biosynthesis, and ECM-receptor interaction (Fig. S4E), and genes in profile 3 were enriched in the B-cell receptor signaling pathway, and histone metabolism (Fig. S4F). These results suggest that there is a decreasing trend in the expression of genes in oocytes induced by multi-superovulation.

Open in a separate window

Fig. 2

Trend analysis of the groups. The expression trends in R1, R3 and R5 were analyzed using ShortTime-series Expression Miner software. Genes were divided into eight profiles (A–H). With the increase of superovulation time, the expression of genes was significantly reduced in profiles 0, 1, and 3 (A, B and D)

Analysis of the possible reasons for the reduced oocyte maturation rate

Our previous studies report that the increase of superovulation times reduces the oocyte maturation rate and early embryo development [16, 18]. To explore the possible reasons for the reduced oocyte maturation rate, we analyzed DEGs involved in oocyte meiosis. We identified 2 DEGs in R3 including anaphase-promoting complex subunit 6 (Apc6) and mitogen-activated protein kinase 3 (Mapk3), and 6 DEGs in R5, including Apc5/8, Mapk3, RING-box protein 1 (Rbx1), serine/threonine-protein phosphatase 2B catalytic subunit (Cna), and F-box protein 5 (Fbxo5, also known as Emi1). As shown in Fig. 3A, the MAPK signaling pathway is crucial for oocyte maturation because it regulates the activation of APC [26]. Activated APC drives the degradation of cyclin B1 and securing, which activates separase to regulate chromosome separation. EMI1 can inhibit the function of APC to regulate oocyte maturation [27]. Although the increased expression of Rbx1 and decreased expression of Emi1 and Mapk3 may activate APC in R5 oocytes, the expression of APC was reduced (Fig. 3B, D). These findings suggest that the altered expression of these genes might not be the main reason for the reduced oocyte maturation rate in R5.

Open in a separate window

Fig. 3

Analysis of DEGs associated with oocyte maturation. A Protein–protein interaction network of several key genes involved in meiotic resumption; colored nodes, query proteins and first shell of interactors; empty nodes, proteins of unknown 3D structure; filled nodes, a 3D structure is known or predicted; color line, known or predicted interaction and others between proteins; B–G expression of related genes is analyzed using PFKM; same letter, there is no significant difference between groups; different letter, there is significantly statistical difference between groups; H, I poly (A) tail size analysis of several genes in GV and MII oocytes; M, DNA marker; numbers on the left of the panel, size of the band (bp)

Disruption of metabolism in oocytes has adverse effects on oocyte quality [21]. In the present study, KEGG analysis showed that DEGs were enriched in metabolic pathways. We identified 28 DEGs in R3 and 49 DEGs in R5, which were associated with energy metabolism, lipid metabolism, carbohydrate metabolism, and amino acid metabolism (Table S2, and Table S3). Li et al. reported that FAO was enhanced in the meiotic resumption of oocytes [21]. Here, we found that the expression of aldehyde dehydrogenase family 7 member A1 (Aldh7A1), which is related to fatty acid degradation, was decreased in R3 and R5 (Fig. 3E). During oocyte maturation, a reduction in PUFAs is crucial. PUFA is catalyzed by phospholipase A2 (PLA2) [21]. Transcriptomic analysis revealed that the expression of PLA2 receptor 1 (Pla2r1) was increased in R5 oocytes (Fig. 3E). These findings suggest that repeated superovulation may induce PUFAs accumulation in R5 oocytes. High levels of PUFAs may reduce the oocyte maturation rate via decreased the expression of NFKB activating protein (NKAP), which inhibits mRNA degradation by decreasing BTG4 [21]. We analyzed the transcriptomic data and found that the expression of Nkap was respectively reduced by approximately 25% and 20% in R3 and R5, respectively, although the difference was not statistically significant (Fig. 3E). Degradation of transcripts is crucial for meiotic resumption, and BTG4 is an important component that triggers oocyte mRNA decay via the deadenylation of poly(A) tails [19, 20]. However, the expression of Btg4 was similar among the groups (Fig. 3F). Nevertheless, we found that the expression of the mRNA surveillance pathway genes, poly(A) binding protein nuclear 1 (Pabpn1) in R3 oocytes and Pabpn1, RNA binding protein with serine rich domain 1 (Rnps1), regulator of nonsense transcripts 1 (Upf1), cleavage stimulation factor subunit 2 (Cstf2) and RNA-binding protein Musashi (Msi) in R5 oocytes, was altered (Fig. 3E-G). During nuclear polyadenylation, Cstf2 plays a key role in initiating polyadenylation and Pabpn1 contributes to limiting the length of the poly(A) tail at 200–250 nt [28]. Here we found that the expression of Cstf2 was increased and that of Pabpn1 was decreased in R5 oocytes (Fig. 3E). Transcriptomic analysis also showed that the expression of Msi2 [29], which mediates cytoplasmic polyadenylation was decreased and that the expression of Upf1 [30] mediating mRNA decay was increased in R5 oocytes (Fig. 3E, G). These results suggest that the poly(A) tail sizes of mRNAs in oocytes may be affected by repeated superovulation via polyadenylation and deadenylation. To address this question, we examined the polyadenylation length of several transcripts in oocytes (Padi6, Cdk1, Eif6, Pabpc11, and Cops3) using the poly(A) tail (PAT) assay (Fig. 3H, I). Consistent with the findings of a previous study, poly (A) tails were shorter in control MII oocytes than in GV oocytes; but the poly (A) size of Pabpc11 in R3 and R5 GV oocyte was greater than that in R1, and it was shorter for Cdk1 in R3 and R5 GV oocytes compared with R1. The deadenylation of Eif6, Cdk1, and Cops3 in R3 and R5 MII oocytes may be blocked compared with those in R1 (Fig. 3H). These data indicate that disturbed metabolism and mRNA degradation might explain the decreased oocyte maturation rate induced by repeated superovulation.

Multi-superovulation alters whole-genome methylation of oocytes

A previous study reported that superovulation at a lower dose of hormones (<10 IU per mouse) doesn’t affect the methylation patterns of imprinted genes [11]. Nevertheless, the effects of multi-superovulation on global genomic methylation in oocytes are still obscure. Here, we investigated the influence of repeat superovulation on the methylation of oocytes at a base-resolution using a whole-genome bisulfite sequencing (WGBS-seq) approach for small samples. 100 MII oocytes from at least 10 mice were pooled together for each sample. We obtained 700 million clean reads, the Q30 value was 90.89%, and the percentage of clean GCs proportion was 25.39%. Clean reads were mapped to the reference genome mm10 using Bismark (v0.16.3), and the mapped ratio and unique mapped ratio were 34.82% and 32.05%, respectively. The genomic methylation levels of mapped C in R3 and R5 was lower compared with those in R1 (Fig. 4A). According to the context, Cs are classed into three types, CG, CHG, and CHH (H=A, T, or G). Compared with R1, there was a slight increase in the genomic CG methylation level in R3 oocytes, but there was a slight decrease in the CG methylation level in R5 oocytes (Fig. 4B). The methylation levels of CHG and CHH was decreased with the multi-superovulation (Fig. 4B). After that, we further analyzed the methylation level of CG in the genome regions, and the results revealed that there was an increase in the methylation level in the 3’-UTRs of R3 and R5 oocytes (Fig. 4C). The methylation level of CGIs in R3 and R5 oocytes was increased compared with those in R1 (Fig. 4D). The methylation levels of repeat regions in R5 oocytes were also lower compared with R1 oocytes (Fig. 4E). These results indicate that multi-superovulation alters the whole genome DNA methylation of oocytes.

Open in a separate window

Fig. 4

Whole-genome DNA methylation in oocytes. A Whole-genome methylation level of all mapped cytosine; B there are three types of C (cytosine) including CG, CHG, and CHH (H=A, T, or G), and the whole-genome methylation of these types of C are presented; C whole-genomic methylation level of CG in the genomic regions of genes; D the genomic methylation level of CpG islands (CGIs) is analyzed; E the CG methylation levels of repeat regions. *p<0.05, **p<0.01

Differentially methylated region (DMR) analysis

To better evaluate the effect of multi-superovulation on the epigenetic landscape, we analyzed the DMRs in oocytes using DSS software. Compared with R1, we respectively identified 2006 DMRs for R3 including 1085 (54.09%) hyper-DMRs and 921 (45.91%) hypo-DMRs (Fig. 5A), and 2208 DMRs for R5 including 1178 (53.35%) hyper-DMRs and 1030 (46.65%) hypo-DMRs (Fig. 5B). We further examined the distribution of DMRs at genomic elements. For R3 and R5, hyper-DMRs were significantly enriched in simple repeats, exons, CGIs and CGI shores, but depleted from LINEs (long interspersed nuclear elements), LTRs (long terminal repeats), intergenic regions, and non-CGIs (Fig. S5). Notably, hypo-DMRs of R3 were enriched in repeats, exons, CGIs and CGI shores, but depleted from LINEs, LTRs, intergenic regions, and non-CGIs (Fig. 5C). For R5, the hypo-DMRs were significantly depleted from DNA transposons, LINEs, intergenic regions and non-CGIs, but enriched in repeats, exons, introns, CGIs and CGI shores (Fig. 5C). These results suggest that DMRs are not randomly distributed throughout the genome in R3 and R5 oocytes.

Open in a separate window

Fig. 5

Differentially methylated region analysis. A and B Compared with R1, the differentially methylated regions (DMRs) are identified in R3 and R5; hyper, hyper-DMRs; hypo, hypo-DMRs; C the distribution of DMRs in genomic regions is analyzed; * p<0.05, **p<0.01, ***p<0.001; D according to the methylation level, CGIs are divided CGIs into three classes: unmethylated (methylation<20%), intermediate (20%≤methylation≤80%), and methylated (methylation>80%). Hyper-DMRs are mainly distributed at unmethylated CGIs, and hypo-DMRs are mainly distributed at intermediate CGIs; E and F the KEGG enrichment of genes with DMRs at promoters is analyzed

CGIs are common at gene promoters in the vertebrate genome, and are typically in an unmethylated state. Hyper-methylated CGIs would lead to gene silencing by inhibiting the binding of transcription factors with promoters [31]. Therefore, we divided CGIs into three classes: unmethylated (methylation<20%), intermediate (20%≤methylation≤80%), and methylated (methylation>80%). We further analyzed the distribution of DMRs in CGIs, and found that hyper-DMRs in R3 and R5 oocytes were distributed mostly in unmethylated CGIs and hypo-DMRs were distributed mainly in intermediate CGIs (Fig. 5D). These results indicate that unmethylated and intermediate CGIs are prone to be affected by multi-superovulation in oocytes.

To better understand the potential biological role, we analyzed the enrichment of genes with DMRs in promoters in KEGG pathways. The results revealed that genes with DMRs in promoters were enriched in insulin resistance, endocrine resistance, and metabolic pathways (such as fatty acid degradation and glycolysis/gluconeogenesis) in R3 and R5 oocytes (Fig. 5E and F) which suggests that the offspring might have a high risk of metabolic disorders.

To further confirm the methylation of DMRs, four DMRs were analyzed using bisulfite sequencing (BS) as described in a previous study [32]: DMR1, DMR2 and DMR3 were hyper-DMRs and respectively located at Chr11: 72,095,768–72,096,262, Chr18: 10,553,105–10,553,724, and Chr1: 85,832,423–85,833,055. DMR4 is a hypo-DMR and located at Chr5: 124,262,048–124,263,233. Compared with R1 oocytes, DMR1 methylation level was significantly increased in R3, but it was similar between R1 and R5 oocytes (Fig. S6). The methylation of DMR2 was not significantly affected by multi-superovulation in oocytes (Fig. S6). In R3 and R5 oocytes, the methylation level of DMR3 was significantly higher than that in R1 oocytes (Fig. S6). The methylation level of DMR4 in R3 and R5 oocytes was significantly higher than that in R1 oocytes (Fig. S6). These results suggest that multi-superovulation has adverse influence on DNA methylation of oocytes.

We acknowledge the use of Center Experimental Platform of Qingdao Agricultural University and College of Life Science and thank Tao Xu for his invaluable help; we thank all the other people in the labs for their invaluable help.