Comparative Transcriptomics Uncover the Uniqueness of Oocyte Development in the Donkey

The donkey is an important domestic animal, however the number of donkeys world-wide is currently declining. It is therefore important to protect their genetic resources and to elaborate the regulatory mechanisms of donkey reproduction, particularly, oocyte development. Here, we adopted comparative transcriptomic analysis and weighted gene co-expression network analysis (WGCNA) to uncover the uniqueness of donkey oocyte development compared to cattle, sheep, pigs, and mice, during the period from germinal vesicle (GV) to metaphase II (MII). Significantly, we selected 36 hub genes related to donkey oocyte development, including wee1-like protein kinase 2 (WEE2). Gene Ontology (GO) analysis suggested that these genes are involved in the negative regulation of cell development. Interestingly, we found that donkey specific differentially expressed genes (DEGs) were involved in RNA metabolism and apoptosis. Moreover, the results of WGCNA showed species-specific gene expression patterns. We conclude that, compared to other species, donkey oocytes express a large number of genes related to RNA metabolism to maintain normal oocyte development during the period from GV to MII.

Introduction

The donkey (Equus asinus) is a descendent of the African wild ass and is a common domestic beast of burden (Beja-Pereira et al., 2004). It can also provide meat and milk, and in particular donkey-hide gelatin, mostly composed of collagen, which is a traditional Chinese medicinal material (Kim et al., 2018). There is a broad consensus that, owing to the high content of linoleic acid, donkey meat is very palatable (Polidori et al., 2009). However, it is worrying that the number of donkeys is currently declining sharply owing to agricultural mechanization and the development of transportation vehicles.

Animal reproductive biotechnology, including in vitro fertilization and embryo transfer (IVF-ET) and its derivative technologies, is a perfect technology for saving endangered species (Bai et al., 2020). However, a limitation of this technology with donkeys is that the rate of in vitro oocyte maturation is relatively low (Goudet et al., 2016). Therefore, a deeper understanding of the characteristics of donkey oocyte development is needed, focusing on the meiotic maturation of oocytes, during the period from germinal vesicle (GV) to metaphase II (MII). The period from GV–MII marks the maturation of an oocyte, both nuclear and cytoplasmic (Eppig, 1996; Zhao et al., 2020). Studying the core regulatory factors of oocyte development during this period, and the precise and meticulous regulatory signaling pathways are of great significance for improving the in vitro maturation rate of donkey oocytes.

Previous studies mainly focused on economic traits, comparing donkeys with other livestock, but little research has addressed gene expression characteristics (Aganga et al., 2003; Polidori et al., 2008; Tian et al., 2020). Currently, research regarding single-cell RNA-seq (scRNA-seq) data of donkey oocytes suggests that the differentially expressed genes (DEGs) during GV–MII are related to the meiotic cell cycle, mitochondrial activity, and so on (Li Z. et al., 2021). The research suggested that wee1-like protein kinase 2 (WEE2) was involved in donkey oocyte development. Moreover, a study in mice indicated that Wee2-deficient caused fertilization failure and female infertility (Sang et al., 2018). WEE2 is indispensable for ensuring exit from meiosis in oocytes and promote pronuclear formation (Oh et al., 2011). However, the function of WEE2 is still unclear so far. In 2013, the donkey genome was been reported for the first time by Orlando et al., and in 2018, they released an updated version (Orlando et al., 2013; Renaud et al., 2018). In 2020, the Dezhou donkey genome was reported by Wang et al. (Wang C. et al., 2020). The reports of these reference genomes provide us with convenient conditions for studying the gene expression characteristics of donkey oocyte development.

RNA sequencing (RNA-seq) is able to determine the sum of all transcribed RNAs of tissues or cells under specific conditions, and has been widely used in animal and plant research (Stark et al., 2019). However, the limitation of this technology is that it requires large number of cells, which is not ideal for precious cell types such as oocytes. Excitingly, the emergence of scRNA-seq solves the problem of how to obtain the transcriptome status of a very small number of cells (Picelli et al., 2013). Here, we collected GV phase and MII phase scRNA-seq data from donkeys, cattle, sheep, pigs, and mice, and applied comparative transcriptomics analysis and weighted gene co-expression network analysis (WGCNA). Surprisingly, including WEE2, we identified 36 hub genes related to donkey oocyte development during GV–MII stages. This study provides novel information about key regulators of donkey oocyte development, which can provide a theoretical basis for the protection of donkey germplasm resources.

