Single-Cell Transcriptome Analysis Reveals Development-Specific Networks at Distinct Antral Follicle Stages in Sheep Oocytes

1. Introduction

Ovarian folliculogenesis is the complex maturation of an ovarian follicle through the primordial, primary, secondary, antral, and preovulatory developmental stages [1]. This developmental complexity is reflected in the high level of synchronization among the different cell types within one follicle that together produce a competent mature oocyte [2]. The generation of healthy oocytes relies on the coordinated development of somatic and germ cells in the ovarian follicle, and the rate of follicular development is determined by a developmental program intrinsic to the oocyte [3]. Activated preantral follicles develop to the antral stage, where follicle growth increases rapidly. The granulosa cells facilitate the accumulation of follicular fluid into a central cavity known as an antrum, giving rise to growing antral follicles [4]. Antral follicles develop in a wave-like pattern during which one or more follicles, “selected” from a group of small follicles (2–3 mm), develop into dominant follicles ≥ 6 mm in diameter. However, in mono-ovulatory bovine and ovine species, only one follicle becomes dominant, and subordinate follicles regress or degenerate [5,6].

The late antral follicle stage, when the oocyte acquires developmental potential, is critical for oocyte maturation [7]. Thus, the quality of an oocyte depends on its follicle, and follicular size correlates with developmental outcomes [8]. Oocytes from larger follicles are more likely to reach the blastocyst stage than those from small follicles [9], but regulation of gene expression is important in determining oocyte quality [10]. Thus, the relationship between the follicle stage and the oocyte transcriptomic landscape should be determined and correlated with developmental competence. The oocyte depends on stored maternal transcripts that play a key role in the development of oocyte meiotic competence during the growth phase in preantral and early antral follicles. Oocytes at the antral follicle stage resume meiosis [11] and support early development [12]. In addition, transcripts involved in energy production and metabolism are present in oocytes [13]. Differences in transcripts in the antral follicle development stages could explain subsequent oocyte maturation and identify biomarkers of oocyte quality. Oocyte development requires continuous intercommunication using highly coordinated paracrine, autocrine, and endocrine regulators, as well as gap junctions [14]. These oocyte-specific factors (OSFs) also promote proliferation of follicular somatic cells and steroidogenesis. Transcription factors (TFs) that are translocated from the cytoplasm to the nucleus mediate signaling and regulate the cellular phenotype [15], and the formation of follicles may be coordinated by a TF expressed by the oocyte. Thus, it is important to 1) identify oocyte-secreted TFs and other proteins involved in oocyte development, 2) classify these proteins by antral follicle stage, and 3) describe the dynamic associations of key factors between the follicle different stages.

Single-cell RNA sequencing (scRNA-seq) is a key tool for characterizing tissue heterogeneity at a cellular level [16], and it identifies changes in gene transcripts at the single-cell level that are undetectable by bulk RNA sequencing. Advances in transcriptome technology have also been applied to characterizing oocytes [17], although most oocyte transcriptome studies use animal models [18,19]. The scRNA-seq transcriptomes of human oocytes at five follicular stages revealed candidate secretory biomarkers of ovarian reserve [20]. The scRNA-seq transcriptome in germinal vesicles (GVs) and in in vitro metaphase II (MII) human oocytes identified maturation differences in oocytes [21,22]. However, the developmental pattern and transcriptional regulatory networks in sheep oocytes during antral follicle development, which differs from human and mouse oocytes [23], is not well characterized.

Here, we analyzed the transcriptome of single oocytes from antral follicles of different stages to identify dynamic changes in gene expression throughout antral follicle growth. We identified unique features of the transcriptional machinery associated with developmental-stage-specific expression patterns. We also identified stage-specific gene signatures in sheep oocytes, including oocyte-specific genes, oocyte-secreted factors, and TFs, which could provide oocyte development-specific profiles. We found differences in transcripts associated with specific biological processes and key signaling pathways that were associated with antral follicle development, providing information on expression patterns of hub genes regulating oocyte maturation. Thus, our characterization of the transcriptional landscape of antral follicle oocytes provides insights into the regulation of follicle growth in sheep.

2. Results

2.1. Overview of single-oocyte sequence data

