Nanopore-based full-length transcriptome sequencing of Muscovy duck (Cairina moschata) ovary.

Unlike mammals, studies on mechanisms that regulate waterfowl ovulation have been rarely reported. To advance our understanding of the ovulation differences in Muscovy duck, we utilized the Oxford Nanopore Technologies (ONT) to generate transcriptome data from 3 groups of female duck ovaries with ovulation differences (i.e., preovulation [PO], consecutive ovulation [CO], and inconsecutive ovulation [IO]). In this study, the full-length transcriptome data qualitative analysis showed that a total of 24,504 nonredundant full-length transcripts were generated, 19,060 new transcripts were discovered and 14,848 novel transcripts were successfully annotated. For the quantitative analysis, differentially expressed genes (DEGs) between the 3 groups were identified and functional properties were characterized. CTNNB1, IGF1, FOXO3, HSPA2, PTEN and SMC4 may be potential hub genes that regulate ovulation. Adhesion-related pathway, mTOR pathway, TGF-β signaling pathway and FoxO signaling pathway have been considered as important pathways that affect follicular development and ovulation. These results provide a more complete data source of full-length transcriptome for the further study of gene expression and genetics in Muscovy duck. The hub genes and potential mechanisms that affect the ovulation of Muscovy duck have been screened out to provide a scientific basis for breeding work to improve the reproduction performance of Muscovy duck.

INTRODUCTION

Muscovy duck (Cairina moschata) is a special waterfowl introduced from abroad to mainland China for purpose of the integration of meat, medicinal and ornamental (Baeza et al., 2002). It is highly suitable for intensive breeding and plays a vital role in the modern waterfowl industry (Aronal et al., 2002). However, female individuals have a strong broodiness, which results in lower reproductive performance and thus severely restricts the development of large-scale breeding. Therefore, improving the female Muscovy duck reproductive performance is one of the most critical approaches to accelerate the process for its industrial development.

The development of poultry follicles is an extremely complex and highly coordinated physiological process (Webb et al., 2003), which is similar to mammals, but has its own unique developmental characteristics. The development of poultry follicles strictly follows a certain level, from small to large and gradually mature until ovulation. The rank arrangement of mature follicles is the unique development characteristic of poultry follicles (Onagbesan et al., 2009; Johnson, 2012). At the peak period of ovulation, the largest follicle (F1) is selected from the follicle pool for subsequent ovulation (Ghanem and Johnson, 2019), but less than 1% of them can reach maturity and ovulation. Unselected follicles will move toward the fate of atresia, the selection of dominant follicles and ovulation will occur at the same time as the atresia of inferior follicles (Zhou et al., 2019).

In recent years, most transcriptome studies focused on the development of mammalian follicles. In contrast, researches on the follicles of waterfowl or poultry are rare, most of them were mainly focused on chickens ( Yin et al., 2020) and geese (Li et al., 2019b; Hu et al., 2020) by using second-generation sequencing technology (SGS). Although SGS is a major milestone breakthrough, it is clear that the short reads generated by the second-generation sequencing platform can no longer satisfy the need of the rapidly developing field of bioinformatics. Over a short period, the third-generation sequencing technology (single-molecule sequencing) has emerged. Full-length transcriptome sequencing based on Oxford Nanopore Technologies (ONT) can measure the change of ion current, according to the size of the current and the change of the current size, the base is judged through the complex algorithm of "Recurrent Neural Network". This process does not need to interrupt RNA fragments, and direct reverse transcription to obtain full-length cDNA, the ultralong read contains the sequence information of a single complete transcript. The transcript sequence can be directly used for further analysis (Jain et al., 2016; Bayega et al., 2018 and Magi et al., 2018). This platform can overcome the shortcoming (i.e., full-length transcripts cannot be accurately obtained) of second-generation sequencing (Koren et al., 2012). It can read from the beginning without assembly, and obtain more transcript information, which is more conducive to comprehensive analysis. Increasing studies utilized ONT sequencers for transcriptome sequencing, such as Arabidopsis thaliana (Cui et al., 2020) and Saccharomyces cerevisiae (Jenjaroenpun et al., 2018).

