Comparative Transcriptome Profiling of Ovary Tissue between Black Muscovy Duck and White Muscovy Duck with High- and Low-Egg Production

The egg-laying rate is an important indicator for evaluating fertility of poultry. In order to better understand the laying mechanism of Muscovy ducks, gene expression profiles and pathways of ovarian tissues in high- and low-laying black (BH and BL) and white Muscovy ducks (WH and WL) during the peak production period were performed by using RNA-seq. The total number of reads produced for each ovarian sample ranged from 44,344,070 to 47,963,328. A total of 113, 619 and 87 differentially expressed genes (DEGs) were identified in BH-vs-WH, BL-vs-BH and BL-vs-WL, respectively. Among them, 54, 356 and 49 genes were up regulated and 59, 263 and 38 genes were down regulated. In addition, there were only 10 up-regulated genes in WL-vs-WH. In the comparison of DEGs in black and white Muscovy ducks, two co-expressed DEG genes were detected between BH-vs-WH and BL-vs-WL and seven DEGs were co-expressed between BL-vs-BH and WL-vs-WH. The RNA-Seq data were confirmed to be reliable by qPCR. Numerous DEGs known to be involved in ovarian development were identified, including TGFβ2, NGFR, CEBPD, CPEB2, POSTN, SMOC1, FGF18, EFNA5 and SDC4. Gene Ontology (GO) annotations indicated that DEGs related to ovarian development were mainly enriched in biological processes of “circadian sleep/wake cycle process,” “negative regulation of transforming growth factor-β secretion,” “positive regulation of calcium ion transport” in BH-vs-WH and “cell surface receptor signaling pathway,” “Notch signaling pathway” and “calcium ion transport” in BL-vs-BH. Besides, “steroid biosynthetic process,” “granulosa cell development” and “egg coat formation” were mainly enriched in BL-vs-WL and “reproduction,” “MAPK cascade” and “mitotic cell cycle” were mainly enriched in WL-vs-WH. KEGG pathway analysis showed that the PI3K-Akt signaling pathway and ovarian steroidogenesis were the most enriched in Muscovy duck ovary transcriptome data. This work highlights potential genes and pathways that may affect ovarian development in Muscovy duck.