To identify changes in oocyte transcripts during the antral follicle stage, we collected oocytes from antral follicles at different developmental stages, including small (≤ 3 mm, SO1, SO2, SO3), medium (3–6 mm, MO1, MO2, MO3) and large (≥ 6 mm, LO1, LO2, LO3) antral follicles for single-oocyte RNA-seq using the Smart-seq2 protocol (Figure 1A–C and Table S1). Following library construction and deep sequencing of oocyte transcriptomes, we obtained a total of 372,027,250 bases from the nine sequence libraries, producing a total of 334,979,082 (90.06 %) clean reads (Figure S1A and Table S2). The mapping efficiency of the clean reads to the sheep genome was high, with an average of 97.55 % for the small antral oocytes (SO), 97.15 % for the medium antral oocytes (MO), and 96.68 % for the large antral oocytes (LO), and a unique mapped rate average of 91.22 % (Figure S1B). The quality of these mRNA libraries was high enough to use them for statistical analyses. We identified an average of 15,341 expressed genes per oocyte, with a total of 138,073 expressed genes; 14,512 were expressed in the SO group, 15,300 in the MO group, and 16,212 in the LO group; however, the LO had the most mRNAs. We found that RNA synthesis was still detectable in the period following oocyte growth arrest and germinal vesicle breakdown (GVBD), and transcription ceased only at GVBD. Therefore, it is possible to resolve differences in mRNA expression levels between stages and identify the factors involved in oocyte developmental competence.

Using principal component analysis (PCA) for the groups of oocytes, we observed three distinct subpopulations, which corresponded to the morphological classification of antral follicular development (Figure 1D, E). Thus, antral follicles at different developmental stages exhibited different gene expression dynamics, as evidenced by distinct transcriptomes.

2.2. Gene expression dynamics and transcriptional characteristics of oocytes during antral follicle development

To determine the global expression profiles during antral follicle development, we characterized the 9,400 differentially expressed genes (DEGs) for the stages of follicular development in oocytes (false discovery rate [FDR, q-value] < 0.05; Figure 2A, B; Figure S2A). To characterize gene expression dynamics between consecutive stages in antral follicle development, we compared DEGs and found 2,141 genes overrepresented in the SO stage, 2,058 genes overrepresented in the MO stage, and 4,486 genes overrepresented in the LO stage. Unsupervised cluster analysis identified nine clusters of genes according to their expression patterns. The genes in cluster 5 (2,122) and cluster 7 (1,304) were underexpressed in the SO stage and were upregulated during antral follicular development, with a peak expression in the LO stage, suggesting that they indicate oocyte maturation. However, the transcripts for genes in cluster 9 (1,137) and cluster 6 (1,004) accumulated before follicle selection and dominance, suggesting that they represent the resting follicle pool. The gene expression patterns of cluster 2 (586), cluster 8 (792), and cluster 4 (678) showed expression at the MO stage (Figure S2B). However, because oocyte mRNA increased substantially during antral follicular development, the putative DEGs might simply reflect changes in the relative abundance of transcripts, with only some of them corresponding to changes in transcriptional state.

To identify transcripts specifically enriched at each developmental stage from the DEG expression profiles, we considered those genes with a log₂ transformed fold change (log₂FC) > |2| and FDR < 0.05 as differentially expressed genes. This yielded 1,299 genes overrepresented in the LO, 246 genes overrepresented in the SO stage (Figure 2C), and 565 genes identified in a comparison between the MO and SO stage (Figure 2D). Compared to the SO stage, 476 genes were upregulated in the MO stage, which accounted for the majority of DEGs, and only 89 genes were downregulated. Using the same criteria, we screened 820 genes, including 195 upregulated genes and 625 downregulated genes in the MO vs. LO oocytes (Figure 2E). The top 10 upregulated or downregulated genes in the SO vs. LO, MO vs. SO, and MO vs. LO are listed in Table S3. Using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of DEGs in oocytes at different stages to characterize their function in antral follicle growth, we identified several significantly enriched biological processes associated with each stage of follicular development (Figure 2F, G; Figure S2 C, D).

2.3. Gene expression signatures of oocytes at different antral follicle stages

Using Venn diagrams, we identified 365 genes that overlapped the 625 genes (MO vs. LO) and 1,299 genes (SO vs. LO) in the LO stage (Figure 3A), and 45 genes overlapped in the MO stage, whereas 33 genes overlapped in the SO stage. Surprisingly, there was no overlap in the three groups for upregulated genes. We clustered the overlapping genes with DEGs (SO vs. MO vs. LO; 9,400) and considered those genes with FDR < 0.05 as the stage-specific genes. We identified 222 stage-specific genes, including genes expressed in the oocytes from each antral follicular stage (Figure 3B and Table S4). Most of these genes (208) were upregulated during antral follicular development, with the highest expression in the large antral follicle stage. In addition, eight and six genes were highly expressed in the MO and SO stages, respectively.

