Pervasive 3′-UTR Isoform Switches During Mouse Oocyte Maturation

Oocyte maturation is the foundation for developing healthy individuals of mammals. Upon germinal vesicle breakdown, oocyte meiosis resumes and the synthesis of new transcripts ceases. To quantitatively profile the transcriptomic dynamics after meiotic resumption throughout the oocyte maturation, we generated transcriptome sequencing data with individual mouse oocytes at three main developmental stages: germinal vesicle (GV), metaphase I (MI), and metaphase II (MII). When clustering the sequenced oocytes, results showed that isoform-level expression analysis outperformed gene-level analysis, indicating isoform expression provided extra information that was useful in distinguishing oocyte stages. Comparing transcriptomes of the oocytes at the GV stage and the MII stage, in addition to identification of differentially expressed genes (DEGs), we detected many differentially expressed transcripts (DETs), some of which came from genes that were not identified as DEGs. When breaking down the isoform-level changes into alternative RNA processing events, we found the main source of isoform composition changes was the alternative usage of polyadenylation sites. With detailed analysis focusing on the alternative usage of 3′-UTR isoforms, we identified, out of 3,810 tested genes, 512 (13.7%) exhibiting significant switches of 3′-UTR isoforms during the process of moues oocyte maturation. Altogether, our data and analyses suggest the importance of examining isoform abundance changes during oocyte maturation, and further investigation of the pervasive 3′-UTR isoform switches in the transition may deepen our understanding on the molecular mechanisms underlying mammalian early development.

Introduction

The proper development of mammalian oocytes is crucial for females to provide high-quality eggs for the generation of new individuals, and improper oocyte maturation usually contributes to female infertility or unhealthy early embryos (Jose-Miller et al., 2007). Oocyte maturation is the last stages of oocyte development, which typically refers to the dynamic process since the first meiotic resumption starting at oocyte germinal-vesicle breakdown (GVBD) (Pan and Li, 2019). Thereafter, oocytes go through metaphase I (MII) and then metaphase II (MII), during which oocyte meiosis is arrested again. The meiosis resumes for the second time at fertilization, and only until this point oocyte maturation is complete (Sen and Caiazza, 2013).

To disentangle the molecular mechanisms during the dynamic process of oocyte maturation, accumulating efforts were attempted from various aspects, including epigenetic regulation (He et al., 2021a), gene expression alteration (Zhao et al., 2020), protein metabolic patterns (Li et al., 2020) as well as their combinations (Gu et al., 2019). It has been shown in literature that the chromatin modification and remodeling in oocytes at the germinal vesicle (GV) stage are largely related to global transcriptional silence before entering meiotic resumption (De La Fuente, 2006). As the pause of transcription continues until zygotic genome activation post fertilization (Bouniol-Baly et al., 1999; Ma et al., 2001), the overall amount of RNA transcripts in oocytes gradually declines, though different genes degrading at various speeds (Zhao et al., 2020). For example, genes related to meiosis I cell cycle process and mitochondrion organization degraded at the top speed during the GV to MII transition in mice (Zhao et al., 2020). With the degradation of such no-longer-needed transcripts, oocytes progressively get ready for fertilization and the sequential early embryonic development (Wu and Dean, 2020).

Recently, researchers have revealed that the RNA degradation machinery is important to alter the transcriptome conformation for developmental processes (Belair et al., 2019; Slobodin et al., 2020; Wu and Dean, 2020), and its regulation is largely dependent on poly(A) tails (Slobodin et al., 2020) and related to 3′-UTR isoforms (Zheng et al., 2018). Although gene-level transcriptomic changes during oocyte maturation have been largely depicted, such changes at the isoform level remain elusive, which would likely impede our advanced understanding in the regulatory mechanisms in RNA degradation and also the functional consequence from fast turnover of selected transcripts. In this study, we focused on the isoform-level transcriptomic differences among main stages of mouse oocyte maturation, demonstrating that the isoform-level changes provided additional information for improving oocyte clustering with respect to their stages. More importantly, we identified 512 genes exhibiting significant changes in 3′-UTR isoform composition between oocytes at the GV and the MII stages, which were linked to essential cellular functions important for the oocyte maturation.

Materials and Methods

Mice

All animal maintenance and experimental procedures were carried out according to guideline of the Animal Experiment Administration Committee of Nanjing Medical University, Jiangsu, China. We used C57BL/6 females to collect oocytes. All mice were given access to food and water, and were maintained on a 12:12 h light-dark artificial lighting cycle, with lights off at 7 p.m.

Collection of Mouse Oocytes

C57BL/6 female mice (6–8 weeks old) were intraperitoneally injected with pregnant mare serum gonadotropin (PMSG, 5 IU, Ningbo Sansheng Pharmaceutical Co., ltd., P.R China), followed by human chorionic gonadotropin (hCG, 5 IU, Ningbo Sansheng Pharmaceutical Co., ltd., P.R China) 48 h later. GV oocytes were isolated from the ovaries 48 h after PMSG injection, and MII oocytes were collected from the ampulla 15 h after hCG injection. For MI oocytes, nude oocytes with clear GVs were cultured in M16 medium in 5% CO2 at 37°C and harvested at 8 h of culture for subsequent analysis as previously described (Jin et al., 2019). All oocytes were collected in hyaluronidase-containing M2 medium (Millipore) drops. Oocytes were selected with morphology, with zona pellucida gently removed by treatment of Acidic Tyrode’s solution (Sigma, cat. no. T1788) and polar body removed by gently pipetting with glass pipette. Thereafter, oocytes were washed with 0.1% BSA/PBS, and transferred into 0.2 ml tubes with the lysis buffer for RNA-seq.