Introduction
Reproductive performance, an important economic trait of poultry, is mainly determined by the egg production. Improving egg production ability is an important breeding objective. With traditional breeding methods, the reproductive performance of laying duck has been gradually improved but more significant improvements have been slow [1]. Ovary is a female reproductive organ, which is closely related to egg-laying characteristics and the ovarian function of poultry directly affects the egg production [2]. Ovarian follicular development is a complex process, involving a variety of endocrine, autocrine and paracrine factors, which control the proliferation and differentiation of oocytes, granulosa cells and thecal cells [3]. The number of primordial follicle pools in the ovary is fixed and the quiescence, survival and activation of follicles in the pools depends on the dynamic balance, which is regulated by many signaling molecules or pathways [4,5]. Genetic or environmental factors affect follicular development, which may lead to follicular inactivation or premature activation and even degeneration and atresia at different stages of development [6]. Therefore, the regulation of ovarian function, including the identification and analysis of reproductive-related genes, has been concerned to improve the reproductive performance of poultry. For example, transforming growth factor-β family (TGF-β), fibroblast growth factors (FGFs), insulin-like growth factor family (IGFs) [7], FSHR [8], GnRH [9], STAR [10] and PRLR [11] have been shown to be associated with poultry ovarian function and affect egg production.
Although candidate genes for avian reproduction have been identified, the specific molecular mechanisms and related signaling pathways are poorly understood. As a nextgeneration sequencing technology, RNA-seq is used to explore the function of genes by detecting the gene expression in animal samples [12]. In recent years, the ovarian transcriptome of livestock, such as pig [13], cattle [14], sheep [15], has been studied by using RNA-seq. Similarly, there are many studies on the transcriptome of poultry ovary. Zhang et al. (2019) found that there were 305 differentially expressed genes (DEGs) identified in the ovaries of hens with relatively greater and lesser egg production. Besides, five candidate genes related to egg production, including ZP2, WNT4, AMH, IGF1 and CYP17A1 genes, were discovered [16]. Tao et al. (2017) performed a comparative transcriptomic analysis of the ovaries of Jinding ducks with high and low egg production using RNA-seq, the results showed that 367 of 843 DEGs were down regulated in high egg production ovaries and 476 DEGs were up-regulated in low egg production ovaries [17]. Luan et al. (2014) identified the DEGs in the ovaries of Huoyan geese between the laying period and the ceased period, 344 up-and 344 down-regulated genes were classified into functional categories and 12 DEGs that were mainly involved in pathways for reproduction regulation [18]. Xu et al. (2016) compared the gene expression of pigeon ovary at different stages, there were 409 DEGs between pre-and post-ovulation. However, there are few reports on the high-and low-yield ovarian transcriptome of Muscovy duck [19].
Muscovy duck, known for its strong adaptability, fertility and high meat yield [20], is native to tropical areas of Central and South America. In China, black Muscovy and white Muscovy are mainly raised. The egg-laying peak period of Muscovy duck is between 35 and 53 weeks, covering most of the egg-laying period [21]. Although these ducks have been raised on a large scale in China, the poor egg production performance and the lack of systematic breeding have resulted in uneven production performance in Muscovy duck and seriously affects the economic benefits of many enterprises. Under the same conditions of nutrition and feeding management, the egg production performance of Muscovy duck is mainly determined by genes and related regulatory factors. Therefore, the purpose of our research was to identify the DEGs and pathways and analyze the similarities and differences of the differential expression in the ovarian tissues of Muscovy ducks with highand low-yield by using RNA-seq. Through a comprehensive analysis of gene expression levels reflecting the ovaries of black and white Muscovy ducks, the genes and molecular mechanisms involved in this process have been found. All results provide theoretical basis for the breeding of Muscovy ducks and are useful for understanding of the molecular regulation mechanism of egg production characteristics in black and white Muscovy duck.

Animal and Sample Collection
The egg production of Muscovy ducks aged 216 to 280 days (65 days in total) was recorded. More than 60 eggs were considered as high-laying black Muscovy duck (BH) and white Muscovy duck (WH) and 15-20 eggs were defined as low-yield duck (BL and WL). Six black Muscovy ducks (3 high-yield and 3 low-yield) and 6 white Muscovy ducks (half each for high-and low-yield laying ducks), raised under the same environmental conditions with free access to feed and water (Table 1), were randomly selected and slaughtered quickly to collect the ovarian tissues, which comprised the whole ovary including the small and large yellow follicles. The samples were immediately frozen in liquid nitrogen and stored at −80 • C until further use. Total RNA was extracted from ovarian tissue of Muscovy duck using TRIzol reagent according to the manufacturer's instructions (Takara, Dalian, China). The purity and concentration of RNA were detected by agarose gel electrophoresis and Nanodrop 2000 (Thermo, Waltham, MA, USA) and RNA integrity was measured using the 2100 Bioanalyzer (Agilent Technologies, San Jose, CA, USA). All samples had an RNA integrity number (RIN) ≥ 7.8 and a RNA concentration ≥ 875.6 ng/µL (Table S1). The tested Muscovy ducks were purchased from Anda Farm in Tongguan County, Shaanxi Province, China and raised according to the standard procedure. Besides this, the Muscovy ducks were known to be unrelated to each other for at least three generations. All animal procedures were performed according to guidelines provided by the China Council on Animal Care and the protocols were approved by the Experimental Animal Management Committee (EAMC) of Northwest A&F University (ethic code: #0726/2018).

