We collected transcriptome data of oocyte development during the period from germinal vesicle (GV) to metaphase II (MII) in different species, including donkeys, cattle, sheep, pigs, and mice. All RNA sequencing (RNA-seq) data was derived from public databases, including sequence read archive (SRA), gene expression omnibus (GEO), and genome sequence archive (GSA) with the following accession numbers: donkey (PRJNA763991) (Li Z. et al., 2021), cattle (CRA005589) (Li M.-H. et al., 2021), pig (GSE160334) (Du et al., 2021), sheep (GSE148022) (Wang J.-J. et al., 2020), and mouse (GSE119906) (Qian et al., 2019).

In order to ensure the accuracy of data analysis, we dealt with raw sequencing data, not processed files. Firstly, the FastQC (v0.11.8) was used to check the raw RNA-seq data and based on the quality control report, unqualified data was eliminated (Andrews, 2010). Next, Fastp (v0.23.1) was used for further quality control, and in this step, low-quality, unqualified reads were removed (Chen et al., 2018). As important sequence alignment software, STAR (v2.7.0f) was selected for sequence alignment to a reference genome, and we directly used the “--outSAMtype BAM SortedByCoordinate” parameter to generate a BAM format file (Dobin et al., 2013). The reference genome of each species was as follows: Equus asinus (assembly ASM1607732v2), Bos taurus (assembly ARS-UCD1.2), Ovis aries (assembly ARS-UI_Ramb_v2.0), Sus scrofa (assembly Sscrofa11.1), and Mus musculus (assembly GRCm38.p5). Finally, FeatureCounts (v1.6.3) was processed to generate gene counts (Liao et al., 2014).

The R package DESeq2 (v1.32.0) was adopted to identify the DEGs of oocytes during GV–MII (Love et al., 2014). The judgment threshold of significantly different DEGs was “|log2fold change| > 2 and p-value < 0.05”. In order to quantify the amount of gene expression, we used a custom R script to calculate fragments per kilobase of exon model per million mapped fragments (FPKM) (Zhao et al., 2021).

To obtain a comprehensive understanding of the uniqueness of donkey oocyte development during GV–MII, compared to other species, we used comparative transcriptomic analysis. Firstly, we used uniform conditions to obtain DEGs during GV–MII. Next, taking gene function annotations into consideration, the gene symbols in donkeys, cattle, sheep, and pigs, were converted into murine homologous gene symbol IDs via the R package gprofiler2 (v0.2.1) (Kolberg et al., 2020). Finally, we compared the similarities and differences between donkey and cattle, sheep, pig, and mouse oocyte development during GV–MII.

PCA of gene expression counts in this study were performed using the R package DESeq2 (v1.32.0), and the top 2 PCs were displayed by scatter plot (Love et al., 2014). The FPKM matrix of gene expression in the examined species, underwent hierarchical clustering analysis and the R package ggtree (v3.0.4) was used for visualization (Yu et al., 2017).

The R package WGCNA (v1.70-3) was used to uncover the correlation between genes (Langfelder and Horvath, 2008). First, FPKM was used to normalize gene expression levels among all species in the study. Next, in order to construct a mixed matrix of the gene expression of different species, the R package gprofiler2 (v0.2.1) was selected to unify the gene symbol IDs (Kolberg et al., 2020). Further, the function pickSoftThreshold() of R package WGCNA to ensure a scale-free network, the soft threshold of β was set to 16. A hierarchical clustering dendrogram, a heat map, and a topology overlap matrix (TOM) were used to show the relationship between the functional modules and genes. Cytoscape software (v3.8.2) was used to exhibit the top 10 topological overlap relationships in functional modules related to donkey oocyte development during GV–MII.