Single Oocyte RNA-Seq Library Preparation and Sequencing

After washed three times with 0.1% BSA/PBS, each oocyte was then transferred into 0.2 ml PCR tubes containing 2 μl lysis buffer. Transcriptome libraries were prepared following the Smart-seq2 protocol. Sequencing libraries were constructed by using KAPA HyperPlus Kit (Kapa Biosystems) according to the manufacturer’s instructions. All libraries were sequenced on the Illumina NovaSeq 6,000 platform with 150-bp paired-end mode at the Annoroad Genomics Company (Beijing, China).

Processing of Single Oocyte RNA-Seq Data

The quality-control metrics of the single oocyte RNA-seq data was checked by FastQC (v0.11.9) (Brown et al., 2017). According to its result, we used Trim Galore (v0.6.4) to trim adapters from reads and remove low quality bases in the data, with arguments length >= 25, stringency = 3 and quality = 25. Then the trimmed reads were then mapped by STAR (v2.7.6a) (Dobin et al., 2013) against the mouse reference genome (mm10) and its corresponding GENCODE gene annotation (VM20). In this process, we also added the information of ERCC spike-ins to the reference. Thereafter, HTSeq (Anders et al., 2015) (v0.12.4) was used to count the reads aligned to each gene for gene-level expression quantification. In parallel, we also used RSEM (Li and Dewey, 2011) (v1.3.1) to align and calculate the expression at both the gene level and the transcript level.

Batch Correction and Cell Clustering

Among a few tested tools, we chose Combat (Zhang et al., 2020) function in the sva (Zhang et al., 2020) (v3.38.0) R package for batch-effect correction according to their performance. In order to prevent overfitting, we used the stage information as the input of the mod parameter. Taking the batch-corrected expression matrix as the input, we used Seurat (Stuart et al., 2019) (v3.2.3) R package for general single oocyte RNA-seq data analysis. In brief, the data were first normalized by the NormalizeData function and the top 2000 highly variable genes (HVG) were selected as the feature genes downstream analysis. Principal component analysis (PCA) was performed using the RunPCA function, and the first five principal components (PCs) were retained for clustering analysis with the FindClusters function in the Seurat package.

Analysis of Differential Expression at Gene and Isoform Levels

Due to the batch-effect issues in our data, we only considered the six GV-stage oocytes and the six MII-stage oocytes for differential expression analysis. The edgeR (Robinson et al., 2010) (v3.26.8) and DESeq2 (Love et al., 2014) (v1.30.1) package was used to identify differentially expressed genes (DEGs) and differentially expressed transcripts (DETs) between the GV and MII oocytes. To achieve this, we took the gene expression matrices calculated by HTSeq and RSEM, and the transcript expression matrix calculated by RSEM as the input of edgeR or DESeq2. Genes or transcripts with false discovery rate (FDR) < 0.05 and expression fold changes >2 were considered as DEGs or DETs. Gene Ontology (GO) enrichment analysis of DEGs and DETs was performed with DAVID (Huang da et al., 2009) (v6.8).

Analysis of Alternative Splicing

SUPPA2 (Trincado et al., 2018) (v2.3) identifies seven types of alternative splicing (AS) events, including skipping exons (SE), alternative 5′ or 3′ splice sites (A5/A3), retained intron (RI), mutually exclusive exons (MX), and alternative first or last exons (AF/AL). GENCODE (VM20) annotation was used as reference to generate all splicing events, and SUPPA2 was used to quantify event inclusion levels with transcript expression (TPM) matrix obtained by RSEM. For each AS event, we define its existence at one stage only when at least half of the oocytes of this stage had valid PSI values (PSI ! = NA). Among the AS events co-existing at both stages, SUPPA2 identified differentially spliced events. rMATS (Shen et al., 2014) (v4.1.0) identifies five types of AS events, including SE, A5, A3, RI, MX. Taking the BAM files of GV and MII oocytes as the input, differentially spliced events were also identified by rMATS.

Detecting Alternative Usage of 3′-UTR Isoforms

Based on the GENCODE gene annotation, DaPars (Xia et al., 2014) (v0.9.1) was adopted to infer unannotated proximal alternative poly(A) sites that result in alternative 3′-UTR isoforms. For its input, we converted bam files to wiggle files by deepTools (Ramirez et al., 2016) (v3.5.0) and bigWigToWig binary downloaded from the UCSC Genome Browser. Then, DaPars calculated the percentage of distal poly(A) site usage index (PDUI) for each oocyte, and tested for significant alternative poly(A) site usage between stages using Fisher’s exact test. The raw p values were subjected to multiple testing adjustment, so that FDR values were calculated. We considered those alternative usage 3′-UTR isoforms between GV and MII oocytes as significant, only if the FDR ≤0.05, the fold change of PDUI ≥1.5, and the absolute differences of PDUI ≥0.3.

Validation of 3′-UTR Isoform Switches