Library Construction and Sequencing
Twelve RNA samples with high quality (average RIN: 8.4 ± 0.251) and concentration (average concentration: 1136.9 ± 3.828 ng/µL) were used to construct the sequencing libraries and high-throughput sequencing was performed on the Illumina novaseq 6000 platform (San Diego, CA, USA) following the manufacturer's recommendations, generating 150 bp paired-end reads. The high-quality clean reads were obtained by filtering the raw reads and removing: (1) the sequences containing adapters; (2) the sequences with more than 10% of N bases; (3) the sequences with more than 50% base quality values less than 10. The clean reads were aligned to the Anas platyrhynchos genome (https:// www.ncbi.nlm.nih.gov/genome/?term=duck) and annotated transcripts (https://www. ncbi.nlm.nih.gov/assembly/GCF_003850225.1) by way of the Tophat 2 software (http: //ccb.jhu.edu/software/hisat2/index.shtml).

Differential Expression Analysis
The number of clean reads for each gene was calculated and the Fragment Per Kilobases per Million Fragment (FPKM) was used to estimate the expression abundance of transcripts from different samples [22]. Differential expression analyses of the 4 comparison groups (BL-vs-BH, WL-vs-WH, BH-vs-WH and BL-vs-WL) were performed using the DESeq R package (http://www.bioconductor.org/packages/release/bioc/html/DESeq. html) [23,24] and genes with an p-value ≤ 0.05 and an expression |log2 Fold | ≥ 1were identified as DEGs.

GO and KEGG Analysis
To annotate the function of these DEGs, Gene Ontology (GO) analysis was conducted by using the GOseq software (https://bioconductor.org/packages/release/bioc/html/ goseq.html) for each of the three main categories: biological process, cellular component and molecular function. Biological pathways enriched for the identified DEGs through KEGG pathway analyses were carried out using the KOBAS software (http://kobas.cbi. pku.edu.cn/kobas3/?t=1).

Verification of RNA-seq by qPCR
To verify the reliability of the RNA-Seq results, 4 DEGs (TGFβ2, FGF18, POSTN, SMOC1) in the ovary of black Muscovy duck and 3 DEGs (POSTN, SMOC1, COL4A1) in white Muscovy duck were randomly selected, respectively. The β-actin gene of duck was used as the control (NC_040060.1), after the stability of β-actin was tested in this system. The RNA extracted from same ovaries with an RIN ≥ 7.8, A260/280 >1.8 and A260/230 > 2.0 were used for cDNA synthesis according to the manufacturer's instructions (Takara, Dalian, China). The primers information were listed in Table 2. qPCR was run in triplicate and performed in a EcoRT48 system (OSA, London, UK), with a 10 µL reaction volume containing 5 µL SYBR Premix Ex Taq II (Takara, Dalian, China), 0.8 µL cDNA (750 ng/µL), 0.2 µL of each primer and 3.8 µL RNase-free ddH 2 O. The qPCR program was 95 • C for 30 s; followed by 40 cycles of 95 • C for 5 s, 60 • C for 30 s and then 95 • C for 15 s, 55 • C for 15 s, 95 • C for 15 s. Gene expression results were calculated using the 2 −∆∆Ct method.

Statistical Analysis
Experiments were repeated 3 times and presented as mean ± standard deviation and conducted using one-way ANOVA with Dunnet's t-test at p < 0.05 probability levels in SPSS18.0. "*" was considered significant difference (p < 0.05); "**" was considered extremely significant difference (p < 0.01).

High-Throughput Sequencing and Read Mapping
In this study, a total of 12 libraries in high-and low-yield of ovarian tissue in Muscovy ducks were established by high-throughput RNA sequencing. The clean reads of each sample were more than 44 million, the mapped reads and unique mapped reads ranged from 60.28% to 67.40% and 57.81% to 64.36%, respectively. In addition, more than 16,276 genes were detected in each sample (Table 3)

Differentially Expressed Genes
According to the correlation analysis of 12 samples, the differences between biological replicates were small and the repeatability was high, which indicated that the selection of experimental samples was consistent and reliable ( Figure S1). A comparison of the gene expression showed that a total of 113 genes were differentially expressed between BH and WH, including 54 up-regulated genes and 59 down-regulated genes. 619 DEGs were discovered in BL-vs-BH and 356 DEGs were up-regulated genes and 263 DEGs were down-regulated genes. In BL-vs-WL, there were 87 DEGs, of which 49 were up-regulated and 38 were down-regulated genes. In addition, compared with WH, there were only 10 up-regulated genes in WL group (Figures 1 and 2).
A comparative analysis of DEGs obtained by black and white Muscovy ducks showed that two co-expressed DEG genes were found between BH-vs-WH and BL-vs-WL and 7 co-expressed DEGs between BL-vs-BH and WL-vs-WH ( Figure 3).

Gene Ontology and KEGG Analysis
To further assess the functional roles of DEGs in ovarian development, a GO

Gene Ontology and KEGG Analysis
To further assess the functional roles of DEGs in ovarian development, a GO enrichment analysis was performed. The most significant biological processes related to ovarian development in BH-vs-WH, included "circadian sleep/wake cycle process," "negative regulation of transforming growth factor-β secretion" and "positive regulation of calcium ion transport." The most relevant terms for cellular component were "actin cytoskeleton" and "focal adhesion." Significant GO terms for molecular function included "extracellular matrix structural constituent" and "calcium channel regulator activity." In BL-vs-BH, DEGs related to ovarian development were mainly enriched in biological processes of "cell surface receptor signaling pathway," "Notch signaling pathway" and "calcium ion transport." "Voltage-gated calcium channel complex" and "proteasome core complex" were mainly enriched in cellular component and "calcium ion binding" and "low-density lipoprotein receptor activity" were mainly enriched in molecular function. The common DEGs related to ovarian development in BL-vs-WL were primarily enriched in biological processes of "steroid biosynthetic process," "granulosa cell development" and "egg coat formation." DEGs were primarily enriched in cellular component, including "voltage-gated calcium channel complex" and "cyclin-dependent protein kinase holoenzyme complex." There were some terms mainly enriched in molecular function of "transforming growth factor β receptor binding" and "calcium ion binding." Besides, several fundamental biological processes related to ovarian development were found to be notably enriched in WL-vs-WH, such as "reproduction," "MAPK cascade" and "mitotic cell cycle" and several cellular component were found to be notably enriched, such as "cyclin-dependent protein kinase holoenzyme complex" and "proteasome complex." Furthermore, DEGs were primarily enriched in molecular function, including "translation regulator activity" and "lipid kinase activity" (Figure 4).  To better understand the biological functions and interaction of genes, a KEGG pathway analysis was performed for the DEGs identified. In the first 20 pathways enriched by DEGs of the comparisons, the KEGG pathways enriched significantly (P-value < 0.05) were mainly concentrated in PI3K-Akt signaling pathway, Rap1 signaling pathway, cell adhesion molecules (CAMs), thyroid hormone signaling pathway, focal adhesion, ECM-receptor interaction, ovarian steroidogenesis, progesterone-mediated oocyte maturation, oocyte meiosis and TGF-β signaling pathway. Among them, PI3K-Akt signaling pathway and ovarian steroidogenesis were involved in the reproductive process of black and white Muscovy ducks ( Figure 5). To better understand the biological functions and interaction of genes, a KEGG pathway analysis was performed for the DEGs identified. In the first 20 pathways enriched by DEGs of the comparisons, the KEGG pathways enriched significantly (Pvalue < 0.05) were mainly concentrated in PI3K-Akt signaling pathway, Rap1 signaling pathway, cell adhesion molecules (CAMs), thyroid hormone signaling pathway, focal adhesion, ECM-receptor interaction, ovarian steroidogenesis, progesterone-mediated oocyte maturation, oocyte meiosis and TGF-β signaling pathway. Among them, PI3K-Akt signaling pathway and ovarian steroidogenesis were involved in the reproductive process of black and white Muscovy ducks ( Figure 5). In order to verify the expression levels of DEGs observed in our RNA-seq analysis, 4 randomly selected transcripts in BL-vs-BH (TGFβ2, FGF18, POSTN and SMOC1) and 3 randomly selected transcripts WL-vs-WH (POSTN, SMOC1 and COL4A1) shown to be differentially expressed on the basis of RPKM values were validated by qPCR. The qPCR results showed the similar regulated trend in the expression of these genes. Therefore, these results indicated the high reliability and accuracy of the RNA-seq data ( Figure 6). In order to verify the expression levels of DEGs observed in our RNA-seq analysis, 4 randomly selected transcripts in BL-vs-BH (TGFβ2, FGF18, POSTN and SMOC1) and 3 randomly selected transcripts WL-vs-WH (POSTN, SMOC1 and COL4A1) shown to be differentially expressed on the basis of RPKM values were validated by qPCR. The qPCR results showed the similar regulated trend in the expression of these genes. Therefore, these results indicated the high reliability and accuracy of the RNA-seq data ( Figure 6).

Discussion
Ovarian development and follicle formation in poultry is a dynamic and complex process, including follicle growth and selection, granulosa cell proliferation and differentiation and so forth. It is very important to understand the gene expression patterns in the ovary to improve the egg production rate in poultry [25]. In this study, the gene expression mechanisms and signal pathways in the ovaries of high-and low-yield Muscovy ducks were investigated by using RNA-seq technology. A total of 555,951,128 clean reads (46,329,260 on average) were discovered and the Mapped reads and Unique mapped reads were more than 60.28% and 57.81%, respectively. Besides, more than 16,276 genes were detected in each sample.

Analysis of DEGs
According to the DEGs in black Muscovy duck with high-and low-laying, several genes that may affect Muscovy duck egg production were identified, such as TGFβ2, NGFR, CEBPD, CPEB2, POSTN, SMOC1, FGF18, EFNA5 and SDC4. Moreover, some studies have shown that these regulatory transcription factors interacted with each other in regulating animal reproduction.
Transforming growth factor-β (TGF-β), one of the large family of multifunctional growth factors, regulates a wide range of biological activities related to morphogenesis, development and differentiation. There is a wide range of biological effects in TGF-β isoforms, including cell differentiation, cell proliferation, cell growth, extracellular matrix formation and immune function [26]. Some studies have found that the TGFβ2 gene may be involved in reproduction, such as Hu sheep [27]. Before maturation of porcine oocytes, the expression of TGFβ2 in large follicles was greater than that in small and medium follicles [28]. There is evidence that the up regulation of TGFβ2 may induce apoptosis

Discussion
Ovarian development and follicle formation in poultry is a dynamic and complex process, including follicle growth and selection, granulosa cell proliferation and differentiation and so forth. It is very important to understand the gene expression patterns in the ovary to improve the egg production rate in poultry [25]. In this study, the gene expression mechanisms and signal pathways in the ovaries of high-and low-yield Muscovy ducks were investigated by using RNA-seq technology. A total of 555,951,128 clean reads (46,329,260 on average) were discovered and the Mapped reads and Unique mapped reads were more than 60.28% and 57.81%, respectively. Besides, more than 16,276 genes were detected in each sample.

Analysis of DEGs
According to the DEGs in black Muscovy duck with high-and low-laying, several genes that may affect Muscovy duck egg production were identified, such as TGFβ2, NGFR, CEBPD, CPEB2, POSTN, SMOC1, FGF18, EFNA5 and SDC4. Moreover, some studies have shown that these regulatory transcription factors interacted with each other in regulating animal reproduction.
Transforming growth factor-β (TGF-β), one of the large family of multifunctional growth factors, regulates a wide range of biological activities related to morphogenesis, development and differentiation. There is a wide range of biological effects in TGF-β isoforms, including cell differentiation, cell proliferation, cell growth, extracellular matrix formation and immune function [26]. Some studies have found that the TGFβ2 gene may be involved in reproduction, such as Hu sheep [27]. Before maturation of porcine oocytes, the expression of TGFβ2 in large follicles was greater than that in small and medium follicles [28]. There is evidence that the up regulation of TGFβ2 may induce apoptosis during the regression of postovulatory follicle in chickens [29]. Ma et al. (2002) found that TGFβ2 was preferentially expressed in the ovaries and testes of chicken embryos, which may play a key role in regulating the development of germ cells of ovaries and testes in chicken embryos [30]. Nerve growth factor receptor (NGFR) is a multifunctional cell surface receptor, which widely exists in many cell types. The gene induces apoptosis and participates in injury, nervous system development and regeneration [31]. NGFR is expressed in female reproductive organs and plays an important role in regulating the growth, development and function of reproductive organs [32]. Studies have shown that NGFR was expressed in the follicles with different sizes in porcine ovaries during estrus, which may have an important impact on the ovarian function of sows [33]. In addition, NGFR is also expressed in the granulosa cells, thecal cells, stromal cells and luteal cells of goat ovarian [34]. The transcription factor C/EBPδ (CEBPD) is a member of the basic leucine zipper family transcription factor, which is related to cell proliferation, differentiation, growth arrest and apoptosis [35]. Becker et al. (2011) found that CEBPD had a higher proliferation rate in calf endometrium under the influence of steroid hormones [36]. All these studies have identified that the DEGs have important functions in the ovary. In this work, TGFβ2, NGFR and CEBPD were significantly expressed in high-laying Muscovy ducks, indicating that they may be very important for ovarian development in Muscovy ducks.
The cytoplasmic polyadenylation element binding protein family (CPEB), a sequencespecific RNA binding protein, regulates the physiological processes of germ cell development, cell division and cell differentiation by regulating the translation process [37,38]. The CPEB2 gene belongs to the CEPB family and has been discovered in mouse germ cells [39]. Its function is to promote meiotic maturation and embryonic development of porcine oocytes [40]. Periostin (POSTN) is a 90-kDa adhesion molecule that is secreted in the extracellular matrix. POSTN binds to cells via the Integrin/FAK/AKT signaling pathway and promotes cell migration, adhesion and survival of various cell types [41,42]. POSTN is associated with follicular granulosa cell proliferation [43,44]. POSTN affects the protein in the follicle and the insufficient protein in the follicular fluid may lead to unique characteristics of follicle development in lambs [45]. SPARC-related modular calcium binding 1 (SMOC1), one of the extracellular glycoproteins, is widely expressed in various tissues with localization to basement membranes. SMOC1 may mediate intercellular signaling and cell type-specific differentiation during gonadal and reproductive tract development of fetal [46]. In this study, the differential expression of CEPB2, POSTN and SMOC1, which was in the ovarian tissues of high-and low-laying black and white Muscovy ducks, may have an important regulatory effect on the egg-laying mechanism of Muscovy ducks.
FGF18, a member of the FGF family, participates in the development of ovarian follicles [47]. FGF18 promotes the mitosis of thecal cells and apoptosis of bovine granulosa cells [48,49] and plays an important role in the regulation of steroid production in fetal bovine ovary [50]. Besides, Zhong et al. (2006) found that FGF18 gene was also presented in mouse thecal cells and oocytes [51]. Ephrin-A5 (EFNA5) is a member of the Ephrin-Eph family, which is highly expressed in granule cells [52]. EFNA5 has pivotal roles in cellular processes such as proliferation, survival, cell migration and invasion. Besides, EFNA5 also mediates cell apoptosis, proliferation and steroid production during follicular development in female mice [53]. In our study, FGF18 and EFNA5 were mainly expressed in the ovarian tissue of high-yield black Muscovy ducks, indicating that FGF18 and EFNA5 may play a key role in ovarian granulosa cell apoptosis and follicular development.
Syndecan4 (SDC4) is a cell surface heparan sulfate proteoglycan, which participates in cell adhesion [54,55] and is widely expressed in most adult tissues [56]. Studies have shown that the expression of SDC4 is related to follicular atresia in mice [57]. In addition, compared with the dominant follicles of cattle, the upregulation of SDC4 gene in ovulatory follicles may affect the ovulation promoting or luteinizing process [58]. In this study, the expression of SDC4 in the ovarian tissue of high-laying white Muscovy ducks was significantly upregulated, which may be related to the high-egg production of Muscovy ducks.

Analysis of GO and KEGG
To better understand the possible functions involved in ovarian development in Muscovy ducks, the DEGs from each comparison group were annotated with the GO database. the results indicated that the DEGs played an important role in regulating duck ovarian development. Functional classification and pathway assignments were based on a KEGG analysis. In this study, 2 major pathways related to ovarian development were identified, namely PI3K-Akt signaling pathway and ovarian steroidogenesis. The PI3K/Akt pathway is an intracellular signaling pathway of great importance in the cell cycle process. Its plays an important role in regulating metabolism and has been widely investigated in cell proliferation, cell transformation, paracrine function and angiogenesis [59]. In the PI3K/Akt signal transduction pathway, with the help of phosphatidylinositol dependent kinase, protein kinase B (Akt) binds to cell membrane and exert biological effects [60]. Phosphatidylinositol-4,5-diphosphate 3-kinase (PI3K) is activated by many genes and it has been shown to be closely related to follicular recruitment and ovulation, selection and development of dominant follicles in mammalian ovaries [61][62][63][64].
Ovarian cell proliferation and apoptosis are important factors affecting ovarian function. Recent studies have shown that the follicular development and granulosa cell survival and proliferation can be improved by regulating the expression of PI3K/Akt signaling pathway [65]. Mishra et al. (2020) found that DEGs associated with high rates of egg production regulate PI3K-Akt signaling pathway at hypothalamic-pituitary-ovarian axis in laying hens [66]. The ovaries with reproductive activity are composed of small prehierarchical and maturing preovulatory follicles showing a hierarchy according to their size [67]. In addition, ovarian function is regulated by gonadotropins, including follicle stimulating hormone (FSH), luteinizing hormone (LH) and ovarian steroids. In animals, steroid hormones are synthesized from cholesterol by LH responsive theca cells and then converted into estrogens [68]. With the maturation of follicles, the concentration of steroid hormones such as estrogen and progesterone increases to exert their reproductive functions. The pathways in this study strongly indicated that PI3K-Akt signaling pathway and ovarian steroidogenesis played an important role in follicle development and productivity of Muscovy ducks.
Finally, we validated the RNA-seq results using qPCR to measure the expression of five DEGs (TGFβ2, FGF18, POSTN, SMOC1 and COL4A1), which showed that the results were reliable. In these five DEGs, the expression of POSTN and SMOC1 were high significantly different in BL-vs-BH and WL-vs-WH, indicating that POSTN and SMOC1 played an important role in regulating ovarian development of Muscovy ducks. Besides, the expression of TGFβ2 and FGF18 in BH were significantly or extremely significantly higher than that in BL and the expression of COL4A1 in WH were extremely significantly higher than that in WL, respectively. These results suggested that TGFβ2, FGF18 and COL4A1 may have relatively more important function in the regulation of egg production in Muscovy ducks.

Conclusions
In this study, several genes that may be important candidate genes involved in Muscovy duck ovarian development, such as TGFβ2, NGFR, CEBPD, CPEB2, POSTN, SMOC1, FGF18, EFNA5 and SDC4, were discussed. Go term annotation and KEGG pathway enrichment analysis showed that a number of significantly enriched DEGs may be related to PI3K-Akt signaling pathway and ovarian steroidogenesis in ovaries of Muscovy ducks. These data provide comprehensive gene expression information at the transcriptional level and will certainly accelerate the study of Muscovy duck ovarian development.

Data Availability Statement:
The datasets generated for this study can be found in the NCBI SRA (Submission: SUB7071226). Bioproject #PRJNA641492 and Biosamples #SAMN14252184-SAMN14252193.