Materials and Methods

Data Collection

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).

Workflow of RNA-Seq Data Processing

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).

Differentially Expressed Genes Analysis

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).

Comparative Analysis of Oocyte Development Among Species

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.

Principal Component Analysis (PCA) and Hierarchical Clustering Analysis

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).

Weighted Gene Co-expression Network Analysis (WGCNA)

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.

Identification of Hub Genes and Protein-Protein Interaction Network Analysis

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.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Analysis

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).

Results

Overview of the Transcriptome Landscape of Oocytes in Different Species

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).

FIGURE 1

FIGURE 1. Transcriptome landscape of oocyte development in different species. (A) The experimental design flow chart of this study. (B) Hierarchical clustering of the transcriptome landscape of oocyte development in different species; different colors represent different species: red, donkey; steel blue, pig; green, mouse; gold, cattle; and purple, sheep (C) Principal component analysis (PCA) of different species; from top to bottom: donkey, cattle, sheep, pig, and mouse.

Identification of the DEGs in Oocytes From GV to MII in Different Species

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).

FIGURE 2

FIGURE 2. Overview of the differentially expressed genes (DEGs) during oocyte development in different species (A–E) The volcano plot exhibits the differential expression levels of oocytes during the period from germinal vesicle (GV) to metaphase II (MII) in different species, including donkey, cattle, sheep, pig, and mouse. (F) The bar plot shows the number of DEGs during the period from GV to MII in different species.

Difference of Oocyte Development in the Donkey Compared With Other Species

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).

FIGURE 3

FIGURE 3. Comparative analysis of donkey oocyte development during the period from germinal vesicle (GV) to metaphase II (MII). (A) The differences in oocyte development between donkeys and cattle during GV–MII. (B) The top5 GO (biological process) terms of donkey unique differentially expressed genes (DEGs) compared to cattle. (C) The differences in oocyte development between donkeys and sheep during GV–MII. (D) The top5 GO (biological process) terms of donkey unique DEGs compared to sheep. (E) The differences in oocyte development between donkeys and pigs during GV–MII. (F) The top5 GO (biological process) terms of donkey unique DEGs compared to pigs. (G) The differences in oocyte development between donkeys and mice during GV–MII. (H) The top5 GO (biological process) terms of donkey unique DEGs compared to mice.

The Uniqueness of Donkey Oocyte Development During GV–MII

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).

FIGURE 4

FIGURE 4. The unique transcriptional characteristics of donkey oocyte development during the period from germinal vesicle (GV) to metaphase II (MII). (A) Venn plot of the up-regulated differentially expressed genes (up-DEGs) in different species during GV–MII. The black arrow refers to the donkey-specific up-DEGs compared to other species. The area size represents the number of genes, and the area marked in blue contains the donkey-specific up-DEGs. (B and C) The GO (biological process) terms and KEGG annotation of donkey-specific up-DEGs, respectively. (D) Venn plot of the down-DEGs in different species during GV–MII. The black arrow refers to the donkey-specific down-DEGs compared to other species. The area size represents the number of genes, and the area marked in rose red indicates the donkey-specific down-DEGs. (E and F) The GO (biological process) terms and KEGG annotation of the donkey-specific down-DEGs, respectively. (G and H) The dot plot shows the GO (biological process) terms and KEGG annotation of the DEGs of donkeys and other species from GV–MII, respectively.

WGCNA of Oocyte Development in Different Species

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).

FIGURE 5

FIGURE 5. Weighted gene co-expression network analysis (WGCNA) of oocyte development during the period from germinal vesicle (GV) to metaphase II (MII) in different species. (A) The heat map shows the correlation between different samples and functional modules identified by WGCNA. (B) The heatmap shows the correlation between different phenotypes and functional modules identified by WGCNA. (C) The correlation between donkey and functional modules identified by WGCNA. The yellow functional module is the key functional module involved in donkey oocyte development.