Total RNA was isolated from oocytes with the RNeasy Mini kit (Qiagen). cDNA was synthesized by using SuperScript II Reverse transcriptase (Invitrogen) with oligo-dT primer at 42°C for 90 min. PCR primers were design for the regions shared between alternative 3′-UTR isoforms and the distal 3′-UTR regions using Primer3 (http://Frodo.wi.mit.edu). PCR was performed with Platinum™ SuperFi II DNA Polymerase (Invitrogen, 12361050). PCR conditions: 30 s at 98°C, 15 s at 60°C, and 30 s at 72°C. PCR products were analyzed on a 2% agarose gel. The primer sequences used are listed in Supplementary Table S1.

Results

Isoforms-Level Expression Provides Extra Information to Distinguish Oocyte Stages

We generated two batches of single oocyte RNA-seq data from 16 oocytes at the GV stage, 15 at the MI stage and six at the MII stage. For most of the 37 oocytes, we sequenced each with more than 15 million paired-end reads (Supplementary Figure S1A), and the reads were then mapped to the mouse reference genome (mm10) and quantified against GENCODE gene models (version VM20) at both the gene and isoform levels. Principal component analysis (PCA) revealed a strong batch effect of our data not only at the gene level (Supplementary Figures S1B,C) but also the isoform level (Supplementary Figures S1D,E), which dominated the first principal component (PC 1). To correct the batch effect in our data we tried a few batch-effect correction algorithms, and it turned out that Combat (Johnson et al., 2007) performed the best with respect to oocyte clustering results (see below), which is consistent with previous reports (Chen et al., 2020; Tran et al., 2020).

After batch-effect correction, we took the gene-level (Figures 1A,B) and the isoform-level (Figures 1E,F) expression profiles, respectively, as the input for oocyte clustering analysis. Although there was heterogeneity among transcriptomes of individual oocytes, oocytes at the MII stage were different dramatically from oocytes at the GV stages after entering the meiotic resumption process, featured by many unnecessary transcripts starting to degrade (Wu and Dean, 2020; Zhao et al., 2020). Therefore, we downgraded the unsupervised clustering results which did not entirely separate GV and MII oocytes into distinct clusters. In terms of their transcriptomes, oocytes at the GV and the MI stage were more similar, as shown in a previous study in human oocyte maturation that there were hardly any differentially expressed genes (DEGs) between the GV and MI stages, but a similar number of DEGs between MI vs. MII and between GV vs. MII (Yu et al., 2020). With this knowledge, we did not consider the clear segregation of GV and MI oocytes as a criterion to evaluate the clustering results. According the above criteria, it is clear that the clustering result from gene-level analysis (Figures 1C,D) was less favorable than the result from isoform-level analysis. Hence, we conclude that rather than aggregating isoforms into gene loci, there is extra transcriptomic information contained in individual isoforms, useful for better understanding the dynamics in the transition from meiotic arrest to resumption during mouse oocyte maturation.

FIGURE 1

FIGURE 1. Isoform expression outperformed gene expression in clustering oocytes. (A–C) PCA visualization of individual oocytes based on batch-effect corrected gene-level expression, colored by batches (A), phenotypic stages (B), and expression clustering results (C). (D) Heatmap showing the contingency table between phenotypic stages and clustering results based on batch-effect corrected gene-level expression. (E–G) PCA visualization of individual oocytes based on batch-effect corrected isoform-level expression, colored by batches (E), phenotypic stages (F), and expression clustering results (G). (H) Heatmap showing the contingency table between phenotypic stages and clustering results based on batch-effect corrected isoform-level expression.

Differentially Expressed Genes and Isoforms During Oocyte Maturation

Since we observed an obvious batch effect in our data (Supplementary Figure S1B,D), we applied two strategies including comparison within one batch and application of a covariant model to account for the batch factor, to compare the transcriptomes of GV vs. MI oocytes (Supplementary Table S2) and GV vs. MII oocytes (Supplementary Table S3). We found that different strategies gave similar results, and consistently transcriptomes of GV vs. MII were more divergent than that of GV vs. MI. Hence, we then focused on the comparison of GV vs. MII based on the first batch of our data, which were generated from 6 GV oocytes and six MII oocytes. At the gene level, there were 1,565 genes specifically expressed at the GV stage, compared to 934 genes uniquely expressed at the MII stage (Supplementary Figure S2A). This is consistent with previous results that there were more GV-specific genes than MII-specific ones (Zhao et al., 2020), as transcriptional silence got started during the GV to MII transition. Similar phenomenon was observed at the isoform level (Supplementary Figure S2B). Of note, we considered here only the relative abundance of genes and transcripts, so that we were able to detect genes and transcripts of absolutely lower abundance at the MII stage. Taking this into account, it is very likely that the number of GV-specific genes was understated while the number of MII-specific genes was overstated.

We then analyzed the expression profiles between the two stages for differentially expressed genes (DEGs) and transcripts (DETs). Again, we considered here the relative expression changes albeit a global shift towards downregulation was observed during the transition (Wu and Dean, 2020). With the edgeR algorithm, we detected 591 genes exhibiting relatively higher expression at the MII stage whereas 800 genes showing lower expression (Supplementary Figure S2C). At the isoform level, 220 transcripts were upregulated while 660 were downregulated (Supplementary Figure S2D). Interestingly, when assigning the DETs to their host genes, we found more than one third of these genes (39.1% or 312 genes) not discovered by DEG analysis (Figure 2A), demonstrating again the additional information provided by isoform-level analysis. Of the 312 genes, 17 genes had both transcripts upregulated and downregulated at the MII stage (Figure 2B). In such cases, the expression differences at the gene level could be canceled out, thus becoming undetectable. There were also quite a number of genes (906) only showing differential expression at the gene level but not the isoform level (Figure 2A), possibly due to two different reasons. First, aggregating isoforms that exhibited weak differences made the differences at the gene level much stronger; second, large uncertainty existed in the isoform expression inference based on short sequencing reads, which resulted in larger variability and attenuated statistical power. Thus, when isoform expression quantification became more accurate, the number of genes showing only differences at the gene level should be reduced.

FIGURE 2

FIGURE 2. Differential expression analysis between the GV and MII oocytes. (A) Venn diagram showing the overlap between DEGs and genes with DETs. (B) Out of the 312 non-DEG genes with DETs, Venn diagram showing the distribution of genes with up- and down-regulated transcripts. (C) Gene ontology (GO) term enrichment analysis of gene with up-regulated transcripts. (D) GO term enrichment analysis of gene with down-regulated transcripts. (E) GO term enrichment analysis of the non-DEG genes yet with DETs.

Gene ontology (GO) enrichment analysis showed that genes with upregulated transcripts or upregulated genes at the MII stage were mainly related to neuronal peptides-mediated processes (Figure 2C and Supplementary Figure S2E), which have been shown in low animals to be related to oocyte maturation (Kato et al., 2009; Quiroga Artigas et al., 2020). By contrast, the genes with downregulated transcripts or downregulated genes were enriched in functions related to ribosome biogenesis and translation (Figure 2D and Supplementary Figure S2F), again consistent with previous observations (Zhao et al., 2020). Intriguingly, the 312 genes unidentifiable by DEG analysis but with DETs were functionally related to macromolecule metabolic process and regulation of gene expression (Figure 2E), coinciding with a recent study on porcine oocytes in vitro maturation where the authors found such genes were downregulated during maturation (Brazert et al., 2020). Organelle organization was another enriched GO term for the 312 genes (Figure 2E), and the role of oocyte organelle has been shown to related to oocyte developmental competence and therefore oocyte quality (Reader et al., 2017). Importantly, these functionally critical genes were only identified by isoform-level analysis, further highlighting the usefulness of considering isoform-specific changes during mouse oocyte maturation.

Isoform Changes Predominantly Contributed by Alternative 3′-UTR Usage

We next sought to dissect the alternative RNA processing events that contributed to isoform changes between GV and MII oocytes. According to a benchmarking study (Mehmood et al., 2020) and our experience (He et al., 2021b), we first adopted SUPPA2 (Trincado et al., 2018) and rMATS (Shen et al., 2014) for comparison of annotated alternative events. SUPPA2 investigated seven types of alternative events including skipped exons, retained introns, mutually exclusive exons, alternative 5′ or 3′ splice sites, and alternative first or last exons; however tandem UTRs were not included. In our data, SUPPA2 identified from the whole transcriptome only four alternative events differentially used between GV and MII oocytes at false discovery rate (FDR) of 0.05 (Supplementary Figure S3A). While SUPPA2 might be over-conservative, rMATS did not identify many alternative events of differential usage between the two stages, either (Supplementary Figure S3B), partially because rMATS only focused on internal alternative splicing events, so that entirely ignored structural differences at both ends of a transcript.

Given that the isoform-level analysis indeed provided useful information in distinguishing the two stages from the above analyses, we wondered which alternative events could be the main reasons responsible for the distinct isoform-level differences. By genome-browser enabled visualization investigation of reads mapped to the genome, we found many genes exhibited differences in read coverage near the 3′-ends. According to it, we finally utilized DaPars (Xia et al., 2014) for analysis of alternative polyadenylation sites (PASs) regardless of gene structure annotation. Not surprisingly, DaPars identified 153 events happening in 128 genes at its default significance level (Figure 3A), which we regarded as a set of very stringent criteria (i.e., FDR ≤0.05, fold change of PDUI ≥1.5, and absolute differences of PDUI ≥0.3). When losing the criteria to FDR ≤0.05 and the absolute differences of PDUI ≥0.15, we found 663 (10.7%) significant alternative isoform pairs differentially used in the two stages out of 6,213 tested pairs, or 512 (13.7%) genes harboring such significant events out of 3,810 tested genes (Supplementary Table S4). Hence, we concluded that the main source of isoform-level expression changes during oocyte maturation was the alternative usage of 3′-UTRs, or more specifically, the alternative usage of PASs.

FIGURE 3

FIGURE 3. Alternative 3′-UTR isoform switches during mouse oocyte maturation. (A) Volcano plot showing the differences in percentage of distal poly(A) site usage index (PDUI) between isoform pairs against the associated statistical significance level (FDR). (B) Pie chart demonstrating the asymmetric 3′-UTR isoform switches towards longer UTRs along with oocyte maturation. (C) GO term enrichment analysis of genes exhibiting significantly alternative usage of 3′-UTR isoforms during the oocyte maturation.

Alternative 3′-UTR Isoforms Used Between GV and MII Oocytes

We then asked whether there was a tendency of 3′-UTR length changes among these genes exhibiting alternative usage of 3′-UTR isoform during mouse oocyte maturation. By looking into the direction of PDUI values, we found a clear preference (85.0%) towards maintaining longer 3′-UTR isoforms against degradation during the process (Figure 3B). Even at the loosed PDUI threshold, 68.3% (478 vs. 222) of the significant isoform pairs showed the preference towards longer ones.

Regardless of the length preference, GO enrichment analysis revealed the genes with 3′-UTR isoform switches during oocyte maturation were functionally linked to glutathione metabolic process, cellular response to oxidative stress, maintenance of location in cells, etc. (Figure 3C). It has been reported that the cellular concentration of glutathione changes during oocyte maturation (Luberda, 2005), and glutathione plays certain roles in protecting oocytes against oxidative damage (Brad et al., 2003) and in maintaining the meiotic spindle morphology of the oocytes (Zuelke et al., 1997). Similarly, cellular location maintenance or relocation/redistribution of particular organelles (e.g., mitochondria) in oocytes was closely related to oocyte quality (Mao et al., 2014). Therefore, the genes with 3′-UTR isoform composition changes are of functional importance in safeguarding the process of oocyte maturation.

Of the 128 genes at the top significance, we randomly picked up six genes for independent experimental validation, including three genes that tended to use longer 3′-UTR isoforms during the maturation and the other three with the opposite tendency. In our single-oocyte RNA-seq data, genes Gkap1, Ninj2, and Ndufa8 tended to switch to the longer 3′-UTR isoforms during the GV to MII transition (Figure 4A; Supplementary Figures S4A,B), which was validated by the independent experiments (Figure 4C): Taking the shared region as a reference, the long UTR-specific region either showed weaker bands at the GV stage (Gkap1 and Ninj2) or a stronger band at the MII stage (Ndufa8). In the contrary, genes Gtf2h1, Ndel1, and Rab1b exhibited a tendency to preferably express shorter 3′-UTR isoforms at the MII stage (Figure 4B; Supplementary Figures S4C,D), and were also clearly recapitulated by the validation experiments (Figure 4D).

FIGURE 4

FIGURE 4. Genome-browser visualization and experiment validation of 3′-UTR isoform switches. (A,B) Genome-browser visualization of the genes Gkap1 (A) and Gtf2h1 (B). In each subplot, showing from top to bottom are genomic coordinates, RNA-seq read coverage of 6 GV oocytes (red), RNA-seq read coverage of six MII oocytes (blue), the RefSeq gene annotation, and possible isoforms with alternative polyadenylation sites. The red arrows point out the approximate location of the proximal polyadenylation sites (i.e., the switch points on read coverage between the two stages). (C) Validation of the genes Gkap1, Ninj2, and Ndufa8 with longer 3′-UTR isoforms preferred during the oocyte maturation by independent experiments. Taking the shared region as a reference, the long UTR-specific region either showed weaker bands at the GV stage (Gkap1 and Ninj2) or a stronger band at the MII stage (Ndufa8). (D) Similar to (C), but validation of genes Gtf2h1, Ndel1, and Rab1b, which showed the opposite trend in length of the 3′-UTR isoform usage.

Discussion

In this study, we showcased the importance of isoform-level expression analysis when applying RNA-seq or similar technologies to developmental and reproductive biology, such as in the maturing mouse oocytes. Along with the line that isoform expression quantification provided essential information in more accurately clustering oocytes at different stages, we revealed pervasive alternative 3′-UTR isoform usage during the GV to MII oocyte transition. Although most of the engaged genes tended to have their longer 3′-UTR isoforms more stable, resulting in higher abundance at the MII stage, GO enrichment analysis highlighted the functional importance of these genes regardless of the direction of the changes.

The isoform-level expression changes could serve as fine-tune regulation of the gene-level abundance during the dynamic process. Specifically, since transcription became silenced after entering the meiotic resumption, translational regulation would play a more predominant role in controlling the protein concentration. Reselecting cis-regulatory elements [such as RNA binding protein binding sites (Yamashita and Takeuchi, 2017) and miRNA target sites (Oliveto et al., 2017)] residing in the UTR regions via degrading particular transcripts might be, however, regarded as a global machinery during the transition, to accommodate the transcriptionally silent mode for a certain time period.

In addition, cytoplasmic polyadenylation would be happening after transcriptional silence in maturing oocytes (Reyes and Ross, 2016), which also contributed to translational regulation for protein synthesis (Radford et al., 2008). Therefore, studying the regulation coordination or synergistic regulation of 3′-UTR isoform selection and cytoplasmic polyadenylation of particular transcripts would uncover the detailed mechanisms in shaping the protein landscape of MII oocytes, and also preparing the sufficient molecular materials for upcoming fertilization and early embryonic development, until zygotic genome activation at approximately the two-cell stage for mice and eight-cell stage for human.

From the technical aspect, the present study still adopted the short-read sequencing technology, which could not usually generate sequencing reads longer than 300 nt. This is the exclusive reason for having to infer the isoform-level expression, instead of direct quantification by counting reads. Due to the identifiability of isoform expression inference (Hiller et al., 2009) and large uncertainty introduced in the process (Jiang and Wong, 2009), isoform-level analysis did not show clear advantages over gene-level analysis in expression changes between the two stages. With the rapid development and wide acceptance of the third-generation long-read sequencing technologies (including PacBio and Nanopore), full-length transcriptome sequencing (Wang et al., 2019; Hu et al., 2020) would offer more direct information for isoform-level analysis. Most recently, methods have been developed to sequence poly(A)-tail inclusive full-length transcriptome using the PacBio platform (Legnini et al., 2019; Liu et al., 2019), which would be extremely helpful for studying the detailed molecular mechanisms during the dynamic process of oocyte maturation.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE178836.

Ethics Statement

The animal study was reviewed and approved by the Committee on the Ethics of Animal Experiments of Nanjing Medical University (IACUC1601220, 2101024).

Author Contributions

XW and YH conceived and designed the study. YH and JZ performed all the experiments with help from JY and MX; QC analyzed and interpreted all the data under the supervision of XW; XW and YH drafted the manuscript with input from QC and other authors; XW revised the manuscript. All listed authors have approved the manuscript.

Fundings

This study was funded in part by the National Natural Science Foundation of China (Project No. 32100681, 82101699, and 32170742), the China Postdoctoral Science Foundation (General Program No. 2019M651894), the Nanjing Medical University Research Start-up Fund (Project No. NMUR2020009), and the Opening Fund of Key Laboratory of the Diagnosis and Treatment Research of Reproductive Disorders of Zhejiang Province (No. 2018002).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Prof. Y. Yang (Nanjing Medical University) and Y. Jiang (Nanjing Medical University) for their technical assistance and insightful comments during performing the experiments and preparing the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2021.727614/full#supplementary-material

References

Anders, S., Pyl, P. T., and Huber, W. (2015). HTSeq--a Python Framework to Work with High-Throughput Sequencing Data. Bioinformatics 31 (2), 166–169. doi:10.1093/bioinformatics/btu638

PubMed Abstract | CrossRef Full Text | Google Scholar

Belair, C., Sim, S., Kim, K.-Y., Tanaka, Y., Park, I.-H., and Wolin, S. L. (2019). The RNA Exosome Nuclease Complex Regulates Human Embryonic Stem Cell Differentiation. J. Cel Biol. 218 (8), 2564–2582. doi:10.1083/jcb.201811148

CrossRef Full Text | Google Scholar

Bouniol-Baly, C., Hamraoui, L., Guibert, J., Beaujean, N., Szöllösi, M. S., and Debey, P. (1999). Differential Transcriptional Activity Associated with Chromatin Configuration in Fully Grown Mouse Germinal Vesicle Oocytes1. Biol. Reprod. 60 (3), 580–587. doi:10.1095/biolreprod60.3.580

PubMed Abstract | CrossRef Full Text | Google Scholar

Brad, A. M., Bormann, C. L., Swain, J. E., Durkin, R. E., Johnson, A. E., Clifford, A. L., et al. (2003). Glutathione and Adenosine Triphosphate Content of In Vivo and In Vitro Matured Porcine Oocytes. Mol. Reprod. Dev. 64 (4), 492–498. doi:10.1002/mrd.10254

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, J., Pirrung, M., and McCue, L. A. (2017). FQC Dashboard: Integrates FastQC Results into a Web-Based, Interactive, and Extensible FASTQ Quality Control Tool. Bioinformatics 33 (19), 3137–3139. doi:10.1093/bioinformatics/btx373

PubMed Abstract | CrossRef Full Text | Google Scholar

Brązert, M., Kranc, W., Nawrocki, M., Sujka-Kordowska, P., Konwerska, A., Jankowski, M., et al. (2020). New Markers for Regulation of Transcription and Macromolecule Metabolic Process in Porcine Oocytes during In Vitro Maturation. Mol. Med. Rep. 21 (3), 1537–1551. doi:10.3892/mmr.2020.10963

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, W., Zhang, S., Williams, J., Ju, B., Shaner, B., Easton, J., et al. (2020). A Comparison of Methods Accounting for Batch Effects in Differential Expression Analysis of UMI Count Based Single Cell RNA Sequencing. Comput. Struct. Biotechnol. J. 18, 861–873. doi:10.1016/j.csbj.2020.03.026

PubMed Abstract | CrossRef Full Text | Google Scholar

De La Fuente, R. (2006). Chromatin Modifications in the Germinal Vesicle (GV) of Mammalian Oocytes. Develop. Biol. 292 (1), 1–12. doi:10.1016/j.ydbio.2006.01.008

CrossRef Full Text | Google Scholar

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). STAR: Ultrafast Universal RNA-Seq Aligner. Bioinformatics 29 (1), 15–21. doi:10.1093/bioinformatics/bts635

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, C., Liu, S., Wu, Q., Zhang, L., and Guo, F. (2019). Integrative Single-Cell Analysis of Transcriptome, DNA Methylome and Chromatin Accessibility in Mouse Oocytes. Cell Res. 29 (2), 110–123. doi:10.1038/s41422-018-0125-4