Intending to gain more insight into the developmental mechanism of Muscovy duck follicles, ONT sequencing technology was used to obtain full-length transcriptome of the follicles in Muscovy duck, full-length transcriptome provide a more complete data source for the further study of gene expression and genetics in Muscovy duck. Especially, we analysis the ovaries transcriptome with ovulation differences, revealed the key pathways and hub genes for ovulation in Muscovy duck follicles. These results are expected to provide new perspective for the molecular mechanism regulating the ovarian development of Muscovy duck even in poultry.

MATERIALS AND METHODS

Ethics Statement

This study complied with institutional and national guidelines for the care and use of animals. It was approved by the Institutional Animal Care and Use Committee at the Northwest A&F University. All efforts were made to minimize animal suffering.

Animals and Sample Collection

Muscovy duck were obtained from Shaanxi Anda Agricultural Development Co., Ltd (Shaanxi, China). We sampled the ovaries of a total of 9 ducks, including 3 ducks that had not ovulation (22 wk) and 6 that were at the peak of egg production (40 wk), in which 3 were consecutive ovulations ducks, 3 were inconsecutive ovulations ducks (i.e., preovulation [PO], consecutive ovulation [CO], and inconsecutive ovulation [IO]). All selected ducks had identical genetic background and appearance, raised under the same feeding management condition, and ducks in each group had a similar body weight. The female duck ovaries with ovulation differences were pooled for 3 biological replicates, respectively. The total RNA was extracted from each ovary sample using TRIzol reagent (Invitrogen, Paisley, UK) according to the manufacturer's protocol. RNA quality was evaluated using the Nanodrop-2000 (Thermo Fisher Scientific, Waltham, MA) and Agilent Bioanalyzer 2100 (Agilent, Santa Clara, CA) by running electrophoresis with 1% agarose gel. After RNA integrity and concentration meet the sequencing requirements, we use it for subsequent analysis.

Oxford Nanopore Technologies Long Read Processing

The ONT processing workflow is summarized in Figure 1. One ug total RNA was prepared for cDNA libraries using cDNA-PCR Sequencing Kit (SQK-PCS109) protocol provided by Oxford Nanopore Technologies (ONT). Briefly, the template switching activity of reverse transcriptases enriches for full-length cDNAs and adds defined PCR adapters directly to both ends of the first-strand cDNA; And following cDNA PCR for 14 circles with LongAmp Tag (NEB). The PCR products were then subjected to ONT adaptor ligation using T4 DNA ligase (NEB). Agencourt XP beads were used for DNA purification according to ONT protocol. The final cDNA libraries were added to FLO-MIN109 flowcells and run on PromethION platform at Biomarker Technology Company (Beijing, China). The original fastq data were filtered through short fragments and low-quality reads (length less than 500 bp, Qscore less than 7), and the total clean data was obtained. The full-length sequence was compared with the reference genome (Anas platyrhynchos PBH1.5) using minimap2 software (Li, 2018), and after clustering by the comparison information, the consensus sequence was obtained using pinfish software. Multiple copies of the same transcript may not be concentrated in the same consensus sequence, resulting in redundant sequences. The consensus sequence of each sample was merged, compared with the reference genome through minimap2, and the comparison results were de-redundant. Sequences with identity less than 0.9 and coverage less than 0.85 were filtered, and 24,504 nonredundant transcript sequences were finally obtained. All raw reads of transcriptome sequences have been submitted to the NCBI Sequence Read Archive (SRA) with identifier PO1: SRR12836353, PO2: SRR12836352, PO3: SRR12836351, IO1: SRR12836350, IO2: SRR12836349, IO3: SRR12836357, CO1:SRR12836356;CO2:SRR12836355;CO3:SRR12836354.

Figure 1

The flow chart illustrating the steps of nanopore-based RNA sequencing.

Structure Analysis and Gene Functional Annotation

