Aberrant spliceosome expression and altered alternative splicing events correlate with maturation deficiency in human oocytes

ABSTRACT

Different strategies of ovarian stimulation are widely used in IVF to retrieve mature metaphase II (MII) oocytes for fertilization. On average, approximately 70% of recovered oocytes are mature, while personalized administration of hCG and/or GnRH agonist trigger and in vitro maturation (IVM) management can further improve the maturation rate. However, even under such conditions, a complete absence of oocyte maturation is still observed sporadically. The probable causes for such maturation-deficient (MD) oocytes – which arrest abnormally at metaphase I (MI) stage – are still under investigation. In the present study, using single-cell transcriptomic RNA sequencing (RNA-seq) and differential expression analysis, we showed that gene expression profiles were aberrant, and alternative splicing (AS) patterns were changed in MD oocytes when compared with normally mature (MN) oocytes. Gene ontology (GO) enrichment demonstrated that the differently expressed genes (DEGs) were mostly correlated with pre-mRNA splicing, RNA transportation, RNA processing, and mRNA regulation. Subsequently, analysis of AS events revealed that genes with altered AS patterns were primarily associated with metabolism and cell cycle. With these findings, we have demonstrated aberrant gene expression in complete maturation-deficient oocytes, and we propose that alterations in post-transcriptional regulation constitute a potential underlying mechanism governing oocyte maturation.

KEYWORDS:

• Oocyte maturation deficiency
• transcriptome
• spliceosome

Introduction

Intrauterine insemination (IUI), in vitro fertilization (IVF), intracytoplasmic sperm injection (ICSI), and other assisted reproductive techniques have helped many infertile couples to have their own babies, and now comprise 1–3% of all births [1]. Different strategies of ovarian stimulation have been used according to a patient’s specific situation – including the long protocol, short protocol, and antagonist protocol. Although with individualized protocols and combined with in vitro maturation (IVM) techniques, we typically achieve 60–80% of oocytes maturing; in some rare cases, none of the recovered oocytes matures. These maturation-deficient (MD) oocytes may arrest abnormally at the germinal vesicle (GV) stage or metaphase of meiosis I (MI) stage, making it impossible to achieve insemination [2]. Unfortunately, only a few genes and their related mechanisms of action have been identified as being responsible for human oocyte maturation arrest; most causes of complete oocyte maturation defects are currently unexplained and require further exploration.

Oocyte maturation is comprehensively regulated at transcriptional, translational, endocrinologic, and metabolic levels [3,4]. Oogenesis is initiated with primordial germ cells (PGCs) at the embryonic stage, and oocytes arrest at the diplotene stage of meiotic prophase I for several years. Meiotic progression then resumes with puberty and at each subsequent menstrual period [5]. Hormonal signaling regulates menstrual cyclicity, while oocyte-gonadal somatic cell interactions control the timing of oocyte meiotic maturation [6]. After first meiosis, oocytes are ovulated and enter a second arrest at metaphase II (MII) [7]. Since transcriptional events are almost completed before full growth of GV oocytes occurs, and are relatively silenced thereafter [8], meiotic resumption and maturation – including meiotic spindle assembly and cytoplasmic maturation – are principally based upon already-existing maternal transcripts [9]. Therefore, integrated transcriptome maintainance is of great importance in oocyte development. Using single-cell RNA sequencing (RNA-seq), we and other groups have previously contributed to an expression atlas of human oocytes and embryos [10–13]. In addition, comparisons between in vivo and in vitro matured oocytes revealed alterations in energy metabolism [14], and comparisons between young women and those at advanced maternal ages showed divergent expression patterns concerning mitochondria, cell cycle kinetics, and cytoskeleton [15]. Recently, a study of failed-to-mature oocytes from rhesus monkeys indicated deficient transcriptional regulation [16]; however, the transcriptomic landscape of meiotically incompetent human oocytes remains largely obscure.