PubMed Abstract | CrossRef Full Text | Google Scholar

He, M., Zhang, T., Yang, Y., and Wang, C. (2021a). Mechanisms of Oocyte Maturation and Related Epigenetic Regulation. Front. Cel Dev. Biol. 9, 654028. doi:10.3389/fcell.2021.654028

CrossRef Full Text | Google Scholar

He, Y., Chen, Q., Dai, J., Cui, Y., Zhang, C., Wen, X., et al. (2021b). Single-cell RNA-Seq Reveals a Highly Coordinated Transcriptional Program in Mouse Germ Cells during Primordial Follicle Formation. Aging Cell 20, e13424. doi:10.1111/acel.13424

PubMed Abstract | CrossRef Full Text | Google Scholar

Hiller, D., Jiang, H., Xu, W., and Wong, W. H. (2009). Identifiability of Isoform Deconvolution from junction Arrays and RNA-Seq. Bioinformatics 25 (23), 3056–3059. doi:10.1093/bioinformatics/btp544

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, Y., Shu, X.-S., Yu, J., Sun, M.-a., Chen, Z., Liu, X., et al. (2020). Improving the Diversity of Captured Full-Length Isoforms Using a Normalized Single-Molecule RNA-Sequencing Method. Commun. Biol. 3 (1), 403. doi:10.1038/s42003-020-01125-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, D. W., Sherman, B. T., and Lempicki, R. A. (2009). Systematic and Integrative Analysis of Large Gene Lists Using DAVID Bioinformatics Resources. Nat. Protoc. 4 (1), 44–57. doi:10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, H., and Wong, W. H. (2009). Statistical Inferences for Isoform Expression in RNA-Seq. Bioinformatics 25 (8), 1026–1032. doi:10.1093/bioinformatics/btp113

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, Y., Yang, M., Gao, C., Yue, W., Liang, X., Xie, B., et al. (2019). Fbxo30 Regulates Chromosome Segregation of Oocyte Meiosis. Cell. Mol. Life Sci. 76 (11), 2217–2229. doi:10.1007/s00018-019-03038-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, W. E., Li, C., and Rabinovic, A. (2007). Adjusting Batch Effects in Microarray Expression Data Using Empirical Bayes Methods. Biostatistics 8 (1), 118–127. doi:10.1093/biostatistics/kxj037