We used GO analysis of genes specific to particular oocyte stages to determine the regulatory functions involved in antral follicle development (Figure 3C). Transcription coregulator binding (representative genes ESR1, JUND, MYC, NR5A2, and SMARCD3), glycolytic process (representative genes IER3, PFKL, PFKP, PGK1, and ENO1), vascular endothelial growth factor receptor signaling pathway (representative genes BCAR1, FGF10, FLT4, NRP1, PTK2B, VEGFA, and VEGFC), positive regulation of miRNA transcription (representative genes AR, JUN, MYC, and NGFR), response to hypoxia (DDIT4, EPO, HIF1A, LONP1, NOS2, P2RX2, PLAU, PLOD1, PLOD2, SOD3, and VEGFA), and cholesterol metabolic process (representative genes CH25H, CLN8, CYP11A1, LCAT, and SMPD1) were enriched in the LO stage, and GO enrichment items with FDR < 0.05 are shown (Figure 3D). Interestingly, regulation of calcium ion import (representative genes CACNG3), gap junction channel activity (representative genes GJB4), and steroid biosynthesis (representative genes HSD17B3) were enriched in the MO stage. We also identified candidate biomarkers uniquely expressed in the medium antral follicle stage: TENT5D, HGD, CRHR1, and CD86, which are involved in mRNA stabilization [24], tyrosine or L-phenylalanine metabolism [25], activation of G protein-coupled receptors [26], and adaptive immunity [27], respectively. In addition, ribosome biogenesis was enriched in the SO stage, with increased expression of ribosome-related genes (representative genes RPL23A and RPL24). GOLT1A plays a critical role in vesicle-mediated protein transport from the endoplasmic reticulum to the Golgi [28]. KBTBD6, as part of the Cul3 (KBTBD6/7)-RING ubiquitin ligase complex [29], is involved in the proteasome-mediated ubiquitin-dependent protein catabolic process and regulates RAC1 signal transduction and downstream biological processes, including the organization of the cytoskeleton, cell migration, and cell proliferation [30]. Members of the F-box protein family, such as FBXL14, are newly recognized E3 ligases that play a critical role in ubiquitylation [31]. EFHD2 is a conserved calcium-binding protein that plays a role in vitro in calcium signaling, apoptosis, actin cytoskeleton, and regulation of synapse formation [32]. Stage-specific gene expression profiles revealed a core hub of genes that regulate specific functions at each stage of oocyte development, including EFHD2 and FBXL14 in the SO stage, CRHR1 and CACNG3 in the MO stage, and PLAU, IL17RC, GPD1, and SKIL in the LO stage (Figure 3D). These signature genes may serve as potential cell-specific markers for each follicle stage. In addition, the subsets of 222 stage-specific signature genes in oocytes could serve as markers of oocyte developmental competency. Together, our results reveal the atlas of stage-specific transcriptional regulatory networks underlying antral follicle development in sheep.

2.4. Characterization of key pathways throughout antral follicle development

Gene Set Enrichment Analysis (GSEA) and KEGG analysis were used for pairwise comparisons of the follicular stages to identify changes in signaling pathways during the transition from the small stage to the medium stage during follicle recruitment (FDR<0.05). Ribosome and phosphatidylinositol signaling systems were upregulated in the SO stage. In contrast, 22 pathways were enriched in the MO stage, including ovarian steroidogenesis, fat digestion and absorption, glutathione metabolism, and neuroactive ligand-receptor interaction (Figure 4A). Notably, follicle selection is a critical step in oocyte maturation. We identified 25 functional pathways in the LO stage (FDR < 0.05), including glycolysis/gluconeogenesis, pentose phosphate pathway, biosynthesis of amino acids, apoptosis, and cellular senescence (Figure 4A). In addition, the nicotine addiction pathway and neuroactive ligand-receptor interaction were overrepresented in the MO stage. These intra-follicular metabolic signaling pathways likely mediate the transition from the medium to the large antral follicular stage, which involves the coordinated interaction of oocytes and granulosa cells for follicular development. Finally, the greatest increase in the number of signaling pathways occurred during the transition of oocytes from small to large antral follicles. Overrepresented pathways in the LO stage included ribosome, oxidative phosphorylation, aldosterone-regulated sodium reabsorption, and thermogenesis (Figure 4A).

To better understand whether these signaling pathways were enriched in sheep oocytes at specific antral follicle stages, we merged or overlapped the results described above to identify pathways that are specific or homogeneous during antral follicular development (Figure 4B). Ribosome genes were expressed predominantly in the SO stage with high expression of several ribosomal protein genes (representative genes RPL24, RPL27A, RPL36AL, and RPS12). Interestingly, the neuroactive ligand-receptor interaction signaling pathway (representative genes CRHR1 and OPRD1) was highly expressed in the MO stage. In addition, we identified 15 co-expressed signaling pathways at the LO stage (Figure 4C), and, among these, biosynthesis of amino acids and various types of N-glycan, arginine and proline metabolism, and protein processing in the endoplasmic reticulum play a key role in oocyte maturation. Because of the critical role of glycans in protein modification and the regulation of diverse cellular functions [33], we further analyzed the N-glycan biosynthesis pathway and found 12 highly expressed genes, including MGAT4B [34], ALG3 [35], and DDOST [36], suggesting that N-glycosylation was one of the key factors affecting sheep oocyte maturation.

