Shan Ma ducks were provided by the Shan Ma Duck Stock Breeding Farm (Longyan, Fujian, China), and kept under the same environmental conditions with free access to feed, water, and commercial corn-soybean-based diets. From hatching to 4 weeks of age, the birds received a starter feed (metabolizable energy [ME], 12.08 MJ/kg; crude protein [CP], 19.0%). From 5 weeks of age to slaughter the birds were fed a grower diet (ME, 11.5 MJ/kg; CP, 16.5%). Six adult female Shan Ma ducks were included in these studies. Samples of ovarian tissue were collected from peak stage of egg production (180 days old) and late stage of egg production (520 days old) (n = 3 each), immediately frozen in liquid nitrogen, and stored at −80°C for later use. Generally, the first egg age of Shan Ma duck was 108 days old, and the peak egg production rate of Shan Ma duck was reached over 95% at 180 days old, and the egg production rate dropped to about 40% to 50% at 520 days old.

Total RNA from ovarian tissue of the 6 ducks was extracted using Trizol reagent and dissolved in diethyl pyrocarbonate water. The concentrations of the RNA samples were determined using a Nanodrop-2000 UV-Vis Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) for nucleic acids and proteins. Each of the peak egg-production with 10 μg RNAs were pooled and each of the late egg-production with 10 μg RNAs were pooled. We froze the 2 pooled RNA samples at −80°C for later use.

We used magnetic beads coupled with poly dT to isolate mRNAs from total RNAs in order to obtain mRNA-enriched RNA samples. mRNA-enriched RNAs were then fragmented by sonication. RNA fragments of 200 to 500 bp were reverse-transcribed using random primers to generate first-strand cDNAs. The first-strand cDNAs were used as templates to synthesize the complementary second strand in a reaction system containing buffer, dNTPs, RNase H, and DNA polymerase I. The double-stranded cDNAs were end-repaired, dA-tailed, and ligated to an adaptor for polymerase chain reaction (PCR) amplification. The cDNA library was then ready for sequencing.