Accordingly, pre-mRNA alternative splicing (AS) was implicated in providing cells with the capacity to respond to stimuli, and AS also contributes to transcriptome complexity [17] – with any dysregulation resulting in aberrant phenotypes and diseases [18]. Splicing is executed by the spliceosome, which is composed of 5 major U-rich small nuclear ribonucleoproteins (U1/U2/U4/U5/U6 snRNPs), NineTeen complex (NTC), NineTeen complex-related (NTR), and many other factors exemplified by serine/arginine-rich (SR) proteins and hetero-nuclear ribonucleoproteins (hnRNPs). In oocyte maturation as well as in spermatogenesis, spliceosomal proteins and AS play crucial roles in regulating stage-specific transcripts [19]. For example, nuclear m6A (N6-methyladenosine) reader YTHDC1 regulates AS during oocyte development [20], and the RNA-binding protein SRSF3 regulates AS to maintain transcriptome integrity [21]. However, most AS regulators and their functions in oocyte maturation remain to be elucidated.

In the current study, by single-cell RNA sequencing analysis, we performed a detailed investigation into the transcriptional profiles and AS patterns of human MD oocytes, and compared them with normally mature (MN) oocytes. High-quality, single-cell RNA-seq analysis provided us with the opportunity to uncover aberrations in expressed genes and pathways, and we then specifically analyzed the altered AS events and splice variants. The present study thus clarifies the post-transcriptional regulation of genes as a mechanism subserving oocyte maturation.

Materials and methods

Ethical approval, participants, and human oocyte collection

The present study was reviewed and approved by the Institutional Review Board of Peking University Third Hospital (No. 2018013). All of the oocytes were voluntarily donated by IVF patients who had been clearly informed of our experimental aims and signed the appropriate informed consent forms.

MD oocytes (Y1, Y2, and Y3) were obtained from individual patients whose oocytes failed to mature completely. These MD oocytes used in this study were collected in the second or third IVF cycle as indicated. All the oocytes were retrieved at GV stage, remaining at MI after germinal vesicle breakdown (GVBD) until the third day of in vitro culture (Supplementary Figure 1). Such immature oocytes are usually discarded, but herein we collected them for transcriptomic analysis. MN oocytes (B1, B2, and B3) were voluntarily donated from individual IVF patients who underwent normal retrieval of MII oocytes that exhibited normal maturation rates. The MII stage oocytes were used for patient IVF treatments, whereas GV stage oocytes were cultured to MII and donated to our study and referred to as MN. No differences were found among patients regarding their demographics and baseline characteristics (Table 1).

Table 1. Primary physiologic indices of patient.

Patient	Maternal age	BMI	IVF cycle numbers	Stimulation protocol	Basal FSH (mIU/ml)	Basal E2(pmol/L)	Basal P(nmol/L)	Basal LH (mIU/ml)	Number of retrieved oocytes	number of mature oocytes
Y1	29	22.27	3	antagonist	1.90	94.30	0.64	1.46	19	0
Y2	33	20.7	2	antagonist	6.60	96.50	1.10	5.81	12	0
Y3	27	20.96	2	antagonist	8.00	250.00	1.41	12.10	10	0
B1	31	19.59	1	antagonist	1.90	73.40	/	0.81	12	8
B2	29	23.04	1	antagonist	2.87	366.00	2.90	6.19	15	12
B3	33	20.80	1	antagonist	9.10	490.00	2.00	5.00	13	10
/, not recorded.

Clinical procedures and oocyte treatment were identical for both MD and MN patients. Ovarian stimulation was performed under an antagonist protocol after hCG (human chorionic gonadotropin) administration, which was then followed by ultrasound-guided needle aspiration of follicles 46 to 48 hours later. Both MN and MD oocytes were retrieved at GV stage (Day 0) and the cumulus cells were denuded 2 hours later. Mature MII oocytes from MN patients were used for their IVF treatments, and their remaining GV oocytes were cultured for another 24–48 hours (Day 1–2) to promote maturation to MII (MN oocytes). For MD patients, all their retrieved oocytes remained unmatured at the MI stage even on Day 3. Both MN and MD oocytes were ultimately collected and prepared for molecular analysis .

cDNA synthesis, library preparation, and sequencing

Cumulus cells were removed and individual oocytes were exposed to hyaluronidase to remove all remaining cumulus cells, and then washed several times with Ca²⁺/Mg²⁺-free PBS. Total RNA was isolated using the PicoPure^TM RNA Extraction kit (Life Technologies, Carlsbad, CA, USA), including RNase-free DNase treatment. We then used oligo (dT) primers and template-switching oligo primers to synthesize first-strand cDNA. In situ PCR primer-guided PCR amplification was performed with a SMARTer Ultra Low Input RNA HV kit (Clontech, USA), following the manufacturer’s protocol. The quality and quantity of the resulting cDNA (after purification by AMPure XP Beads [A63882; Beckman Coulter Inc., USA]) was evaluated with a Qubit 3.0 fluorometer (Thermo Fisher Scientific, USA).