PubMed Abstract | CrossRef Full Text | Google Scholar

Jose-Miller, A. B., Boyden, J. W., and Frey, K. A. (2007). Infertility. Am. Fam. Physician 75 (6), 849–856.

PubMed Abstract | Google Scholar

Kato, S., Tsurumaru, S., Taga, M., Yamane, T., Shibata, Y., Ohno, K., et al. (2009). Neuronal Peptides Induce Oocyte Maturation and Gamete Spawning of Sea Cucumber, Apostichopus Japonicus. Develop. Biol. 326 (1), 169–176. doi:10.1016/j.ydbio.2008.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Legnini, I., Alles, J., Karaiskos, N., Ayoub, S., and Rajewsky, N. (2019). FLAM-seq: Full-Length mRNA Sequencing Reveals Principles of Poly(A) Tail Length Control. Nat. Methods 16 (9), 879–886. doi:10.1038/s41592-019-0503-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., and Dewey, C. N. (2011). RSEM: Accurate Transcript Quantification from RNA-Seq Data with or without a Reference Genome. BMC Bioinformatics 12, 323. doi:10.1186/1471-2105-12-323

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Zhu, S., Shu, W., Guo, Y., Guan, Y., Zeng, J., et al. (2020). Characterization of Metabolic Patterns in Mouse Oocytes during Meiotic Maturation. Mol. Cel 80 (3), 525–540. doi:10.1016/j.molcel.2020.09.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Nie, H., Liu, H., and Lu, F. (2019). Poly(A) Inclusive RNA Isoform Sequencing (PAIso−seq) Reveals Wide-Spread Non-adenosine Residues within RNA Poly(A) Tails. Nat. Commun. 10 (1), 5292. doi:10.1038/s41467-019-13228-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2. Genome Biol. 15 (12), 550. doi:10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Luberda, Z. (2005). The Role of Glutathione in Mammalian Gametes. Reprod. Biol. 5 (1), 5–17.