Identification of Hub Genes Involved in Donkey Oocyte Development

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).

FIGURE 6

FIGURE 6. Functional exploration of hub genes related to donkey oocyte development identified by weighted gene co-expression network analysis (WGCNA). (A) The GO and KEGG annotation of genes in the yellow functional module. (B) The top 10 topological overlap relationship in the yellow functional module. (C) The scatter plot shows the relationship between module membership and gene significance (MM & GS). The large blue dots represent hub genes related to donkey oocyte development. (D) The heat map shows the expression levels of hub genes selected by the yellow functional module. (E) The protein-protein interaction network of hub genes. The connection between proteins represents the degree of protein interaction. (F) The GO (biological process) terms of hub genes related to donkey oocyte development.

Discussion

Donkey (Equus asinus), was domesticated around 3000 B.C. and became a traditional domestic animal whose ancestors are believed to include the African wild ass (Beja-Pereira et al., 2004; Rossel et al., 2008). With agricultural and industrial development and mechanization, the donkey has become somewhat redundant as a beast of burden and the global donkey population has notably reduced. China has the largest number of donkeys in Asia, estimated at more than three million in 2017 (Tian et al., 2020). However, as the demand for donkey-hide gelatin products increases, the number of donkeys in China is decreasing. There is now a need to address aspects of donkey breeding in order to preserve the species. Unfortunately, there is little research regarding the reproductive performance of donkeys (Tian et al., 2020; Li Z. et al., 2021). Hence, in this study, comparative transcriptomic and WGCNA were carried out to reveal the developmental characteristics of donkey oocytes, which would provide a theoretical basis for the protection of donkey germplasm.

The current study made an in-depth exploration of gene expression characteristics by thoroughly comparing the similarities and differences in oocyte development during GV–MII between donkeys and cattle, sheep, pigs, and mice (Qian et al., 2019; Wang et al., 2020b; Du et al., 2021; Li et al., 2021a; Li et al., 2021b), Furthermore, the release of the newly assembled donkey reference genome and the update of bioinformatics analysis tools (Wang C. et al., 2020), allowed us to annotate donkey genes more accurately, and to obtaining more refined results. Interestingly, in terms of transcriptome characteristics, we found that the GV and MII phases formed two groups of cells with a clear dividing line (Figure 1). Consistent with the results of previous studies, we found that in the studied species there were far fewer up-DEGs than down-DEGs (Figure 2). This phenomenon may be due to the gene transcription repression that occurs from the GV phase to the MII phase (De La Fuente et al., 2004); previous studies report that maternal mRNA is selectively degraded during development from GV to MII (Su et al., 2007; Zhao et al., 2020). Many studies suggested that non-coding RNAs (ncRNAs) were related to oocyte development (Barragán et al., 2017; Bouckenheimer et al., 2018). Interestingly, a study reported that ncRNAs were related to age and ovarian reserve in human (Barragán et al., 2017). But, the reports of ncRNAs in donkey oocyte development is few.

Unexpectedly, comparative transcriptomic analysis showed that donkey-specific DEGs are involved in the immune response, compared with other species. IL10 is an anti-inflammatory cytokine that is involved in apoptosis (Dhingra et al., 2009), and in this study, we found that it was a donkey specific up-DEG (Figure 4C). It has been reported in the literature that the expression of IL10 in granulosa cells is related to the success rate of IVF (Çavus and Deger, 2020). Moreover, the expression of IL10 attenuates the apoptosis rate of cells (Dhingra et al., 2009). However, few studies reported the role of IL10 in donkey oocytes. Interestingly, donkey specific down-DEGs, such as CXCL10 (Figure 4F), also participated in immune processes. It is reported that, CXCL10 serves as a marker of apoptosis (Sui et al., 2004), and the overexpression of CXCL10 increases the apoptosis of cells (Sui et al., 2006). The up-regulation of anti-apoptotic factors and the down-regulation of pro-apoptotic factors in donkey oocytes may be related to the maintenance of cell survival during oocyte development.