2.5. Expression patterns of oocyte-specific genes at antral follicle developmental stages

We integrated our sheep oocyte transcriptomes with the human and mouse transcriptomes [37] to characterize the global expression profiles of oocyte-specific genes. To study the expression patterns of meiosis-related genes in oocytes, we obtained 38 genes involved in female meiosis I and nine genes involved in female meiosis II from the Human Protein Atlas (HPA) database and found 23 meiosis I genes and six meiosis II genes (FDR < 0.05) that were differentially expressed across different stages of oocyte maturation (Table S5). Based on the expression patterns, unsupervised cluster analysis identified 14 clusters of genes (Figure 5A). Most of these genes (including ESPL1, PSMC31P, PLK1, HSF2BP, MEIKIN, and UBB) showed upregulation as antral follicular development proceeds, with the highest expression at the large antral stage. Analysis of 44 genes from the Mouse Genome Informatics (MGI) and HPA databases involved in oocyte maturation revealed that some of these genes (including GRB14, TUT4, SIRT2, INSL3, FOXO3, and BMP15) increased gradually during oocyte development whereas others (such as REC8, FBXO5, WEE2, AURKA, and H3F3A) were expressed at high levels at the medium antral follicle stage (Figure 5B and Table S6). Among 62 maternal-effect genes, which play a critical role in mammalian development before zygote genome activation [38], 38 genes were expressed during antral follicle development (Table S7). The genes in cluster 5 (IHH, MBD2, LDOC11, and NLRP5), cluster 7 (ANGPT2 and KPNA6), cluster 8 (WT1), and cluster 11 (GET1, SDHD, and MELK) were upregulated during follicular maturation (Figure 5C). However, the genes in cluster 6 (HAVCR2 and PHLDA2), cluster 10 (PRDX3), cluster 12 (RGS2), and cluster 13 (BTG4) were expressed before follicle selection. Genes OOEP, AGRP, ZAR1, GRB10, UBE3A, LNPEP, NR2F2, BMPR2, and BICC1 were underexpressed in oocytes at all stages of follicular development. Among the maternal-to-zygotic transition genes, BTG4 and PABPN1L were relatively highly expressed, increasing in abundance during oocyte development (Table S8).

2.6. Identification of oocyte-secreted factors throughout antral follicle development

We identified 38 genes encoding sheep oocyte secretory proteins that were differentially expressed during antral follicular development (Table S9). The expression patterns of the genes from each cluster and their functions are shown in Figure 6A and 6B. Cluster 1 comprised nine genes expressed only in the LO stage, suggesting their involvement in the subsequent developmental competence of oocytes. Cluster 2 included nine genes expressed only in the MO stage, suggesting their involvement in the initiation of follicle selection. Cluster 3 consisted of nine genes that represented follicle recruitment. Protein-protein interaction (PPI) analysis of the oocyte-specific protein-coding genes showed that the PPI network was reconstructed from 18 gene symbols (Table S10). Among these, secretory factors GDF9, BMP15, NLRP5, ZP2, ZAR1, and Rec8 formed the core of the PPI network (Figure 6C, D; Table S11). These genes were related to oocyte development and were involved in progesterone-mediated oocyte maturation, the HIF-1 signaling pathway, and cytokine-cytokine receptor interactions. These oocyte-specific secretory factors reveal the diverse functions of oocytes in antral follicle development. Two important oocyte-secreted factorsare growth differentiation factor 9 (GDF9) and bone morphogenetic protein 15 (BMP15). They activate signaling pathways in cumulus cells (CCs) that regulate the differentiation and the distinctive phenotype of CCs. Also, the interaction of GDF9 and BMP15 with follicle-stimulating hormone (FSH) potentiates the expression of IGF-2 in human granulosa cells [39]. Nlrp5 is a maternal lethal effect gene in mammals, and loss of NLRP5 in mouse oocytes triggers premature activation of the mitochondrial pool, causing mitochondrial damage and apoptosis [40]. The zona pellucida surrounding all mammalian oocytes plays a vital role in fertilization and early development. Knockout of zona pellucida protein ZP2 in female mice produces oocytes without a ZP, resulting in fertilization failure and infertility [41]. Zygote arrest 1 (ZAR1) is the first oocyte-specific maternal-effect gene known to function in the oocyte-to-embryo transition [42]. Rec8-containing cohesin, bound to Smc3/Smc1a or Smc3/Smc1b, maintains bivalent cohesion in mammalian meiosis [43]. Sister chromatid cohesion mediated by the cohesin complex is essential for chromosome segregation in mitosis and meiosis [44]. Here, we found that the REC8 transcript decreased in oocytes from large antral follicles, which may indicate oocyte meiotic resumption. In addition, Wee1-like protein kinase 2 (WEE2), a key regulator of meiosis during prophase I, is involved in maintaining meiotic arrest [45]. Further evaluation of the expression levels of OSFs in the follicular fluid may identify novel biomarkers for assessing oocyte quality.