PubMed Abstract | Google Scholar

Ma, J., Svoboda, P., Schultz, R. M., and Stein, P. (2001). Regulation of Zygotic Gene Activation in the Preimplantation Mouse Embryo: Global Activation and Repression of Gene Expression1. Biol. Reprod. 64 (6), 1713–1721. doi:10.1095/biolreprod64.6.1713

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, L., Lou, H., Lou, Y., Wang, N., and Jin, F. (2014). Behaviour of Cytoplasmic Organelles and Cytoskeleton during Oocyte Maturation. Reprod. BioMed. 28 (3), 284–299. doi:10.1016/j.rbmo.2013.10.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehmood, A., Laiho, A., Venäläinen, M. S., McGlinchey, A. J., Wang, N., and Elo, L. L. (2020). Systematic Evaluation of Differential Splicing Tools for RNA-Seq Studies. Brief Bioinform. 21 (6), 2052–2065. doi:10.1093/bib/bbz126

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliveto, S., Mancino, M., Manfrini, N., and Biffo, S. (2017). Role of microRNAs in Translation Regulation and Cancer. WJBC 8 (1), 45–56. doi:10.4331/wjbc.v8.i1.45

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, B., and Li, J. (2019). The Art of Oocyte Meiotic Arrest Regulation. Reprod. Biol. Endocrinol. 17 (1), 8. doi:10.1186/s12958-018-0445-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Quiroga Artigas, G., Lapébie, P., Leclère, L., Bauknecht, P., Uveira, J., Chevalier, S., et al. (2020). A G Protein-Coupled Receptor Mediates Neuropeptide-Induced Oocyte Maturation in the Jellyfish Clytia. Plos Biol. 18 (3), e3000614. doi:10.1371/journal.pbio.3000614