WGCNA was used to find co-expressed gene modules and to explore the relationship between gene networks and the phenotype of interest, as well as the core genes in the network; it is an effective method for analyzing the regulatory relationship of complex gene expression patterns (Pei et al., 2017). There is a consensus that the mechanism of oocyte development and maturation is species-specific (Yamashita et al., 2000). However, WGCNA data in the current study showed similar results for all studied species (Figure 5C). Subsequently, we explored the genes in the modules most related to donkeys and the results suggested that, consistent with the results of comparative transcriptomic analysis, these genes are closely related to the meiotic cell cycle and immune processes. Moreover, we selected 36 hub genes including WEE2, which was consistent with previous reports. Previous research reported that the overexpression of WEE2 is able to postpone the reentry of oocytes into meiosis in both mice and monkeys (Hanna et al., 2010). Furthermore, Wee2-deficient mice exhibit fertilization failure and female infertility (Sang et al., 2018).

Conclusion

Overall, our research brought new perspectives regarding the development of donkey oocytes. Compared with other species, donkey oocytes express a large number of genes related to RNA metabolism to maintain normal oocyte development during GV–MII. Our study offers a theoretical basis for improving oocyte maturation in the donkey.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Ethics Statement

The animal study was reviewed and approved by Qingdao Agricultural University.

Author Contributions

F-LZ and WS designed the study and WS managed funding. F-LZ did the bioinformatics. S-EZ and Y-JS collected the data, F-LZ and J-JW drafted the manuscript. The final manuscript approved by all authors.

Funding

This work was supported by the Start-up Fund for High-level Talents of Qingdao Agricultural University, Major Agricultural Application Technology Innovation Projects of Shandong Province (SD2019XM008), Modern Agricultural Industrial Technology System of Shandong Province (SDAIT-27), and Taishan Scholar Construction Foundation of Shandong Province (ts20190946) of China.

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.

The reviewer TZ declared a past co-authorship with several of the authors JJW, WS to the handling editor.

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.

Supplementary Material

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

Supplementary Figure S1 | (A) The evaluation of WGCNA genes of oocyte development during GV–MII in different species. Left panel, analysis of the scale-free fit index for various soft thresholding powers. The blue line represents a scale-free fit index of 0.9. Right panel, analysis of the mean connectivity for various soft-thresholding powers. (B) The hierarchical cluster tree shows functional modules distinguished by WGCNA. Bottom part, different colors represent different modules. (C) The heat map plot shows the topological overlap matrix (TOM) among 500 randomly selected genes. Light purple color depict a small overlap, and the red color indicates a greater overlap.

References

Aganga, A., Aganga, A., Thema, F., and Obocheleng, K. (2003). Carcass Analysis and Meat Composition of the Donkey. Pakistan J. Nutr. 2, 138–147. doi:10.3923/pjn.2003.138.147

CrossRef Full Text | Google Scholar

Andrews, S. (2010). “FastQC: A Quality Control Tool for High Throughput Sequence Data,” in Babraham Bioinformatics (Cambridge, United Kingdom: Babraham Institute).

Google Scholar

Bai, F., Wang, D. Y., Fan, Y. J., Qiu, J., Wang, L., Dai, Y., et al. (2020). Assisted Reproductive Technology Service Availability, Efficacy and Safety in mainland China: 2016. Hum. Reprod. 35, 446–452. doi:10.1093/humrep/dez245

PubMed Abstract | CrossRef Full Text | Google Scholar

Barragán, M., Pons, J., Ferrer-Vaquer, A., Cornet-Bartolomé, D., Schweitzer, A., Hubbard, J., et al. (2017). The Transcriptome of Human Oocytes Is Related to Age and Ovarian reserve. MHR: Basic Sci. Reprod. Med. 23, 535–548. doi:10.1093/molehr/gax033

CrossRef Full Text | Google Scholar