cDNA library preparation began with cDNA fractionation by ultrasonic waves using a Bioruptor® sonication system (Diagenode Inc., Belgium) followed by terminal repair, adaptor addition, and purification using AMPure XP Beads. Then, 200–300 bp DNA fragments were recovered with a CWBIO Gel Extraction Kit (CW2302A; Cwbiotech Inc., China) after 2% agarose gel electrophoresis. The insert size of the constructed cDNA library was evaluated with an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) and quantified by quantitative real-time PCR using Bio-Rad CFX 96. Finally, Illumina paired-end 150 bp sequencing was performed on a HiSeq 2500 (Illumina, USA). All sequencing data are available at the Gene Expression Omnibus (GSE129123).

Single-cell transcriptome analysis: expression profiling and pathway enrichment

After adaptor-moving and quality-control filtration of raw reads with NGSQC [22], the remaining clean reads for each sample were mapped to the annotated human genome (hg19.ucsc.fa) with HISAT under default settings [23]. We estimated mapping distribution, homogeneity, and saturation with online Picard (http://broadinstitute.github.io/picard/). Afterward, fragment counts normalized per kilobase of exon per million mapped fragments (FPKM) were estimated using Cuffquant and Cuffnorm in Cufflinks (v2.0.2) [24]. Unsupervised hierarchical clustering was conducted using Cluster 3.0 and visualized using Java TreeView. Expression levels were estimated for all isoforms of a gene, and we determined differentially expressed genes (DEGs) using Cuffdiff in Cufflink, with |log2-fold change| ≥ 1 and a Q-value <0.05. The resulting DEGs were visualized by volcano plot.

For the annotation and enrichment of DEGs, Gene Ontology (GO) categories were classified by 2-sided Fisher’s exact probability test and Chi-squared test using DAVID (https://david.ncifcrf.gov/). Pathway enrichment analysis was performed using a similar method with reference to different databases, including the Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome, and BioCyc. The false discovery rate (FDR) was calculated to adjust the p-value, and GO terms or KEGG pathways were considered to be statistically significant with a FDR value <0.05.

Alternative splicing (AS) analysis and enrichment of differential AS genes

Alternative splicing events and global expression patterns of different exon-exon junction reads were analyzed using ASprofile [25] for each sample. Briefly, all pairs of transcripts within a gene were compared to determine exon-intron structural differences that indicated an AS event. The FPKM of these events were calculated and splicing ratios across samples were compared, indicating the dynamics of AS profiles. Accordingly, 12 major types of AS events were implemented.

The classification of AS and evaluation of differentially alternatively spliced genes were analyzed using rMATS [26] between MD and MN oocytes. We estimated all of the following 5 major types of AS: skipped exon, SE; alternative 5ʹ splice site, A5SS; alternative 3ʹ splice site, A3SS; mutually exclusive exons, MXE; and retained intron, RI. The inclusion level (IncLevel) of splicing sites was quantified by the mapped transcript ratio, and regarded as differential AS (DAS) events according to the likelihood-ratio test with a FDR value <0.05.

Pathway enrichment was performed for genes with at least 1 DAS, and the enrichment was similar for the DEGs described previously.

Results

Transcriptional profiles of single human MD and MN oocytes

Although we found no differences in patients’ physiologic indices and although patients were treated with similar IVF procedures, MD oocytes remained arrested at the MI stage while sharing similar sizes and morphologies with MN MII oocytes (Table 1, Supplementary Fig. S1A). From the 3 MD and 3 MN oocytes analyzed in our study, we obtained 6 individual, high-quality RNA-seq datasets after total RNA extraction, subsequent cDNA library construction, and sequencing. Approximately 67 million reads were detected for each sample, with a clean rate exceeding 99.7% (Supplementary Fig. S1B), and nearly 96% were mapped to the human genome on average (Supplementary Fig. S1C) – with over 60% constituting exonic reads (Supplementary Fig. S2A), thus displaying a high coverage. Alignment statistics also showed high homogeneity and saturation of sequencing (Supplementary Fig. S2B, S2C).

In total, 29,911 transcripts were detected as expressed, in which 62% were annotated as protein-coding genes, followed by pseudogenes (15%) and long non-coding RNA (lncRNAs) genes (11%) (Figure 1(a)). Among these transcripts, 4531 were only detected in MD oocytes, 4532 were only expressed in MN oocytes, and 20,848 were jointly detected in samples from both groups (Figure 1(b)). We displayed all 29,911 genes among 6 samples in a heatmap, and processed them by Spearman’s correlation analysis and hierarchical clustering (Figure 1(c)). Altogether, 2107 genes were considered to be differentially expressed genes (DEGs) (|log2-fold change| ≥1 and FDR value <0.05) between the MD and MN groups. We considered 496 genes to be significantly downregulated and 1611 to be upregulated, and of these, 389 and 1151 were protein-coding genes, respectively (Figure 1(d–e) and Supplementary Tables S1 & S2). We also included some pseudogenes and lncRNA transcripts (Supplementary Tables S1 & S2).

Figure 1. Characteristics and differences between transcriptomes from human MD and MN oocytes. (a) Statistics of expressed transcripts, with the majority being protein-coding genes and pseudogenes. Long non-coding RNAs (lncRNAs), processed transcripts, antisense RNAs and some other transcripts were also detected. (b) Venn diagram indicating the numbers of overlapping genes between MD and MN groups. (c) Unsupervised hierarchical clustering based on transcriptome characteristics. (d) A volcano plot depicting significantly upregulated or downregulated genes between the 2 groups. (e) Compared with the MN group, 1611 and 496 transcripts were upregulated (shown in red) and downregulated (shown in green) in MD oocytes, respectively. Dark colors represent protein-coding genes and light colors represent non-coding genes.

Dramatic changes in expression of spliceosomal proteins in MD oocytes

To elucidate the potential impact of DEGs, we performed pathway enrichment analysis by referring to the GO, KEGG, Reactome, BioCyc, and some other repositories. Of the 64 GO terms that differed between MD and MN groups, 19/19/26 were with cellular components, molecular functions, and biological processes, respectively (Figure 2(a) and Supplementary Table S3). For the database-combined pathway enrichment analysis, 72 pathways (both upregulated and downregulated) were enriched, and we delineate the top 10 in Figure 2(b). Notably, spliceosome (involving 33 DEGs) was the pathway that was most significantly different between MD and MN oocytes.

Figure 2. Enriched GO terms and pathways based on DEGs. (a) Distribution of 64 GO terms according to DEGs. Twenty-six terms of biological process, 19 terms of cellular component, and 19 terms of molecular function are colored blue, green, and red, respectively. Full bars represent the number of DEGs in the recent study; dashed bars represent the numbers of total expressed background genes. Right vertical axis, the numbers of DEGs and background genes; left vertical axis, the percentage of annotated GO-term genes relative to background genes. (b) The top-10 significantly enriched pathway terms (both upregulated and downregulated) according to FDR values.

Pre-mRNA splicing – especially pre-mRNA alternative splicing (AS) – is significantly elaborated in eukaryotes and contributes to increases in the complexity of the proteome, introducing additional layers for regulating gene expression in different cell types and during development [27]. Composed of both proteins and RNAs, the spliceosome executes canonical pre-mRNA splicing by 2 sequential transesterification and splicing reactions (Figure 3(a)). In addition to the 130 genes collected in KEGG, we searched over 200 spliceosomal proteins recorded in the literature [28,29]. In 236 detected genes, 49 were differentially expressed (nearly 21%), with 40 genes upregulated and 9 genes downregulated (Figure 3(b) and Supplementary Table S4).

Figure 3. Pre-mRNA splicing and differentially expressed spliceosomal proteins. (a) Cartoon of a representative splicing cycle. (b) DEGs belong to both core snRNP sub-complexes and regulatory splicing factors. Upregulated genes are presented in red and downregulated genes in green. (c-d) Expression of LSm8 and BRR2 (c), and of SR and hnRNP genes (d). For a clearer display, the vertical axis was adjusted to log₂ FPKM, and the same adjustment was used below. Error bar, standard deviation. *, p <0.05. (e-f) Expression pattern of SR proteins (e) and hnRNPs (f) during early embryonic development. 2 C, 2-cell; 4 C, 4-cell; 8 C, 8-cell; M, morula; B, blastocyst. Expression values were calculated as RPKMs (reads per kilobase transcriptome per million reads), which were also used in the original literature.

Among core components in constitutive splicing, LSm8/NAA38 was expressed at a 4-fold lower level in MD oocytes (Figure 3(c)). LSm8 forms a hetero-heptameric ring with the other 6 LSm proteins in U6 snRNP, and prevents the degradation of the U6 small nuclear RNA (snRNA) by binding to its 3ʹ-end [30]. Since U6 snRNA is the principle catalytic component of the splicing reaction, reduced LSm8 expression likely disturbs the stability of U6 snRNA and impedes the splicing rate. The expression of BRR2 – also termed SNRNP200, the core component of U5 snRNP – was also significantly attenuated (4-fold) in MD oocytes. As a major regulator of the constitutive splicing reaction, BRR2 ATPase unwinds base pairs between U4/U6 snRNAs to promote the activation of the spliceosome [31]; thus, its downregulation would also eliminate the splicing process. In summary, the downregulation of constitutive spliceosomal components in the MD group indicated a dysregulation of pre-mRNA splicing at a global level.

The major regulators of AS are serine-arginine–rich (SR) proteins and heteronuclear ribonucleoproteins (hnRNPs). SR proteins bind to exonic and intronic splicing enhancers (ESEs and ISEs) to promote splice site selection, while hnRNPs recognize exonic and intronic splicing silencers (ESSs and ISSs) and function in an opposite manner [32]. The reciprocal antagonistic actions of SR proteins and hnRNPs regulate the pattern and efficiency of AS. In our MD oocytes, SRRM2, SRSF1/5/6/9/11, SREK1, and TRA2A were expressed aberrantly (Figure 3(d)). Interestingly, TRA2A and TRA2B (2 human homologs of Drosophila Tra) were expressed in an inverse fashion (Figure 3(d)). In Drosophila, Tra is required for spermatogenesis in the male germ line [33], while mammalian Tra2 proteins are sequence-specific splicing activators that participate in the control of cell-specific splicing patterns [34]. However, whether there are any differences between TRA2A and TRA2B remains to be determined. Our data indicated a potential interoperable role for TRA2A and TRA2B in splicing regulation in oogenesis. Therefore, to further explore the probable functions of these splicing regulators during early embryonic development, we assessed their expression tendencies based upon their published expression profiles [13] (Figure 3(e–f)). Surprisingly, the expression of both TRA2A and TRA2B increased dramatically from the 4-cell (4 C) to 8-cell stages (8 C) during human embryonic development (Figure 3(e)). This was also the case for hnRNPR (Figure 3(f)) – coincident with zygotic genomic activation (ZGA) – suggesting its regulatory role in embryogenesis. AS regulation has been reported to occupy a crucial function in spermatogenesis, although a role in oogenesis has not yet been reported. By discovering the association between spliceosomal components/events and oogenesis, our study provided novel insights into the mechanisms underlying developmental regulation, and opened new areas for further investigation.

Genes with altered AS patterns are specifically associated with cell cycle and metabolic processes

Protein diversity generated from differentiated AS (DAS) is highly correlated with cellular identity [35], and anomalous AS events always portend pathologic outcomes [36]. We thereby speculated that altered expression of spliceosomal proteins led to AS aberrance. Thus, we investigated this hypothesis by analyzing the DAS events between MD and MN oocytes based on exon-exon junction reads, and observed 12 detailed types of AS events across samples (Figure 4(a)). Five basic types, i.e., skipped exon (SE), mutually exclusive exon (MXE), alternative 3ʹ splice site (A3SS), alternative 5ʹ splice site (A5SS), and retention intron (RI), were counted and proportioned (Figure 4(b–c)). The 2329 genes were differentially spliced into more than 1 alternative isoform, suggesting that at least 1 DAS event occurred for each gene (Supplementary Table S5). Of these, 17 genes possessed more than 10 AS events (Figure 4(d)). DMKN (Dermokine) is known for its extensively differential splicing [37], and this is in agreement with our analysis. Three genes, FBXW12, GAS5, and MYL6, contained all 5 basic AS types (Figure 4(e)). The lncRNA gene GAS5 (growth arrest specific 5) contained 8 exons from which at least 15 alternatively spliced transcripts were identified, and functional studies have shown an emerging understanding of this gene in cell-cycle balancing and apoptotic regulation [38]. The DAS we observed in our study indicated that it might be associated with oocyte arrest as well.

Figure 4. Changed patterns of alternative splicing (AS) between normal and arrested oocytes. (a) Statistics of 12 detailed AS events for each sample. AE, alternative exon ends; XAE, approximate AE; SKIP, skipped exon; XSKIP, approximate SKIP; MSKIP, multi-exon SKIP; XMSKIP, approximate MSKIP; IR, intron retention; MIR, multi-IR; XMIR, approximate MIR; XIR, approximate IR; TSS, alternative transcriptional start site; TTS, alternative transcriptional terminal site. (b-c) Number (b) and percentage (c) of basic DAS events for MD and MN oocytes. SE, skipped exon; MXE, mutually exclusive exon; A3SS, alternative 3’splice site; A5SS, alternative 5ʹ splice site; RI, retention intron. (d) Genes with 10 or more AS event numbers. (e) FBXW12, GAS5, and MYL6 contained all 5 basic AS types.

To systematically summarize the impact of altered AS events on gene products and subsequent signaling, we conducted KEGG pathway enrichment analysis for the genes containing DAS events. Nine phenotypes/pathways were enriched, including Fanconi anemia pathway, cell cycle, mRNA surveillance pathway, ubiquitin-mediated proteolysis, protein processing in endoplasmic reticulum, amino sugar and nucleotide sugar metabolism, endocytosis, toll and IMD signaling pathway, and oocyte meiosis (Figure 5). The disorders in mRNA turnover and DNA/RNA synthesis may indicate that nucleotide metabolism is of great importance before oocyte maturation. Fanconi anemia pathway appears to be a master regulator of genomic integrity that repairs DNA inter-strand crosslinks (ICLs) in the genome, functions in genome maintenance, and coordinates several repair processes [39]. Owing to the evident enrichment of this pathway, we assumed that the genomic integrity of MD oocytes was disturbed to some extent, and that this was probably caused by irregular expression of spliceosomal proteins.

Figure 5. Enriched pathways of genes by DAS. Bubble chart indicating KEGG pathway enrichment for the genes featured by DAS and calculated with FDR values.

Discussion

Initiated in the embryonic stage, oogenesis consists of oocytogenesis, ootidogenesis, and final maturation so as to ultimately form an ovum. Primary oocytes remain arrested in the late diplotene (dictyate) stage of meiotic prophase I for up to 50 years, while transcription is completed and transcripts accumulate. The LH surge during each menstrual cycle then stimulates final oocyte maturation, and mature eggs are ovulated and become competent to be fertilized. This process is so complicated and time consuming that any disturbance would lead to maturation deficiency and infertility. Hereditary and/or environmental factors are currently suspected to lead to oocyte maturation failure, which reflects genomic and/or transcriptomic disorders. In the present study, by performing single human oocyte polyadenylated-RNA sequencing, gene expression profiling, and pathway enrichment analysis, we demonstrated for the first time that some fundamental genes – characterized by spliceosomal proteins and related factors – were differentially expressed between MD and normal MN oocytes; and that the AS modes of many genes were altered as well in MD oocytes.

The MD oocytes we analyzed in this study were from patients who showed complete oocyte maturation failure following standardized or even personalized clinical procedures. This failure in oocyte maturation may be caused by either hereditary or environmental factors, and we herein provided a possible explanation for aberrant gene expression and splicing. Nevertheless, limitations did exist for our study. First, even with the assistance of the most advanced single-cell sequencing, our sample size was relatively small. Extant challenges may be that the maturation block was rather specific in each of the 3 cases and might not be followed by similar expression patterns. This might explain why the RNA expression profile looked different between Y1 and Y2, Y3. Second, we did not retain any cumulus cells or additional oocytes for further experimental validation. Moreover, all the oocytes we saved were to be discarded and no MI oocytes were collected. Therefore, we could not exclude the possibility that the differences we observed might be caused by oocytes being at different stages (MI vs. MII), and this would be done in our following investigation.

Despite these limitations, our current study did provide some transcriptional insights into the process of oocyte meiotic maturation. Because transcription is quiescent and splicing events are comparatively “turned off” when the oocyte approaches the fully grown GV stage and during maturation [40], the aberrant expression of 2107 DEGs was probably from an accumulation at earlier steps such as primordial follicle activation. Among downregulated DEGs, many mitochondrial proteins – including cytochrome c oxidase subunit COX6B1 (subunit 6B1) and COX7C (subunit 7 C), ubiquinol-cytochrome-c reductase complex subunit UQCC2 (assembly factor 2), and ATP synthase complex subunit C14orf2 (6.8 kDa mitochondrial proteolipid) – were most prominently downregulated in MD oocytes (Supplementary Fig. S3A), indicating a high correlation between dysfunction of mitochondrial energy metabolism and abnormal oocyte arrest. Another dramatically decreased gene in MD oocytes was RPL36 (60S ribosomal protein L36, a component of the large ribosomal subunit [41]), whose expression level was only ~2% of that in MN oocytes (Supplementary Fig. S3B). In contrast, expression of the CBP80/20-dependent translation initiation factor CTIF was augmented (Supplementary Fig. S3B), representing abnormal functioning of ribosomal translation in oocyte maturation deficiency.

To expand the number of DEGs, we performed another round of enrichment, observing fold changes >1.5 (or <0.67) and FDR values <0.05 for upregulated and downregulated DEGs separately, and the spliceosome was enriched for both upregulated and downregulated DEGs (Supplementary Tables S6 and S7). Interestingly, other RNA processing pathways, including RNA transport, RNA degradation, and mRNA surveillance, were also enriched, indicating that these pathways were moderately regulated. Owing to the close connection between pre-mRNA splicing and RNA transport [42], we speculated that deficient transcriptional control – especially of pre-mRNA splicing – was intimately related to oocyte maturation incompetence.

Several genes have been previously reported to be responsible for oocyte maturation and fertilization. For example, PATL2 (protein PAT1 homolog 2) is an oocyte-specific translation repressor whose mutations cause oocyte maturation arrest at the GV stage [43,44]. Mutations in TUBB8 (tubulin beta 8 class VIII) impaired meiotic spindle assembly and were responsible for oocyte maturation arrest at the metaphase I (MI) stage [45]. Oocyte-specific kinase WEE2 has been identified as being mutated in the case of fertilization failure [46]. Deficiency in sperm-borne oocyte activation factor PLC_Z (phospholipase C zeta) or PAWP (post-acrosomal sheath WW domain-binding protein) accounted for oocyte activation deficiency [47]. Most recently, mutations in hemichannel protein PANX1 (pannexin1) were reported to cause oocyte death [48]. In our study, we demonstrated that transcriptional expression levels differed between MD and MN oocytes. However, without genomic examination, we were not able to conclude whether or not this phenomenon was caused by genetic mutations, and we are therefore currently pursuing this aspect further.

Recently, a case of unexpected repeat immature oocyte response was reported [49], similar to our cases, and herein we therefore tried to provide some clues for this type of primary infertility. Pre-mRNA splicing functions are fundamental to gene expression and regulation, but it is currently difficult to ascertain the primary mechanism(s) underlying oocyte maturation deficiency. Owing to the immense significance of DAS regulation in spermatogenesis, we hypothesize that it is also indispensable for oogenesis. Further investigations are thus needed to confirm how aberrant splicing is connected to oocyte maturation deficiency. Nevertheless, our data illustrate a significant aspect of gene regulation in oocyte maturation, and provide further opportunity for identifying biomarkers of oocyte quality.

Author contributions

This work was performed by all authors. J.H. and J.Q. conceived and designed the study; J.L. and M.L. collected the oocytes and performed the experiments; J.H. and P.Z. performed the data analysis and wrote the manuscript; and E.H., X.L., X.X., Z.W., and Y.F. contributed to sample collection and data analysis.

Disclosure statement

The authors assert that they have no potential conflicts of interest.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work received funding support from the under Grant Nos. 2018YFC1003702, 2016YFC1000601, and 2018YFC1004403; the National Natural Science Foundation of China under Grant Nos. 31800624 and 81501322; the Key Clinical Program of Peking University Third Hospital under Grant BYSY2017031; and the State Key Laboratory of Molecular Developmental Biology under Grant 2018-MDB-KF-16.