PubMed Abstract | CrossRef Full Text | Google Scholar

Radford, H. E., Meijer, H. A., and de Moor, C. H. (2008). Translational Control by Cytoplasmic Polyadenylation in Xenopus Oocytes. Biochim. Biophys. Acta Gene Regul. Mech. 1779 (4), 217–229. doi:10.1016/j.bbagrm.2008.02.002

CrossRef Full Text | Google Scholar

Ramírez, F., Ryan, D. P., Grüning, B., Bhardwaj, V., Kilpert, F., Richter, A. S., et al. (2016). deepTools2: a Next Generation Web Server for Deep-Sequencing Data Analysis. Nucleic Acids Res. 44 (W1), W160–W165. doi:10.1093/nar/gkw257

PubMed Abstract | CrossRef Full Text | Google Scholar

Reader, K., Stanton, J.-A., and Juengel, J. (2017). The Role of Oocyte Organelles in Determining Developmental Competence. Biology 6 (3), 35. doi:10.3390/biology6030035

PubMed Abstract | CrossRef Full Text | Google Scholar

Reyes, J. M., and Ross, P. J. (2016). Cytoplasmic Polyadenylation in Mammalian Oocyte Maturation. WIREs RNA 7 (1), 71–89. doi:10.1002/wrna.1316

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 26 (1), 139–140. doi:10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Sen, A., and Caiazza, F. (2013). Oocyte Maturation A story of Arrest and Release. Front. Biosci. S5, 451–477. doi:10.2741/s383

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, S., Park, J. W., Lu, Z.-x., Lin, L., Henry, M. D., Wu, Y. N., et al. (2014). rMATS: Robust and Flexible Detection of Differential Alternative Splicing from Replicate RNA-Seq Data. Proc. Natl. Acad. Sci. USA 111 (51), E5593–E5601. doi:10.1073/pnas.1419161111