Beja-Pereira, A., England, P. R., Ferrand, N., Jordan, S., Bakhiet, A. O., Abdalla, M. A., et al. (2004). African Origins of the Domestic Donkey. Science 304, 1781. doi:10.1126/science.1096008

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouckenheimer, J., Fauque, P., Lecellier, C. H., Bruno, C., Commes, T., Lemaître, J. M., et al. (2018). Differential Long Non-Coding RNA Expression Profiles in Human Oocytes and Cumulus Cells. Sci. Rep. 8, 2202. doi:10.1038/s41598-018-20727-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Çavus, Y., and Deger, U. (2020). Investigation of ADAMTS-1, BAX and IL-10 Expressions in Granulosa Cells Determinanting Fertilization Success in IVF. Int. J. Morphol. 38, 427–434. doi:10.4067/S0717-95022020000200427

CrossRef Full Text | Google Scholar

Chen, S., Zhou, Y., Chen, Y., and Gu, J. (2018). Fastp: An Ultra-fast All-In-One FASTQ Preprocessor. Bioinformatics 34, i884–i890. doi:10.1093/bioinformatics/bty560

PubMed Abstract | CrossRef Full Text | Google Scholar

De La Fuente, R., Viveiros, M. M., Burns, K. H., Adashi, E. Y., Matzuk, M. M., and Eppig, J. J. (2004). Major Chromatin Remodeling in the Germinal Vesicle (GV) of Mammalian Oocytes Is Dispensable for Global Transcriptional Silencing but Required for Centromeric Heterochromatin Function. Develop. Biol. 275, 447–458. doi:10.1016/j.ydbio.2004.08.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhingra, S., Sharma, A. K., Arora, R. C., Slezak, J., and Singal, P. K. (2009). IL-10 Attenuates TNF- -Induced NF B Pathway Activation and Cardiomyocyte Apoptosis. Cardiovasc. Res. 82, 59–66. doi:10.1093/cvr/cvp040

PubMed Abstract | 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, 15–21. doi:10.1093/bioinformatics/bts635

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, Z. Q., Liang, H., Liu, X. M., Liu, Y. H., Wang, C., and Yang, C. X. (2021). Single Cell RNA-Seq Reveals Genes Vital to In Vitro Fertilized Embryos and Parthenotes in Pigs. Sci. Rep. 11, 14393. doi:10.1038/s41598-021-93904-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Eppig, J. (1996). Coordination of Nuclear and Cytoplasmic Oocyte Maturation in Eutherian Mammals. Reprod. Fertil. Dev. 8, 485–489. doi:10.1071/rd9960485

PubMed Abstract | CrossRef Full Text | Google Scholar

Goudet, G., Douet, C., Kaabouba-Escurier, A., Couty, I., Moros-Nicolás, C., Barrière, P., et al. (2016). Establishment of Conditions for Ovum Pick up and IVM of Jennies Oocytes toward the Setting up of Efficient IVF and In Vitro Embryos Culture Procedures in Donkey (Equus A). Theriogenology 86, 528–535. doi:10.1016/j.theriogenology.2016.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanna, C. B., Yao, S., Patta, M. C., Jensen, J. T., and Wu, X. (2010). WEE2 Is an Oocyte-Specific Meiosis Inhibitor in Rhesus Macaque Monkeys1. Biol. Reprod. 82, 1190–1197. doi:10.1095/biolreprod.109.081984

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J.-S., Kim, D., Kim, H.-J., and Jang, A. (2018). Protection Effect of Donkey Hide Gelatin Hydrolysates on UVB-Induced Photoaging of Human Skin Fibroblasts. Process Biochem. 67, 118–126. doi:10.1016/j.procbio.2018.02.004

CrossRef Full Text | Google Scholar

Kolberg, L., Raudvere, U., Kuzmin, I., Vilo, J., and Peterson, H. (2020). gprofiler2--an R Package for Gene List Functional Enrichment Analysis and Namespace Conversion Toolset G: Profiler. F1000Research 9, 709. doi:10.12688/f1000research.24956.1

CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: An R Package for Weighted Correlation Network Analysis. BMC bioinformatics 9, 559–613. doi:10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M.-H., Niu, M.-H., Feng, Y.-Q., Zhang, S.-E., Tang, S.-W., Wang, J.-J., et al. (2021a). Establishment of lncRNA-mRNA Network in Bovine Oocyte between Germinal Vesicle and Metaphase II Stage. Gene 791, 145716. doi:10.1016/j.gene.2021.145716

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Song, X., Yin, S., Yan, J., Lv, P., Shan, H., et al. (2021b). Single-Cell RNA-Seq Revealed the Gene Expression Pattern during the In Vitro Maturation of Donkey Oocytes. Genes 12, 1640. doi:10.3390/genes12101640

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, Y., Smyth, G. K., and Shi, W. (2014). FeatureCounts: An Efficient General Purpose Program for Assigning Sequence Reads to Genomic Features. Bioinformatics 30, 923–930. doi:10.1093/bioinformatics/btt656

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, 550. doi:10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Oh, J. S., Susor, A., and Conti, M. (2011). Protein Tyrosine Kinase Wee1B Is Essential for Metaphase II Exit in Mouse Oocytes. Science 332, 462–465. doi:10.1126/science.1199211

PubMed Abstract | CrossRef Full Text | Google Scholar

Orlando, L., Ginolhac, A., Zhang, G., Froese, D., Albrechtsen, A., Stiller, M., et al. (2013). Recalibrating Equus Evolution Using the Genome Sequence of an Early Middle Pleistocene Horse. Nature 499, 74–78. doi:10.1038/nature12323

PubMed Abstract | CrossRef Full Text | Google Scholar

Pei, G., Chen, L., and Zhang, W. (2017). WGCNA Application to Proteomic and Metabolomic Data Analysis. Methods Enzymol. 585, 135–158. doi:10.1016/bs.mie.2016.09.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Picelli, S., Björklund, Å. K., Faridani, O. R., Sagasser, S., Winberg, G., and Sandberg, R. (2013). Smart-seq2 for Sensitive Full-Length Transcriptome Profiling in Single Cells. Nat. Methods 10, 1096–1098. doi:10.1038/nmeth.2639

PubMed Abstract | CrossRef Full Text | Google Scholar

Polidori, P., Vincenzetti, S., Cavallucci, C., and Beghelli, D. (2008). Quality of Donkey Meat and Carcass Characteristics. Meat Sci. 80, 1222–1224. doi:10.1016/j.meatsci.2008.05.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Polidori, P., Cavallucci, C., Beghelli, D., and Vincenzetti, S. (2009). Physical and Chemical Characteristics of Donkey Meat from Martina Franca Breed. Meat Sci. 82, 469–471. doi:10.1016/j.meatsci.2009.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Qian, Y., Liao, J., Suen, A. H. C., Lee, A. W. T., Chung, H. S., Tang, N. L. S., et al. (2019). Comparative Analysis of Single-Cell Parallel Sequencing Approaches in Oocyte Application. Int. J. Biochem. Cel Biol. 107, 1–5. doi:10.1016/j.biocel.2018.12.003

CrossRef Full Text | Google Scholar

Renaud, G., Petersen, B., Seguin-Orlando, A., Bertelsen, M. F., Waller, A., Newton, R., et al. (2018). Improved De Novo Genomic Assembly for the Domestic Donkey. Sci. Adv. 4, eaaq0392. doi:10.1126/sciadv.aaq0392

PubMed Abstract | CrossRef Full Text | Google Scholar

Rossel, S., Marshall, F., Peters, J., Pilgram, T., Adams, M. D., and O'connor, D. (2008). Domestication of the Donkey: Timing, Processes, and Indicators. Proc. Natl. Acad. Sci. 105, 3715–3720. doi:10.1073/pnas.0709692105

CrossRef Full Text | Google Scholar