Transcripts were validated against known reference transcript annotations with gffcompare (gffcompare: classify, merge, tracking and annotation of GFF files by comparing to a reference annotation GFF.). Simple Sequence Repeat (SSR) of the transcriptome was identified using MISA (Thiel et al., 2003). CDS were predicted by TransDecoder (https://github.com/TransDecoder/TransDecoder). Gene function was annotated based on the following databases: the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2004), nonredundant protein sequences (NR) ( Deng et al., 2006), Gene Ontology (GO) (Ashburner et al., 2000), Protein Family (Pfam) (McKenna et al., 2010), Eukaryotic Orthologous Groups/Clusters of Orthologous Groups (KOG/COG) (Tatusov et al., 2000; Koonin et al., 2004), and reviewed protein sequence database (Swiss-Prot) ( Rolf et al., 2004).

Transcription Factors Prediction and lncRNA Analysis

Animal transcription factors were identified from Animal TFDB (Hu et al., 2019). Four computational approaches include CPC (Kong et al., 2007), CNCI (Sun et al., 2013), CPAT (Wang et al., 2013) and Pfam (Finn et al., 2014) were combined to screen out nonprotein coding RNA candidates from putative protein-coding RNAs in the transcripts. Putative protein-coding RNAs were filtered out using a minimum length and exon number threshold. Transcripts with lengths more than 200 nt and more than 2 exons were selected as lncRNA candidates and further screened using CPC/CNCI/CPAT/Pfam that can distinguish the protein-coding genes from the noncoding genes.

Differentially Expressed Gene Analysis

Full length clean reads were mapped to the reference transcriptome sequence. Reads with match quality above 5 were further used to quantitative analysis. Expression levels were estimated by reads per gene/transcript per 10,000 reads mapped. Differential expression analysis of 2 conditions/groups was performed using the DESeq R package (1.18.0). DESeq provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution (Anders and Huber, 2010). The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with a P-value <0.01 and fold change ≥1.5 found by DESeq were assigned as differentially expressed.

Protein-Protein Interaction (PPI)

Align the target gene with the protein in the database to find the homologous protein, and construct an interaction network based on the interaction relationship of the homologous protein (the protein-protein interaction of which exists in the STRING (Franceschini et al., 2013) database: http://string-db.org/) Then the PPI of these DEGs were visualized in Cytoscape (Shannon et al., 2003).

Quantitative Real-Time PCR (qRT-PCR)

Reactions were conducted under the following conditions: predenaturation at 95°C for 5 min, followed by 40 cycles of denaturation at 95°C for 15 s and annealing/extension at corresponding temperature of each primer set for 30 s. The no-template controls and negative controls without reverse transcriptase were also included in all qPCR runs. All dissociation curves exhibited a single peak, indicating that all the primers were suitable for qRT-PCR. The qRT-PCR primers of selected genes are listed in Supplementary Table 1. The threshold cycle (Ct) was determined using the default settings, all data were analyzed using the 2−△△Ct method.

Statistical Analysis

The data were analyzed with one-way ANOVA using the GraphPad Prism 8.0 software. The data were presented as mean ± SEMA value of P < 0.05 was considered statistically significant.

RESULTS

Overview the Full-Length Sequences

In this study, samples were sequenced for full-length transcriptome. The clean data reached 6.02 GB. The number of full-length sequences obtained for each sample ranges from 4,839,392 to 7,213,668 (Supplementary Table 2). The consensus isoform is compared to the reference genome by the software minimap2 and then removed redundancy. Finally, 24,504 nonredundant transcript sequences were obtained. The deredundancy of all samples was compared with the known annotations of the reference genome, and 19,060 novel transcripts were found.

Structure Analysis and Functional Annotation of Novel Transcript

All the obtained novel transcript sequences are aligned with 7 databases: NR, Pfam, KOG, COG, Swiss-Prot, KEGG and GO. In total, 14,796 isoforms were annotated; there were 12,177 isoforms annotated in the GO database and 10,256 isoforms annotated in the KEGG database. The NR database had the highest number of isoform annotations (14,796 isoforms), whereas the COG database (3,671 isoforms) had the lowest isoform number (Table 2).

Table 2

Information of function annotation.

Annotated databases	New isoform number
COG	3,671
GO	12,177
KEGG	10,256
KOG	10,718
Pfam	10,386
Swiss-Prot	9,085
eggNOG	14,044
nr	14,796
All	14,848

SSR and TF Analysis

A total number of 20,540 sequences with a total size of 31,469,080 bp were subjected to SSR prediction. As a result, the total number of identified SSRs is 8,040, and the number of SSRs containing sequences is 5,797 (Table 3). Seven types of SSRs were identified: Mononucleotide, dinucleotide, Tri-nucleotide, Tetra-nucleotide, Penta-nucleotide, Hexa-nucleotide, and compound SSR (Figure 2). SSR was analyzed using MISA software. Animal transcription factor identification predicts a total of 3,235 transcription factors from the new transcripts obtained (Figure 3A).

Table 3

Type and number of SSRs identified.

Searching item	Numbers
Total number of sequences examined	20,540
Total size of examined sequences (bp)	31,469,080
Total number of identified SSRs	8,040
Number of SSR containing sequences	5,797
Number of sequences containing more than 1 SSR	1,613
Number of SSRs present in compound formation	701
Mono nucleotide	5,588
Di nucleotide	749
Tri nucleotide	1,554
Tetra nucleotide	85
Penta nucleotide	51
Hexa nucleotide	13

Figure 2

Distribution of SSR type. Abbreviations: C, compound SSR; p1, perfect Mono-nucleotide SSR; p2, perfect Di-nucleotide SSR; p3, perfect Tri-nucleotide SSR; p4, perfect Tetra-nucleotide SSR; p5, perfect Penta-nucleotide SSR; p6, perfect Hexa-nucleotide SSR.

Figure 3

The predictive analysis of novel transcript. (A) Number and family of transcription factors predicted. (B) Length distribution of the coding sequence of complete ORFs. (C) Venn diagram of lncRNA transcripts identified from CPC, CNCI, CPAT, and Pfam. (D) The type of lncRNA.

CDS and lncRNA Prediction

A total of 15,810 ORFs were predicted, including 7,903 complete ORFs. The predicted CDS results are shown in Figure 3B. A total of 1,815 lncRNA transcripts were predicted by CPC/CNCI/CPAT/Pfam analysis. The noncoding transcripts identified by the 4 methods were intersected with the results of lncRNAs for subgroup analysis. All detected lncRNAs were subdivided into the following 4 types: lincRNAs (long intergenic noncoding RNAs, 1,353, 74.5%), antisense-lncRNAs (186, 10.2%), intronic-lncRNAs (179, 9.9%), and sense-lncRNAs (97, 5.3%) (Figures 3C and D).

Identification and Functional Profiles of Differentially Expressed Genes

We identified genes with a Fold Change≥1.5 and P-value <0.01 in a comparison as significant DEGs. The number of DEGs in the IO-CO group comparison was 179, including 53 up-regulated genes and 126 down-regulated genes, while a total of 3,046 DEGs (1,046 up-regulated and 2,000 down-regulated) were identified in the PO compared with the CO group (Table 4). Volcano Plot can be used to more intuitively view the differences in the expression levels of DEGs in different groups and the statistical significance of the differences (Figures 4A and B). As shown in Figure 4C, there were 135 DEGs present between pairwise comparisons. Six DEGs were randomly selected to verify using qRT-PCR (Supplementary Figure 1). Regardless of the fold change, the expression of almost all of these selected mRNAs determined by qRT-PCR showed a similar trend to the results of RNA-seq, indicating the true reliability of our ONT sequencing methods.

Table 4

Number of different expressed genes in three groups.

DEG set	DEG numbers	Upregulated	Downregulated
IO1_IO2_IO3_vs _CO1_CO2_CO3	179	53	126
PO1_PO2_PO3_vs _CO1_CO2_CO3	3,046	1,046	2,000

Figure 4

Identification and functional analysis of differentially expressed genes. (A, B) Volcano plots showing significantly up and down regulated genes between PO vs. CO and IO vs. CO, respectively. (C) Venn diagram indicating the number of DEGs between three groups. (D, E) Annotation of Gene Ontology (GO) functional between PO vs. CO and IO vs. CO, respectively. Red represents BP, green represents CC and blue represents MF.

Gene Ontology (GO) Analysis

The GO database describes the functional properties of genes and gene products in an organism. The database has 3 broad categories: molecular function (MF), cellular component (CC), and biological processes (BP), each describing the molecular function that a gene product might perform, as well as the cellular environment and the biological processes involved. By annotating the 1,737 DEGs into GO database, we found that 1,131 DEGs were classified into biological process, 1,226 DEGs were classified into cellular component and 986 DEGs were classified into molecular function (Figure 4D). The DEGs most significantly enriched in the CC category in PO-CO group were cell, followed by cell part and organelle, while the MF of the GO category encompassed binding, catalytic activity, molecular transducer activity, signal transducer activity, transporter activity, nucleic acid binding transcription factor activity. Cellular process, single-organism process, biological regulation, metabolic process, and response to stimulus category were enriched among the top 5 GO terms in the BP. Further, there were 68 DEGs significantly enriched in GO terms between IO-CO groups (Figure 4E), the biological process (48 genes) was the most significantly enriched, followed by cell composition (47 genes) and molecular function (41 genes).

KEGG Pathway Enrichment Analysis

To further identify the major biochemical, metabolic, and signal transduction pathways of the DEGs, we performed a KEGG pathway enrichment analysis. Comparison of PO and CO, there were 1,466 DEGs were mapped to the reference pathways recorded in the KEGG database, yielded a total of 57 DEGs mapped to known KEGG pathways between IO-CO group. The top 20 pathways with the lowest Q values were shown in Figures 5A and B. DEGs from PO-CO group were mainly distributed in the Ribosome, Focal adhesion, Regulation of actin cytoskeleton, Oxidative phosphorylation and ECM-receptor interaction pathways, while IO-CO group focus on Focal adhesion, ECM-receptor interaction and Phagosome pathways.

Figure 5

Top 20 KEGG pathways enriched by DEGs between PO vs. CO (A) and IO vs. CO (B).

DEG Protein-Protein Interaction Analysis and Hub Gene Identification

To further explore the interactions among DEGs, a novel interaction network was constructed using Cytoscape software. Network Analysis can quickly and intuitively understand the complex gene regulatory network of follicle development. DEGs were selected in reproduction related-GO term were selected (Supplementary Table 3), then the Protein-Protein Interaction network was constructed. Among them, nodes with the greatest number of neighbors were considered as hub genes, CTNNB1, IGF1, FOXO3, HSPA2, PTEN, and SMC4 being the most significant node genes (Figure 6).

Figure 6

Protein-Protein Interaction analysis. The network shows genes up-regulated (yellow) and down-regulated (blue).

DISCUSSION

Egg production is one of the most important indicators that affect the economic benefits of the poultry industry. High-yielding Muscovy duck will ovulate continuously every day during the peak period of egg production to achieve higher egg production. In this study, we conducted ONT technology to obtain complete transcriptome of Muscovy duck ovary. Finally, 24,504 transcript sequences were obtained, of which 19,060 new transcripts were discovered. A total of 8,040 SSRs were obtained after SSR prediction. Analysis of the sequence structure of the newly discovered transcripts, a total of 7,903 complete ORF sequences and 1,815 lncRNAs were predicted, and 14,848 new transcripts were functionally annotated. Transcription factor prediction was performed on all transcripts, and a total of 3,235 transcription factors were predicted. The linkage relationship between microsatellite DNA and trait genes can be used to improve certain economic traits. In the breeding process, low-heritability quantitative traits can be targeted to improve according to genotypes (Bakhtiarizadeh et al., 2012). In the future, we will use the transcriptome of the Muscovy duck; highly polymorphic SSRs were screened to provide theoretical basis for genetic diversity research, germplasm resource protection and genetic selection, and breeding of Muscovy duck. The current research on lncRNAs related to animal reproduction shows that their mechanism of regulation of important economic traits such as livestock reproduction is involved in the process of animal oogenesis (Ryu and Macdonald 2015), embryonic development (Keniry et al., 2012), and even gonadal development (Kaneko et al., 2013). The predicted lncRNA from the transcriptome can be helpful for future research on the reproduction mechanism of Muscovy ducks.

Based on the quantitative analysis, we obtain the DEGs in the 3 groups with differential ovulation, and screened the key genes and pathways through KEGG and GO enriched functional annotations. For the KEGG pathway analysis, 5 adhesion-related pathways: Focal adhesion (ko04510), Extracellular matrix (ECM) (ko04512), cell adhesion molecules (CAMs) (ko04514), tight junctions (TJ) (ko04530), and adhesion junctions (AJ) (ko04520) are significant enrichment and make them stand out of the first 20 upregulated pathways (Supplementary Table 4). Focal adhesion is the second-highest enrichment pathway, the ECM receptor pathway acts upstream of the focal adhesion pathway, while the CAMs can achieve intercellular interactions as well as cell-ECM interactions through a variety of mechanisms. Extracellular matrix mainly composed of fibronectin, laminin and collagen (Rodgers and Irving-Rodgers, 2003), and is a highly dynamic structure that will be continuously refactored (Smith et al., 2002), this periodic remodeling is closely related to the development of the ovary, including the formation of follicles, the selection of dominant follicles and the final ovulation (Wang et al., 2019). Previous studies have shown that the outer layer of granular cell (GC) is always above the basement membrane formed by collagen and laminin, and is protected by the follicular basement membrane from the primordial follicle stage to the dominant follicle selection stage of the final ovulation. The presence of laminin in the basement membrane may promote the growth and ovulation of the dominant follicle in follicular selection by enhancing the proliferation of GC (Monniaux et al., 1999; Ożegowska et al., 2019). Cadherins, integrins, CAMs of the immunoglobulin superfamily, and selectins are considered the foremost CAMs classes (Senderoff, 2013), among which N-cadherin plays an indispensable role in follicle formation, and its expression has been increasing during phase of developing follicles. In addition to these 3 pathways, the AJ pathway and TJ pathway have also been significantly enriched. Cadherin, catenin, and α-actin are the main components of AJ. The adhesion junction becomes a strong connection, linking the actin cytoskeleton of adjacent cells (Dowland et al., 2016). AJ and TJ are connected to each other, and AJ regulates TJ components (Suzuki, 2013). TJ protein is involved in follicle and corpus luteum development and also changes with follicle size (Zhang et al., 2018a). In poultry, OCLN participates in the transport of nutrients and regulates the transport of yolk matter through the intercellular signals between granulosa cells. The speed of yolk accumulation promotes the occurrence of follicles to ovulation (Dupont and Scaramuzzi, 2016). These 5 pathways play an important role in ovulation, and the Focal adhesion and ECM-receptor interactions pathways are also significantly upregulated in the IO-CO group. The signal pathways related to cell adhesion undergo significant changes in the follicle selection process. These pathways are closely connected and interact with each other and play an indispensable role in follicle ovulation and ovarian development (Figure 7).

Figure 7

Interaction networking of five adhesion-related KEGG pathways.

In addition, The PI3/Akt, mTOR, FoxO and TGF-β signaling pathway are significant enrichment pathway in the PO-CO (Figure 8). The PI3/Akt signaling pathway plays a key role in the regulation of follicle formation and oogenesis (John et al., 2009). In mammals, IGF-1 stimulates the PI3K/AKT pathway to promote follicle survival, by reducing DNA fragmentation and simultaneously stimulating granulosa cell proliferation. In poultry, IGF-1 stimulates the release of progesterone and affects egg production ( Tosca et al., 2008; Kim et al., 2004). PTEN is the phosphatase of phosphatidylinositol (3,4,5)-trisphosphate (PIP3), which regulates cell proliferation cycles and inhibits cell migration (Leslie and Downes, 2004). The mTOR pathway responds strongly to several growth factors and hormones, so it not only regulates the maturation of oocytes, but also plays a role in different stages of follicular development (Zhang et al., 2017). At the initial stage, maintaining the state of primary follicles depends on the mTOR activity of oocyte-GC communication, while the mTOR dysregulation can lead to abnormal activation of the primary follicle (Zhang et al., 2018b; Makker et al., 2014). At the follicle selection stage, mTOR, as a positive regulator of GC proliferation, regulates the development of follicles (Wen et al., 2015). p70S6K is defined as the downstream signal protein of mTOR complex 1 (mTORC1) ( Liu et al., 2018). The mitogenic signal transmitted to p70S6K is through Akt and mTOR (Laplante and Sabatini, 2013). mTOR complex 2 (mTORC2) promotes actin cytoskeleton rearrangement. Mice with oocyte-specific ablation of the key component of mTORC2 results in abnormal secretion of gonadal hormones, death of a large number of follicles, and even loss of functional follicles. SGK1 is a member of the SGK family. Among them, SGK1 activated by mTORC2 phosphorylates FOXO3a at Ser314, and then FOXO3a inhibits the expression of CDK1 and ultimately promotes cell proliferation (Mori et al., 2014). FOXO3a is one of the most important target genes in the downstream of IGF-1/PTEN/Akt signaling pathway (Huang et al., 2016). FOXO3a can regulate various cellular processes, such as cell cycle, apoptosis, DNA damage repair (Liu et al., 2015). In mammalian ovaries, FOXO3a protein regulates atresia and follicular growth by promoting the apoptosis of ovarian granulosa cells (Matsuda et al., 2011). While in poultry, FOXO3a protein can be used to inhibit the excessive activation of primitive follicles and chicken ovarian granulosa cells ( Cui et al., 2019). These results indicate that FOXO3a plays an indispensable role in animal ovaries. TGF-β signaling pathway is divided into 2 branches through Smad2 and Smad3, TGF-b1 is expressed in oocytes of different stages of follicles, granulosa cells, and follicular membrane cells, and it affects steroid production and regulates follicle development. In recent years, many studies have reported Smad protein transduce the TGF-β signaling pathway, by affecting the proliferation of granulosa cells, regulating the survival of ovarian tissue follicles and maintaining the homeostasis of the follicles (Shen and Wang, 2019). The BMP branch signal through Smads 1/5/8, which play an active role in promoting granulosa cell proliferation and follicle survival (Knight and Glister, 2006). The selection of dominant follicles may depend on the different sensitivity to FSH. Changes in the expression of AMH and BMP in the follicles may regulate FSH-related pathways in granulosa cells, and then participate in the selection of dominant follicles (Chu et al., 2018). Overall, PI3/Akt, mTOR, FoxO, and TGF-β signaling pathway are inextricably related to the development of follicles, and its potential molecular mechanism needs further experiments to clarify

Figure 8

The signal pathway pattern diagram with heat map representation of differentially expressed genes (yellow, higher; blue, lower expression).

Moreover, GO enrichment analysis enhanced our understanding of DEGs; especially, we focused on the GO terms which related to reproduction. There were 78 DEGs involved in reproductive process, 73 DEGs involved in multiorganism reproductive process, and 59 DEGs involved in sexual reproduction. Protein-Protein Interaction network is constructed by the DEGs from reproduction related-GO terms, and we obtain the hub genes (CTNNB1, IGF1, FOXO3, HSPA2, PTEN and SMC4) through the PPI network. In addition to the genes mentioned above, catenin beta 1 (CTNNB1) is a downstream effector molecule of the Wnt pathway and one of the key factors determining sex (Chassot et al., 2014). It promotes the synthesis of estrogen and participates in the meiosis of germ cells and angiogenesis, and promotes the development of preovulatory follicles (Terakawa et al., 2019). Heat shock related protein 2 (HSPA2), although previous studies have mostly focused on its effect on spermatogenesis, recent studies have shown that HSPA2 is involved in the development of primordial follicles to primary follicles by participating in endoplasmic reticulum protein processing (Xu et al., 2017). Structural maintenance of chromosomes 4 (SMC4) participates in the cell cycle and reaches its highest concentration during mitosis, which may have a close relationship with DNA repair (Griese et al., 2010; Wei-Shan et al., 2019).

Note that here, we also found that some DEGs were upregulated significantly only in the IO-UO group. Cysteine and glycine rich protein 1 (CSRP1) is mainly expressed in vascular and visceral smooth muscle cells. It has been proven to play an important role in cattle muscle development (He et al., 2013). WNT works by binding to Frizzled (FZD) receptors. The expression of FZD1 is significantly induced in granulosa cells of follicles before ovulation (Lapointe et al., 2012), while FZD4 is preferentially expressed during ovulation and luteinization, and is highly expressed in granulosa cells of small follicles of postpartum ovaries. FZD4 is not only involved in retinal angiogenesis but also expressed in blood vessels and the surrounding matrix of the embryo (Hsieh et al., 2005). Secreted protein acidic and cysteine rich (SPARC) contributes to oogenesis in ovaries by helping to maintain the nuclear divisions and cytoskeleton of the follicular cells (Irles et al., 2017). Thrombospondin (THBS) family belongs to the group of ECM proteins. Among them, THBS1 and THBS2 have a high degree of structural homology. THBS2 is the form expressed in luteinized granular cells, and has been shown to be involved in regulating angiogenesis (Berisha et al., 2016). Tubulin beta 1 class VI (TUBB1) is an important component of microtubules, it participates in mitosis and cell movement, and plays a key role in maintaining cell stability (Zhao et al., 2020). These genes may be the key genes affecting continuous ovulation in Muscovy duck; we need to explain through further experiments.

CONCLUSIONS

Taken together, this study obtained the full-length transcript by using ONT, by the way of qualitative analysis; we obtain the structural outcomes and gene functional annotation of novel transcript. More importantly, in accordance with GO and KEGG enrichment analysis, a large number of ovulation related genes and pathways were obtained. We hypothesize that CTNNB1, IGF1, FOXO3, HSPA2, PTEN, and SMC4 genes may be the key genes affecting reproduction and play an important role in ovulation in Muscovy ducks. In addition, we preliminarily explored the potential relationship and role of Adhesion-related pathway, mTOR pathway, TGF-β signaling pathway and FoxO signaling pathways in the follicular development of Muscovy duck. These genes and pathways may play a decisive role in the follicular development, selection and ovulation of Muscovy duck. This research will help to gain a deeper understanding of the selection of follicles in Muscovy duck and even poultry and the regulatory mechanism of ovarian development.

ACKNOWLEDGMENTS

This study was supported by funds from the National Waterfowl Industrial Technology System (CARS-42) and the Research on the Effect of New Breeding of the Black Muscovy Duck (2019QYPY-130).

DISCLOSURES

The authors declared that they have no conflicts of interest to this work. The authors declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

Table 1

Information of clean reads mapped with the reference transcriptome.

Sample	Total reads	Mapped reads	Mapped rates %
CO1	4,842,304	3,940,400	81.37%
CO2	3,630,592	2,942,843	81.06%
CO3	5,254,548	4,166,762	79.30%
IO1	6,477,736	5,693,539	87.89%
IO2	4,237,737	3,459,506	81.64%
IO3	4,842,807	3,903,137	80.60%
PO1	4,743,008	3,843,865	81.04%
PO2	5,432,493	4,462,321	82.14%
PO3	4,599,417	3,761,808	81.79%

Abbreviations: CO, consecutive ovulation; IO, inconsecutive ovulation; PO, preovulation.