2.7. TF regulatory networks in the oocytes

To construct regulatory networks of TFs for the SO and LO stages, we obtained all 1,376 sheep TFs from the Animal Transcription Factor Database (TFDB v4.0), a comprehensive database for predicting and classifying genome-wide TFs. During the development from the small to large antral stages, we found 73 differentially expressed TF genes (Figure 7A); 63 TF genes were highly expressed in the LO stage, and the others were highly expressed in the SO stage. This suggests that these TFs play a critical role in the follicular cavity and oocyte maturation. During the transition from the small to the medium antral stage, 24 TF genes were upregulated, and five TF genes (NOTO, NR1I3, LMX1B, SOX2, and SNAPC4) were downregulated (Figure 7B). Genes for 37 TFs were overexpressed in the LO stage vs. the MO stage and, thus, are candidate regulators of the medium to large antral stage transition. Genes MIER3, PAX7, HELT, PPARD, and FOXJ2 were more highly expressed in the MO stage (Figure 7C). To identify the master regulators of sheep antral follicle development, we overlapped DEGs for TFs at different stages (Figure 7D). A large number of TFs were expressed only in the LO stage, indicating possible regulatory roles in the cytoplasmic and nuclear maturation of oocytes (Figure 7E). Searching for recurring instances of neighboring TFs with STRING provided a predicted PPI network for 50 TFs as hub genes for the LO stage. The cellular tumor antigen p53 (TP53) TF, along with Kruppel-like factor 6 (KLF6) and Jun proto-oncogene (JUN), are hub genes for the positive regulation of RNA polymerase genes (Figure 7F and Figure S3A; Table S12). Among 21 TFs that were co-expressed the MO and LO stages (Figure 7G and Figure S3B; Table S13), the gene for mothers against decapentaplegic homolog (SMAD4), the central mediator of TGF-β signaling, is indispensable for oocyte formation and development [46].

Next, we speculate that L3MBTL3, MYRF, SALL1, USF3, NOTO, NR1I3, LMX1B, and ZNF200 most likely initiate the transcription network in small antral follicles, whereas DLX4, ZNF169, and TBX15 could be the unique regulators of follicle selection. To determine the potential roles of TFs in regulating the expression of DEGs, a TF-gene regulatory network was established based on predicted stage-specific target genes. In the SO stage, key TFs (MYRF, SALL1, and LMX1B) showed a significant correlation with multiple DEGs (Figure 7H). Alkylglycerone phosphate synthase (AGPS), nuclear import carrier for Hsp70s (HIKESHI), and ectodermal-neural cortex 1(ENC1) encode proteins that promote oocyte stability through lipid metabolism [47], protein transport [48], and the cytoskeleton[49], respectively. In addition, we identified significant TF-target gene pairs ZNF169-ITGB3BP, ZNP169-SNRPD2, DLX4-ITGB3BP, and DLX4-SNRPD2 that may play a role in mRNA processing at the medium antral follicle stage (Figure 7I). Further characterization of these TFs will provide insights into the transcriptional control of oocyte-induced antral follicle development.

3. Discussion

Understanding the mechanisms underlying oocyte maturation is essential for understanding reproductive biology. From the GV stage, oocytes reduce transcription and rely on stored mRNA and the translational machinery to control gene expression [50]. Thus, to understand oocyte maturation, it is important to determine the transcriptional profile of genes that regulate this process. In this study, we removed the oocytes from single follicles and generated scRNA-seq data from oocytes at the small, medium, and large antral follicle stages for a comprehensive assessment of gene expression. We identified DEGs by pairwise comparisons of the three stages and found that the LO stage had more specifically expressed genes than the SO and MO stages. The oocyte is transcriptionally active when the antral follicle first forms, but stored transcripts are critical for the initial stages of embryonic development until transcription is activated. A total of 565 DEGs were found in the MO stage vs. the SO stage, representing protein transport, glycolytic processes, and regulation of mRNA stability. Circulating FSH supports the growth of a group of antral follicles until the largest follicle produces estradiol, which occurs in sheep when the largest follicle is ~6 mm. GO analysis revealed that the biological processes that were enriched in the LO vs. the MO stage included regulation of transcription, positive regulation of fibroblast proliferation, and response to hypoxia, which may indicate intrinsic differences in these follicle stages. These signature genes may serve as cell-specific markers for each follicle stage. Also, EFHD2 and FBXL14 genes in oocytes from the SO stage, CRHR1 and CACNG3 in oocytes from the MO stage, and the top 15 DEGs in oocytes from the LO stage could serve as markers of oocyte maturation and developmental competency and may provide essential genetic tools for stage-specific labeling.

Although the transition from the small to the large antral is a key step in oocyte development, the molecular mechanisms and signaling pathways that drive this transition are unknown. Here, we identified signaling pathways associated with oocyte development and maturation based on the expression patterns of maternal transcripts and genes for key pathways. We found that ribosome biogenesis was regulated in the small antral follicle stage, producing sufficient ribosomes to support early embryo development until the initiation of zygotic transcription [51]. Thus, the oocyte translation program plays a critical role in the preparation for oocyte development. Oocyte-granulosa cell communication, mediated by several ligand-receptor pairs, is essential for oocyte development [52], and ligand-receptor pairs transmit signals between granulosa cells and oocytes to coordinate oocyte-granulosa cell functions (e.g., regulation of transcription, translation, apoptosis, cell differentiation, and transport). We found that signaling between the oocyte and CCs is regulated by the neuroactive ligand-receptor interactions pathway. Moreover, the key active pathways in large antral follicles are primarily in protein synthesis and glucose metabolism, although we also found that the thyroid hormone signaling pathway was upregulated in the oocytes from large antral follicles. Interestingly, thyroid hormones may directly affect oocytes, as specific binding sites for thyroxin are found in mouse and human oocytes. Thyroid hormone status may also affect the estrogenic sensitivity of the ovary [53], which is required for follicular dominance. Our results indicate that regulation of mRNA translation is important in follicle recruitment and that the neuroactive ligand-receptor interaction pathway is critical for the activation of oocyte maturation.

Oocyte quality determines the ability of the oocyte to overcome meiotic arrest and progress to the MII stage [54]. Thus, we identified the gene expression signatures for meiosis-related and oocyte maturation-associated genes and found that MLH1 and AURKA, which are associated with spindle building and chromosome segregation [55], were among the genes highly expressed in medium antral follicles. Lodde et al. [56] demonstrated a temporal relationship between chromatin remodeling and the main morpho-functional events that characterize the final growth phase in bovine oocytes. MEIOC, together with YTHDC2, promotes a meiotic cell cycle program via post-transcriptional control [57]. In addition, changes in the expression of maternal-effect genes revealed different hub genes at each oocyte stage. The maternal-effect genes, which are the mRNAs produced and stored in the cytoplasm of growing oocytes [58], are the sole substrates for protein synthesis after transcription is silenced in oocytes. The transcript levels of ANGPT2, KPNA6, CBS, CRHBP, and BLCAP decreased sharply with antral follicle development. CBS is required for the acetylation of α-tubulin for proper spindle assembly during meiotic maturation, which affects oocyte quality [59]. Here, we provide evidence that meiotic and developmental competence is acquired gradually during the antral follicular process. This involves sequential events regulated by maternal transcripts within the oocyte, including transcriptional regulation, spindle assembly, and chromosomal recombination. Our high-quality data reveal extensive changes in gene expression associated with oocyte development at the single-cell level.

OSFs are critical for the regulation of oocyte maturation and follicular development, and the regulation of the ovarian follicular microenvironment through OSFs may be important for oocyte quality. Two oocyte-specific growth factors, growth differentiation factor 9 (GDF9) and bone morphogenetic protein 15 (BMP15), are critical players in the follicular development of sheep [60]. Interestingly, we found expression of the gene for natriuretic peptide receptor (NPR2) in the MO stage. NPR2 in sheep, unlike in mice and pigs, is not only in CCs but also in oocyte membranes, indicating that C-type natriuretic peptide (CNP) may activate intra-oocyte cGMP production via the membrane NPR2 in addition to the CC-mediated pathway [61]. In meiosis, Aurora kinase A (AURKA) regulates spindle building and chromosome segregation in oocytes [62], whereas thyroid hormone receptor interactor 13 (TRIP13), which is involved in recombination, is required early in meiotic recombination [63]. We showed here that meiotic resumption is activated at the medium antral follicle stage, whereas chromosome alignment and segregation occur at the large antral follicle stage, suggesting that oocytes from large antral follicles are competent to resume meiosis. These secretory protein-coding genes might serve as biomarkers for certain outcomes of interest.

TFs, which regulate mammalian development, bind to the regulatory regions of genes, forming a complex system that guides the expression of the genome [64]. We found expression of key TFs during antral follicle development, including Homeobox, helix-loop-helix (bHLH), and zf-C2H2, which are essential for transcriptional regulation of developmental processes [65,66]. We also determined that LMX1B, MYRF, and SALL1 play critical roles in the small antral follicle stage, whereas ZNF169 and DLX4 are critical in regulating follicle selection. TP53 and KLF6 serve as regulatory hubs for many TFs in large antral follicles. We found that the regulation of TFs is strongly associated with the follicle development stages, and we characterized the interactions between TFs and their target genes to identify the transcriptional regulators of DEGs. The functions of these TFs in oocyte development are only partially known; thus, understanding the role of these TF-DEG pairs in developmental regulation and the underlying molecular mechanisms requires further oocyte studies.

Our work revealed unique features in the transcriptional machinery and gene signatures of sheep oocytes at different antral follicular stages. These oocyte-specific genes and TFs may provide valuable clues for future functional studies. First, we identified the gene expression profiles for three antral follicle stages of oocytes and the functions of stage-specific transcribed genes. Second, we identified key regulatory networks and signaling pathways associated with oocyte development, including those for ribosome biogenesis and a neuroactive ligand receptor. Third, we identified oocyte-specific changes in gene expression that elucidated the regulation of oocyte maturation by maternal transcripts at different stages. Fourth, we identified subsets of oocyte-specific secretory factors that might serve as biomarkers for specific outcomes of interest. Fifth, we identified stage-specific transcriptional regulators and constructed a regulatory network of core transcriptional regulators and their target genes for the stages of antral follicle development. Together, these discoveries provide key information on oocyte development and broaden our understanding of sheep folliculogenesis.

4. Materials and Methods

4.1. Ovarian stimulation and oocyte collection

A total of 12 healthy adult Chinese Hu sheep (1.5–2.5-year-old), procured in April 2023 from the Inner Mongolia JinLai Animal Husbandry Technology Co., Ltd. (Hohhot, Inner Mongolia), were subjected to hormonal treatment as described by Bragança [67]. The animals received a progestogen-releasing device for 6 d before the initiation of superovulation by FSH soon after normal ovulation, i.e., during the emergence of the first follicular wave. We synchronized estrus in all ewes using the intravaginal progesterone device (CIDR, Pharmacia & Upjohn, Hartwell, Australia) for 6 days. One day prior to the removal of the CIDR, intramuscular injections of 300 IU eCG (Sansheng, Ningbo, China) and 0.12 mg cloprostenol sodium (Sansheng, Ningbo, China) were administered, and then 36 h after removing the CIDR 0.025 mg lecirelin (Sansheng, Ningbo, China) was administered. On the day of ovarian stimulation initiation (D0), each ewe was treated with a new intravaginal progesterone device (CIDR). Thereafter, the ewes received 360 IU of FSH (Sansheng, Ningbo, China) in six decreasing doses (75, 75, 60, 60, 45, and 45 IU) every 12 h (Figure 1A).

Oocyte retrieval was performed as described previously [68]. Briefly, cumulus-oocyte complexes (COCs) were recovered from each donor and classified according to antral follicle size: small (< 3 mm; representing early antral follicle stage), medium (3–6 mm, indicative of the follicle selection stage), and large (˃ 6 mm; indicating the follicular dominance stage). The oocytes and GCs were mechanically stripped by pipet several times. Subsequently, the naked oocytes were collected from small, medium, and large follicles and were labeled as small oocytes (SO), medium oocytes (MO), and large oocytes (LO). The oocytes from each stage were separated into different samples (three samples for each stage, SO1, SO2, SO3; MO1, MO2; MO3; LO1, LO2, LO3; 10–20 oocytes/each sample) and total RNA was extracted and purified as described previously [69].

4.2. Single-cell RNA sequencing and data analysis

Single-cell full-length cDNA libraries using Smart-seq2 were prepared as previously described [70] using an Illumina Novaseq^TM 6000 sequencing platform (Illumina, San Diego, CA, USA). Paired-end reads were processed using an RNA-seq pipeline. The pipeline first performed quality control on the raw reads using Cutadapt (v1.9) and FastQC (v0.11.9), then aligned the reads to the ovine genome sequence (Ovisaries.Oarv3.1) using the HISAT2 (v2.21) package. An average of 92.5% of the reads per sample were mapped to the reference genome. Next, the mapped reads of each sample were assembled by StringTie (v2.1.6) using the default parameters. All the transcriptomes from all samples were merged to reconstruct a comprehensive transcriptome using gffcompare software (v0.9.8). For the final transcriptome, StringTie and Ballgown were used to estimate the expression levels of all transcripts, and the abundance of the mRNAs was calculated from the fragment per kilobase mapped reads (FPKM).

4.3. PCA and identification of DEGs

We applied the Seurat method to analyze the single-cell data to create the whole clustering profile [71]. Only highly variable marker genes (coefficient of variation > 0.5) were used as inputs for PCA. The marker genes in the PCA plots were plotted by the FeaturePlot function in the Seurat package (v4.0.3). For better accuracy, we complemented the PCA oocyte clustering separately using the FactoMineR (v2.8) package in R (v4.1.0). The multiple t-test in the DESeq2 software was used to obtain the statistical significance of DEGs in each stage (one stage versus all other stages). Only genes with significant p-values, an FDR < 0.05, and a log₂FC transformed FPKM (fragments per kilobase of exon per million fragments mapped) larger than |2| were considered to be DEGs. The Enhanced Volcano R package was used to visualize the significance of each gene based on its log₂FC, and the top-ranked 20 markers for each development stage were plotted.

To identify unique gene expression functions and regulatory pathways of oocytes during antral follicle development, we performed GO and KEGG enrichment analysis using the DEGs generated by the R package ClusterProfiler [72]. Gene numbers for terms in GO and KEGG were calculated, and a hyper-geometric test was used to find significantly enriched GO terms (biological process) and KEGG pathways. We defined a term or pathway with an FDR < 0.05 as significantly enriched.

4.4. Identification of stage-specific DEGs

Only genes with (log₂FC) greater than |2|, an FDR < 0.05, and expressed in both groups were considered as stage-specific DEGs. GO analysis of stage-specific DEGs was performed with ClusterProfiler. Next, we analyzed stage-specific processes. GO enrichment items for the three samples were obtained based on GSEA (v4.1.0). We selected GSEA terms with p-values < 0.05 and listed the top-ranked 10 terms based on the normalized p-value enrichment score. We used GSEA and the KEGG pathway enrichment analysis to identify the significantly enriched pathways in oocytes in the transition from medium to large antral stages. The gene sets with FDR < 0.05 were considered to be enriched.

4.5. Identifying expression patterns of oocyte-specific genes from DEGs

To characterize the expression patterns of oocyte-specific genes and identify candidate markers for enhancing oocyte developmental competence in vitro, we focused on oocyte-specific functions. The meiosis-related genes and maternal-effect genes were downloaded from the HPA database (https://www.proteinatlas.org/). The genes involved in oocyte maturation were downloaded from the MGI database (https://www.informatics.jax.org/) and incorporated into the analysis. Sheep oocyte DEGs were overlapped with genes expressed in human or mouse oocytes to find sheep oocyte candidate genes. The average FPKM value of genes at each time point was imported into STEM software (v1.3.13; School of Computer Science, Carnegie Mellon University) to generate gene expression clusters.

4.6. Identifying OSFs at sequential antral follicle developmental stages

To identify candidate biomarkers for predicting oocyte development, we focused on 103 genes encoding oocyte-derived secretory proteins that we downloaded from the HPA database and selected those that overlapped the 9,400 sheep oocyte-secreted protein-coding DEGs. These candidate genes were analyzed using the online Search Tool for Recurring Instances of Neighboring Genes (STRING, https://cn.string-db.org/) to construct a PPI network. We set the cut-off as a combined score > 0.4 based on the number of associations with other genes in the PPI network and visualized the results using CytoScape software (v3.9.1).

4.7. Construction of transcriptional regulatory factor network in oocytes

We identified SO vs. LO stage DEGs (log₂(FC) > |2| and FDR < 0.05) and selected those that overlapped with 1,376 sheep TFs (AnimalTFDB, v4.0) as candidates TF to construct the network of transcriptional regulators and TF families. To identify the TFs associated with specific biological processes, we used STRING to visualize PPI interaction networks for key TFs. The active interaction sources for our search criteria were limited to Experiments, Databases, Co-expression, Neighborhood, Gene Fusion, and Co-occurrence, with a minimum interaction score set to the medium confidence value of 0.400. We also characterized the TF regulatory network and their differential target genes at different antral follicle stages. We used the Gene Transcription Regulation Database (GTRD) [73] and the Human Transcription Factor Targets Database (hTFtarget) [74] to find possible targets (within the 9,400 DEGs) that could bind the stage-specific TFs. Finally, the TF-target regulatory network was visualized using CytoScape software.

4.8. Statistical analysis

The experimental data were analyzed statistically using a two-tailed t-test to compare differences between different stages with PRISM software (GraphPad 6 Software). An FDR < 0.05 was considered statistically significant. The boxplot represents the median, the first quartile, and the third quartile values.

5. Conclusions

In summary, our results highlight the complexity and dynamic nature of oocyte development. We identified stage-specific genes, OSFs, and transcriptional regulators that contribute to the functional and morphological changes that occur during oocyte maturation. Overall, these findings provide insights into the molecular mechanisms that regulate oocyte maturation, with important implications for improving assisted reproductive technologies in domestic animals.