PubMed Abstract | CrossRef Full Text | Google Scholar

Slobodin, B., Bahat, A., Sehrawat, U., Becker-Herman, S., Zuckerman, B., Weiss, A. N., et al. (2020). Transcription Dynamics Regulate Poly(A) Tails and Expression of the RNA Degradation Machinery to Balance mRNA Levels. Mol. Cel 78 (3), 434–444. doi:10.1016/j.molcel.2020.03.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Stuart, T., Butler, A., Hoffman, P., Hafemeister, C., Papalexi, E., Mauck, W. M., et al. (2019). Comprehensive Integration of Single-Cell Data. Cell 177 (7), 1888–1902. doi:10.1016/j.cell.2019.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Tran, H. T. N., Ang, K. S., Chevrier, M., Zhang, X., Lee, N. Y. S., Goh, M., et al. (2020). A Benchmark of Batch-Effect Correction Methods for Single-Cell RNA Sequencing Data. Genome Biol. 21 (1), 12. doi:10.1186/s13059-019-1850-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Trincado, J. L., Entizne, J. C., Hysenaj, G., Singh, B., Skalic, M., Elliott, D. J., et al. (2018). SUPPA2: Fast, Accurate, and Uncertainty-Aware Differential Splicing Analysis across Multiple Conditions. Genome Biol. 19 (1), 40. doi:10.1186/s13059-018-1417-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., You, X., Langer, J. D., Hou, J., Rupprecht, F., Vlatkovic, I., et al. (2019). Full-length Transcriptome Reconstruction Reveals a Large Diversity of RNA and Protein Isoforms in Rat hippocampus. Nat. Commun. 10 (1), 5009. doi:10.1038/s41467-019-13037-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, D., and Dean, J. (2020). EXOSC10 Sculpts the Transcriptome during the Growth-To-Maturation Transition in Mouse Oocytes. Nucleic Acids Res. 48 (10), 5349–5365. doi:10.1093/nar/gkaa249

PubMed Abstract | CrossRef Full Text | Google Scholar

Xia, Z., Donehower, L. A., Cooper, T. A., Neilson, J. R., Wheeler, D. A., Wagner, E. J., et al. (2014). Dynamic Analyses of Alternative Polyadenylation from RNA-Seq Reveal a 3′-UTR Landscape across Seven Tumour Types. Nat. Commun. 5, 5274. doi:10.1038/ncomms6274

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamashita, A., and Takeuchi, O. (2017). Translational Control of mRNAs by 3'-Untranslated Region Binding Proteins. BMB Rep. 50 (4), 194–200. doi:10.5483/bmbrep.2017.50.4.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, B., Doni Jayavelu, N., Battle, S. L., Mar, J. C., Schimmel, T., Cohen, J., et al. (2020). Single-cell Analysis of Transcriptome and DNA Methylome in Human Oocyte Maturation. PLoS One 15 (11), e0241698. doi:10.1371/journal.pone.0241698

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Parmigiani, G., and Johnson, W. E. (2020). ComBat-seq: Batch Effect Adjustment for RNA-Seq Count Data. NAR Genom Bioinform 2 (3), lqaa078. doi:10.1093/nargab/lqaa078

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Z.-H., Meng, T.-G., Li, A., Schatten, H., Wang, Z.-B., and Sun, Q.-Y. (2020). RNA-seq Transcriptome Reveals Different Molecular Responses during Human and Mouse Oocyte Maturation and Fertilization. BMC Genomics 21 (1), 475. doi:10.1186/s12864-020-06885-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, D., Wang, R., Ding, Q., Wang, T., Xie, B., Wei, L., et al. (2018). Cellular Stress Alters 3′UTR Landscape through Alternative Polyadenylation and Isoform-specific Degradation. Nat. Commun. 9 (1), 2268. doi:10.1038/s41467-018-04730-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zuelke, K. A., Jones, D. P., and Perreault, S. D. (1997). Glutathione Oxidation Is Associated with Altered Microtubule Function and Disrupted Fertilization in Mature Hamster Oocytes1. Biol. Reprod. 57 (6), 1413–1419. doi:10.1095/biolreprod57.6.1413

PubMed Abstract | CrossRef Full Text | Google Scholar