The hub genes, referring to highly interconnected nodes in functional modules, are considered as functionally important genes in WGCNA (Langfelder and Horvath, 2008). The hub genes were selected by the module membership and gene significance (MM & GS) method, and the conditional threshold was set as MM > 0.98 and SG > 0.8. Next, the hub genes were placed into the PPI network analysis through STRING (v11.5) (https://string-db.org/) to observe the interaction between genes.

DEGs were processed for GO and KEGG analysis using clusterProfiler (v4.0.5) (Wu et al., 2021). Considering that different function annotation libraries may lead to differences in results, Metascape, a web-based and timeously updated biological annotation database was also used for the same analysis of hub genes (Zhou et al., 2019).

Committed to obtaining key regulatory factors and signal pathways for donkey oocyte development, we designed the experimental program shown in Figure 1A. First, we collected the single cell RNA sequencing data of GV and MII phase oocytes in the studied species (the detail information in Supplementary Table S1). To achieve improved analysis results, we took a strict bioinformatics analysis process (as described in Materials and Methods). Next, we focused on the DEGs of the different species, and adopted a plan of pairwise comparison and overall comprehensive comparison. Moreover, we carried out WGCNA of all uniformed detected genes of oocyte development during GV–MII in different species (Figure 1A). Unsurprisingly, it uncovered that the transcriptome profile of oocyte development during GV–MII between different species was species-specific, and in the same species, the hierarchical clustering indicated that there was a clear difference between GV and MII stages (Figure 1B). Furthermore, in all species, the GV and MII stages showed a clear dividing line, which indicated that the transcription characteristics of these two periods were clearly separated (Figure 1C).

There is a broad consensus that DEGs tend to serve as key factors with the potential to regulate the transformation of cell fate (Stark et al., 2019; Zhang et al., 2019). Consequently, we collected detailed statistics of all the DEGs in oocyte development during GV–MII in different species. For the donkey, the number of down-regulation DEGs (down-DEGs) was 2,966, and the number of up-DEGs was 674 (Figure 2A and detail in Supplementary Table S2). For cattle, the number of down-DEGs was 1,526, and the number of up-DEGs was 990 (Figure 2B and detail in Supplementary Table S3). For sheep, the number of down-DEGs was 927, and the number of up-DEGs was 314 (Figure 2C and detail in Supplementary Table S4). For pigs, the number of down-DEGs was 1783, and the number of up-DEGs was 280 (Figure 2D and detail in Supplementary Table S5). For mice, the number of down-DEGs was 5,247, and the number of up-DEGs was 1,311 (Figure 2E and detail in Supplementary Table S6). In general, the number of down-DEGs was much greater than the number of up-DEGs; moreover, the number of DEGs detected in sheep was the least among all species (Figure 2F).

Organisms of different species have species-specific development (Yamashita et al., 2000); we may therefore ask whether oocyte development in different species also has species-specificity? To address this question, we made a detailed comparison of DEGs during oocyte development in donkeys and other species. For donkeys and cattle, when comparing the DEGs of oocyte development during GV–MII, we found a total of 402 identical DEGs; and for donkeys, there were 1,500 donkey specific DEGs (Figure 3A). On the 1,500 donkey specific DEGs, we performed GO analysis, and found that the top5 terms focused on non-coding RNA (ncRNA) metabolic process (Figure 3B). For donkeys and sheep, we found a total of 215 identical DEGs, and 1,687 donkey specific DEGs (Figure 3C). GO analysis of the 1,687 donkey specific DEGs found that the top5 terms focused on the ncRNA metabolic process (Figure 3D). For donkeys and pigs, we found a total of 298 identical DEGs, and 1,604 donkey specific DEGs (Figure 3E) of which GO analysis indicated that the top5 terms focused on the ncRNA metabolic process (Figure 3F). For donkeys and mice, we found 815 identical DEGs, with 1,087 of donkey specific DEGs (Figure 3G) of which GO analysis, indicated that the top5 terms included the homeostasis of number of cells (Figure 3H).

In order to better understand the intricacies of donkey oocyte development, we comprehensively analyzed the similarities and differences of DEGs in oocyte development between the studied species. For up-DEGs, compared with all other tested species, there were 257 donkey-specific genes (Figure 4A). Furthermore, results of GO and KEGG analysis suggested that the 257 donkey specific up-DEGs were involved in cell differentiation and immune process (Figures 4B,C). Interleukin 10 (IL10) is a cytokine involved in inflammation and immunosuppression, and is engaged in the regulation of cell growth and differentiation. Moreover, it participates in apoptosis (Dhingra et al., 2009). Interestingly, we found that IL-10 was a donkey specific up-DEG (Figure 4C). For down-DEGs, there were 555 donkey specific genes (Figure 4D). GO and KEGG analysis indicated that the 555 donkey specific down-DEGs participated in nucleoside biosynthetic metabolic process and cytosolic DNA-sensing pathway (Figures 4E,F). To our surprise, we found that C-X-C motif chemokine 10 (CXCL10) was a donkey specific down-DEG (Figure 4F). Because CXCL10 is an inflammatory cytokine and involved in apoptosis, it stands to reason that these donkey specific down-DEGs are also involved in apoptosis. Moreover, we constructed a list of DEGs of all detected species, including up-DEGs and down-DEGs, and used clusterProfiler’s compareCluster function to perform GO and KEGG function analysis on them. For GO analysis, compared with other species, the donkey was mainly reflected in the process of RNA metabolism regulation (Figure 4G). For KEGG analysis, not surprisingly, compared with other species, donkeys had significant enrichment in their ribosomal regulation process, which was well combined with the ncRNA and rRNA metabolism annotation in GO terms (Figure 4H).

WGCNA, increasingly used in bioinformatics analysis, is able to accurately pick out the hub genes related to traits of interest in a complex gene expression profile (Pei et al., 2017). In order to make the correlation between genes conform to the scale-free network distribution, the appropriate soft threshold must be set correctly. Here, we set the soft threshold at 16 to assure the downstream analysis of WGCNA (Supplementary Figure S1A). Next, 20 gene functional modules were identified and exhibited as a hierarchical clustering dendrogram and a module eigengene adjacency heat map (Supplementary Figure S1B). The turquoise module contained the greatest number of genes (Supplementary Figure S1B). Moreover, the heat map of module-sample relationships elucidated that the different functional modules were species-specific, and donkey-specific functional modules were yellow, blue, and light cyan (Figure 5A). The heat map of module-trait relationships also suggested that the different functional modules were species-specific, and donkey-specific functional modules were yellow, blue, and light cyan (Figure 5B). In order to identify the core functional modules that regulate the development of donkey oocytes, we constructed a hybrid matrix pool between the traits of donkeys and functional modules, and performed cluster analysis. The results showed that the yellow modules were closely related to the development of donkey oocytes (Figure 5C).

The TOM plot indicated the interactions between gene functional modules, and the results revealed an obvious strong interactive relationship between different functional modules, which was related to the species specificity between functional modules (Supplementary Figure S1C). In order to further explore the regulatory relationship between genes in the yellow functional module, first the GO analysis of genes in the yellow functional module showed that these genes were involved in meiotic cell division such as DNA repair, and apoptosis process such as transcriptional regulation by the TP53 pathway (Figure 6A). Moreover, the top 10 genes with the strongest topology overlap relationship were exhibited as a circle graph (Figure 6B). Notably, the hub genes were selected by the MM & GS method (Figure 6C); finally, we obtained a total of 36 hub genes. After observing the expression characteristics of these hub genes in different species, we found that these genes, including wee1-like protein kinase 2 (WEE2) that is a notable meiotic gene (Hanna et al., 2010), were highly expressed in donkey oocytes, compared to other species (Figure 6D). In order to systematically understand the interaction of these hub genes, protein-protein interaction (PPI) network analysis showed that the proteins translated by these hub genes did not directly interact strongly, which suggested that these hub genes were involved in different biological pathways to regulate the specificity of donkey oocyte development (Figure 6E). Furthermore, GO analysis showed that these hub genes were related to the negative regulation of cell development (Figure 6F).