cDNAs were sequenced by Beijing Biomarker Technologies Co., Ltd, using the Illumina HiSeq 2000 system. The original sequence reads were aligned to the reference duck genome (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF_000355885.1_BGI_duck_1.0/GCF_000355885.1_BGI_duck_1.0_genomic.fna.gz) using TopHat2 software [11], allowing up to two base mismatches. The gene expression level was then calculated using the fragments per kilobase of exon per million fragments (FPKM) method by CUFFLINKS V2.1.1 [12]. The expression levels of a gene in different samples were subject to Chi-square test using IDEG6 software (http://telethon.bio.unipd.it/bioinfo/IDEG6). The p values obtained were corrected with the false discovery ratio. The genes with corrected p values less than 0.01 and a ratio of the FPKM of the 2 samples equal or greater than 2 (fold-change≥2) were considered to be significantly differentially expressed.

The differentially expressed genes were analyzed in the NR database, and the annotation of the homologous gene in the NR database was used to annotate the differentially expressed genes identified. BLASTX was used to align differentially expressed genes with the orthologous genes from the clusters of orthologous groups of proteins (COG) database to obtain the COG annotation and classification of the genes. BLAST2GO [13] was used to align the differentially expressed genes with genes in the gene ontology (GO) database to obtain the GO annotation of the genes. The obtained annotation was enriched and refined using TopGo (R package version 1.16.2) with the “elim” method and Kolmogorov–Smirnov test. The differentially expressed genes were also aligned in BLASTX with genes in the Kyoto encyclopedia of genes and genomes (KEGG) database to obtain pathway annotation for the genes. KEGG pathways were enriched using the right sided Fisher’s exact test and p value ≤0.05 was taken as a threshold of significance.

To validate the reliability of Illumina analysis, five differentially expressed genes in the significantly enriched pathways were selected and detected by real-time quantitative PCR (qRT-PCR). The relative expression level of each gene was normalized to the control gene β-actin (NM_001310421.1) using the 2^−ΔΔCt method [14]. The qRT-PCR primers were designed using Primer Premier 6.0 and Beacon designer 7.8. The primer pairs used for qRT-PCR are listed in Table 1. All primers were synthesized by Invitrogen Biotechnology Co., Ltd (Carlsbad, CA, USA). The qRT-PCR analysis was conducted on the Bio-Rad CFX384 real-time PCR detection System (BIO-RAD, City, CA, USA). The qRT-PCR reaction contained 1 μL of cDNA template, 10.0 μL of Power SYBR Green Master Mix (Applied Biosystem, Foster, CA, USA), 8.0 μL of sterile water, and 0.5 μL of each gene-specific primer. Thermal cycling conditions were 1 cycle at 95°C for 2 min, 40 cycles of 95°C for 15 s and 63°C for 25 s. There were three replicates for each amplification.

The quality of RNA preparations used for transcriptome sequencing was examined (Table 2). The optical density 260/optical density 280 (OD₂₆₀/OD₂₈₀) ratios were between 1.8 and 2.2, the RNA integrity number values were greater than 8.4, and the 28S/18S ratios were greater than 1. All of these results indicated that the quality of the RNAs used in this study was satisfactory for transcriptome sequencing.

The quality control results of the ovarian transcriptomic sequencing data from ducks at their peak or late stages of egg production are summarized in Table 3, 4. Table 3 shows that the two samples obtained 42.78 M (million) and 45.32 M reads, respectively, and more than 4.3 Gb of base pairs. The Q20 value was greater than 93% and the Q30 value was greater than 84%. These results indicated that the quality of the transcriptome sequencing data was reliable and could be used for subsequent analysis.

Table 4 shows that the total numbers of reads from sequencing were 42,782,676 and 45,316,166, respectively, for the two samples. Among the total number of reads 31,101,460 and 32,324,534 reads were mapped when compared to the reference genome, accounting for 72.70% and 71.33% of the total reads, respectively. The percentage of mapped reads was high for both samples.

The FPKM values of the ovarian tissue from ducks at the peak and late stages of egg production were calculated, as shown in Table 5. At the peak stage, 52.00% of the genes had FPKM values less than 5, and 8.10% of the genes had FPKM values greater than 100. At the late stage, 53.34% of the genes had FPKM values less than 5, and 7.95% of the genes had FPKM values greater than 100. In order to identify the differentially expressed genes, we used IDEG6 software. The genes were evaluated by both fold-change and significance level (p value). Only genes with a fold-change equal to or greater than 2 and a significance level of p< 0.05 were considered differentially expressed. As shown in the volcano plot (Figure 1), there were many differentially expressed genes between ducks at the peak and late stages of their egg-producing performance. We identified 2,002 differentially expressed genes, including 1,212 down-regulated genes and 790 up-regulated genes (Figure 2). These genes were further analyzed in the NR database, and we noted that 1,645 genes were already annotated, among which 999 genes were down-regulated and 646 genes were up-regulated. The top ten up- and down-regulated DEGs from the total of 1,645 DEGs identified between peak and late stage were: homeobox protein Hox-A3 (HOXA3), dihydrofolate reductase (DHFR), cerebellar degeneration related protein 2 like (CDR2L), forkhead box O1 (FOXO1), ER membrane protein complex subunit 3 (EMC3), G protein-coupled receptor 20 (GPR20), stem-loop binding protein (SLBP), acyl-CoA synthetase long-chain family member 3 (ACSL3), cytochrome P450 family 7 subfamily A member 1 (CYP7A1), ATP synthase, alpha subunit 1 (ATP5A1), sphingosine-1-phosphate receptor 3 (S1PR3), and TNF receptor associated protein 1 (TRAP1), transmembrane protein 50B (TMEM50B), nebulin-related anchoring protein (NRAP), glycine receptor, beta (GLRB), zinc finger, CCHC domain containing 24 (ZCCHC24), ribosomal protein L37a (RPL37A), secretion regulating guanine nucleotide exchange factor (SERGEF), solute carrier family 29, member 3 (SLC29A3), THAP domain containing 11 (THAP11), Phosphatidylinositol glycan anchor biosynthesis, class V (PIGV), respectively (Table 6).

The differentially expressed genes annotated in NR database were analyzed in the COG database for functional annotation. The results showed that 696 genes received COG function annotations in 25 functional categories (Figure 3). Among the 25 categories, the one with the greatest number of genes was R (general function prediction only), which had 171 genes; followed by categories K (transcription), L (replication, recombination and repair), and T (signal transduction mechanisms), with 60, 60, and 55 genes each, respectively.

The differentially expressed genes annotated in NR database were analyzed using BLAST2GO for functional annotation and genomic dataset analysis in the GO database. GO analysis showed that 1,204 genes received GO functional annotation in 3 domains: cellular component, molecular function, and biological process, and 56 categories (Figure 4). There were 17 categories in the cellular component domain, 17 categories in the molecular function domain, and 22 categories in the biological process domain. In the cellular component domain, the category with the highest number of reads was cell part, which had 804 genes, accounting for the highest percentage—66.78%; the next categories were cell and organelle, with 791 and 610 genes each, accounting for 65.70% and 50.66%, respectively; and the category with the least number of genes was virion, which had only 1 read. In the molecular function domain, the category with the highest number of genes was binding, which had 681 genes, accounting for the highest percentage—56.56%; the next category was catalytic activity, with 430 genes, accounting for 35.71%; and the categories with the fewest genes were translation regulatory activity, chemorepellent activity, nutrient reservoir activity and chemoattractant activity, with only 1 gene each. In the biological process domain, the category with the highest number of genes was cellular process, which had 794 genes, accounting for the highest percentage—65.95%; the next categories were biological regulation and metabolic process, with 587 and 585 genes each, accounting for 48.75% and 48.59%, respectively; and the category with the fewest genes was cell killing, with only 1 gene. In addition, as shown in Figure 5, in GO categories, biological processes related to reproduction and development included reproduction, reproductive process and developmental process, involving 451 corresponding genes.

The differentially expressed genes annotated in the NR database were analyzed in the KEGG database. The result showed that 446 genes were annotated into 175 KEGG pathways, in which 5 were significantly enriched pathways (p<0.05). As shown in Table 7, significantly enriched pathways are mainly involved in peroxisome proliferator-activated receptor (PPAR) signaling pathway, insulin signaling pathway, fructose and mannose metabolic pathways, GnRH signaling pathway, and transforming growth factor beta (TGF-β) signaling pathway.

To validate the reliability of Illumina analysis, we further used qRT-PCR to examine expression levels of five DEGs (acyl-CoA synthetase long-chain family member 3 [ACSL3], protein phosphatase 1 regulatory subunit 3A [PPP1R3A], myotubularin related protein 7 [MTMR7], follicle stimulating hormone beta [FSHβ], and SMAD specific E3 ubiquitin protein ligase 2 [SMURF2]) and in the significantly enriched pathways from peak and late stages of ovary. As shown in Table 8 and Figure 5, expression of these five genes can be detected, and the patterns of differential expression are consistent with the results from Illumina sequencing. This result suggested that the Illumina sequencing data was credible.

High-throughput transcriptomic sequencing technology has in recent years played an important role in revealing gene expression and regulation in animals and plants. It attracted considerable attention in studies focusing on gene expression and differential expression in particular aspects of animal development [15–17]. In this study, we used an Illumina HiSeq 2000 system for high-throughput sequencing of the ovarian transcriptomes from Shan Ma ducks at their peak and late stages of egg production, and obtained a great deal of transcriptional information. High-throughput sequencing results showed that sample quality was good and high-quality sequencing data Q20 accounted for 93% of the data. Transcriptomic sequencing data from ducks at their peak and late stages of egg production received 42,782,676 and 45,316,166 reads, respectively. The percentages of aligned reads to the reference genome were high, which were 72.70% and 71.33%, respectively. This suggests that the sequencing data have high coverage and the cDNA library was created successfully.

Analysis for differentially expressed genes between ovarian transcriptomes of ducks at their peak or late stages of egg production showed that there were 2,002 differentially expressed genes, with 790 genes upregulated and 1,212 genes downregulated. Further analysis of these genes in the NR database received 1,645 annotated genes, with 646 genes upregulated and 999 genes downregulated. There were 357 differentially expressed genes that did not have annotations. Genes that are highly expressed in reproductive tissues or cells may also indicate that the gene itself is activated, as shown in previous studies [18–21]. Therefore, in our study, we identified the top 10 up- and down-regulated DEGs in the ovaries of peak and late stage, respectively (Table 6). Some of them may be involved in specific reproductive processes. For example, previous research has shown that bone morphogenetic protein 15 gene is a member of TGF-β superfamily, and its function is to regulate the ovary development and to promote the proliferation and differentiation of follicular cells [22,23]. ATP5A1 gene involves in oxidative phosphorylation and ATP synthesis with an important role in aerobic respiration [24], and the study indicated that the expression level of ATP5A1 in the ovarian tissue during egg laying period was significantly higher than that in the early laying [25]. In oocytes, some related factors are coupled to the SLBP, which contains the only RNA binding region and interaction with stem-loop structure [26]. The pervious study showed SLBP played an important role in oocyte maturation, involving in translation regulation, processing and transport of histone mRNA [27,28]. The other most DEGs may play other important roles in biological processes. These DEGs may play an important role in regulating the function of the ovary and regulation of the egg production. So it should be further investigated the function of these genes.

From our GO analysis, the 1,645 genes with annotation in the NR database were analyzed in the GO database for functional classification. We observed that 1,204 of the genes received GO annotations and 441 did not. From our KEGG analysis, the 1,645 genes with annotations in the NR database were also analyzed in the KEGG database for pathway analysis. The 466 genes were annotated to KEGG pathways and 1,199 were not. We hypothesize that the lack of annotations of some genes was because studies on the duck transcriptome were still at early stages, and because there were not enough annotated genes available. Various types of biological information databases are constantly being updated and improved, and as a result, some of the sequencing data do not currently have corresponding annotations.

COG analysis results showed that 696 of the differentially expressed genes received COG functional annotations in 25 functional categories. This suggests that COG annotation is more comprehensive and covers the majority of life activities. GO analysis results showed that the differentially expressed genes involved 17 categories in the cellular component domain, of which cell-related (including cell and organelle) categories accounted for the majority, and membrane-related categories also accounted for a significant proportion. This result implies that gene products involved in the physiologic activities of the Shan Ma ovary are not only in the cells, but also widely present on membranes, which is consistent with studies on the ovarian transcriptomes of yaks [29], and suggests that the cell membrane plays an important role in the physiologic and reproductive activities of the ovary. GO analysis results showed that the differentially expressed genes also involved 17 categories in the molecular function domain, of which the binding category accounted for the greatest number. We further confirmed that binding proteins accounted for a high percentage in proteins expressed in oocytes. Other researchers have also shown that a high percentage of the proteins expressed in oocytes are binding proteins [30,31]. The fact that we confirmed that binding proteins were abundantly expressed in the ovary suggests that binding proteins play important roles in the physiologic and reproductive activities of the ovary. GO analysis results showed that the differentially expressed genes also involved 22 categories in the biological process domain, of which the cellular process category accounted for the largest percentage, followed by the biological regulation and metabolic process categories. In addition to the above biological processes, there were genes in the categories related to development and reproduction, including reproduction, reproductive process, and developmental process categories, involving 451 genes. This shows that these 451 reproduction- and development-related genes play different roles in Shan Ma ducks at their peak and late stages of egg-production.

In vivo, different genes are coordinated to perform biological functions. Pathway enrichment analysis will contribute to further comprehension and better annotation of differentially expressed genes, and can help to identify the main metabolic and signal-transduction pathways in which these genes are involved. In the present study, we performed KEGG pathway enrichment analysis, and our results showed that the differentially expressed genes identified involved 175 pathways, of which the PPAR signaling pathway, the insulin signaling pathway, fructose and mannose metabolism, the gonadotropin releasing signaling pathway, and the TGF-β signaling pathway are the 5 enriched pathways. Insulin is a polypeptide with multiple functions, and its binding to the insulin receptor can cause a series of cascading amplifications of its signal. It has been confirmed that insulin is an important factor in promoting the proliferation of ovarian granulosa cells and facilitating their physiologic functions, and its signaling pathway is involved in regulating granulosa cell proliferation, cellular metabolism and steroid hormone biosynthesis [32,33]. Previous study has showed that insulin signaling plays important roles in the physiologic and reproductive activities within the yak ovary; and is involved in cellular differentiation, development, and apoptosis, and in the integration of energy metabolism and reproductive activity in the ovary [29]. In our study, 14 of the differentially expressed genes we identified (suppressor of cytokine signaling 2 [SOCS2], lipase E [LIPE], forkhead box O1 a [FOXO1A], phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit gamma [PIK3CG], protein kinase cAMP-dependent type I regulatory subunit alpha [PRKAR1A], protein phosphatase 1 regulatory subunit 3B [PPP1R3B], phosphorylase kinase gamma 1 [PHKG1], hexokinase 2 [HK2], protein phosphatase 1 regulatory subunit 3A [PPP1R3A], phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit alpha [PIK3CA], phosphoenolpyruvate carboxykinase 2 [PCK2], SHC adaptor protein 3 [SHC3], 3-phosphoinositide dependent protein kinase 1 [PDPK1], and AKT serine/threonine kinase 2 [AKT2]), are involved in insulin signaling. GnRH is a neurohormone secreted by the hypothalamus that promotes the synthesis and release of gonadotropins by the pituitary (luteinizing hormone [LH] and follicle-stimulating hormone [FSH]), which in turn play important roles in animal reproduction [34], including sex hormone biosynthesis [35,36]. In this study, 5 of the differentially expressed genes we identified (follicle stimulating hormone beta subunit [FSHβ], phospholipase A2 group IIE [PLA2G2E], G protein subunit alpha q [GNAQ], mitogen-activated protein kinase kinase [MAP2K4], and calmodulin 3b [CALM3B]), are involved in GnRH signaling. TGF-β signaling is one of the most important signal-transduction pathways, and plays important roles in the development of various organs [37], and in regulating animal reproduction, especially ovarian functions [38,39]. We identified 4 differentially expressed genes (E1A binding protein p300 [EP300], protein phosphatase 2 scaffold subunit Aalpha [ppp2r1a], bone morphogenetic protein receptor, type 1A [BMPR1A], and SMAD specific E3 ubiquitin protein ligase 2 [SMURF2]) as being involved in TGF-β signaling. Collectively, we have confirmed that insulin-, GnRH-, and TGF-β-signaling pathways are important for physiologic and reproductive activities of Shan Ma ducks. In addition, other pathways such as the PPAR-signaling pathway and the fructose and mannose metabolic pathways were also enriched in our study. The functions of these two latter pathways in physiologic and reproductive activities within the ovary are not clear and require further elucidation.