Sang, Q., Li, B., Kuang, Y., Wang, X., Zhang, Z., Chen, B., et al. (2018). Homozygous Mutations in WEE2 Cause Fertilization Failure and Female Infertility. Am. J. Hum. Genet. 102, 649–657. doi:10.1016/j.ajhg.2018.02.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Stark, R., Grzelak, M., and Hadfield, J. (2019). RNA Sequencing: The Teenage Years. Nat. Rev. Genet. 20, 631–656. doi:10.1038/s41576-019-0150-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, Y.-Q., Sugiura, K., Woo, Y., Wigglesworth, K., Kamdar, S., Affourtit, J., et al. (2007). Selective Degradation of Transcripts during Meiotic Maturation of Mouse Oocytes. Develop. Biol. 302, 104–117. doi:10.1016/j.ydbio.2006.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Sui, Y., Potula, R., Dhillon, N., Pinson, D., Li, S., Nath, A., et al. (2004). Neuronal Apoptosis Is Mediated by CXCL10 Overexpression in Simian Human Immunodeficiency Virus Encephalitis. Am. J. Pathol. 164, 1557–1566. doi:10.1016/s0002-9440(10)63714-5

CrossRef Full Text | Google Scholar

Sui, Y., Stehno-Bittel, L., Li, S., Loganathan, R., Dhillon, N. K., Pinson, D., et al. (2006). CXCL10-Induced Cell Death in Neurons: Role of Calcium Dysregulation. Eur. J. Neurosci. 23, 957–964. doi:10.1111/j.1460-9568.2006.04631.x

CrossRef Full Text | Google Scholar

Tian, F., Wang, J., Li, Y., Yang, C., Zhang, R., Wang, X., et al. (2020). Integrated Analysis of mRNA and miRNA in Testis and Cauda Epididymidis Reveals Candidate Molecular Markers Associated with Reproduction in Dezhou Donkey. Livestock Sci. 234, 103885. doi:10.1016/j.livsci.2019.103885

CrossRef Full Text | Google Scholar

Wang, C., Li, H., Guo, Y., Huang, J., Sun, Y., Min, J., et al. (2020a). Donkey Genomes Provide New Insights into Domestication and Selection for Coat Color. Nat. Commun. 11, 6014–6015. doi:10.1038/s41467-020-19813-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J.-J., Niu, M.-H., Zhang, T., Shen, W., and Cao, H.-G. (2020b). Genome-Wide Network of lncRNA-mRNA during Ovine Oocyte Development from Germinal Vesicle to Metaphase II In Vitro. Front. Physiol. 11, 1019. doi:10.3389/fphys.2020.01019

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, T., Hu, E., Xu, S., Chen, M., Guo, P., Dai, Z., et al. (2021). ClusterProfiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data. The Innovation 2, 100141. doi:10.1016/j.xinn.2021.100141

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamashita, M., Mita, K., Yoshida, N., and Kondo, T. (2000). Molecular Mechanisms of the Initiation of Oocyte Maturation: General and Species-Specific Aspects. Prog. Cel. Cycle Res. 4, 115–129. doi:10.1007/978-1-4615-4253-7_11

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Smith, D. K., Zhu, H., Guan, Y., and Lam, T. T. Y. (2017). Ggtree : An R Package for Visualization and Annotation of Phylogenetic Trees with Their Covariates and Other Associated Data. Methods Ecol. Evol. 8, 28–36. doi:10.1111/2041-210x.12628

CrossRef Full Text | Google Scholar

Zhang, F.-L., Li, N., Wang, H., Ma, J.-M., Shen, W., and Li, L. (2019). Zearalenone Exposure Induces the Apoptosis of Porcine Granulosa Cells and Changes Long Noncoding RNA Expression to Promote Antiapoptosis by Activating the JAK2-STAT3 Pathway. J. Agric. Food Chem. 67, 12117–12128. doi:10.1021/acs.jafc.9b05189

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, 475–511. doi:10.1186/s12864-020-06885-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Li, M.-C., Konaté, M. M., Chen, L., Das, B., Karlovich, C., et al. (2021). TPM, FPKM, or Normalized Counts? A Comparative Study of Quantification Measures for the Analysis of RNA-Seq Data from the NCI Patient-Derived Models Repository. J. Translational Med. 19, 1–15. doi:10.1186/s12967-021-02936-w

CrossRef Full Text | Google Scholar

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape Provides a Biologist-Oriented Resource for the Analysis of Systems-Level Datasets. Nat. Commun. 10, 1523–1610. doi:10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar