Hypothalamic and ovarian transcriptome profiling reveals potential candidate genes in low and high egg production of white Muscovy ducks (Cairina moschata).

In China, the low egg production rate is a major challenge to Muscovy duck farmers. Hypothalamus and ovary play essential role in egg production of birds. However, there are little or no reports from these tissues to identify potential candidate genes responsible for egg production in White Muscovy ducks. A total of 1,537 laying ducks were raised; the egg production traits which include age at first egg (days), number of eggs at 300 d, and number of eggs at 59 wk were recorded. Moreover, 4 lowest (LP) and 4 highest producing (HP) were selected at 59 wk of age, respectively. To understand the mechanism of egg laying regulation, we sequenced the hypothalamus and ovary transcriptome profiles in LP and HP using RNA-Seq. The results showed that the number of eggs at 300 d and number of eggs at 59 wk in the HP were significantly more (P < 0.001) than the LP ducks. In total, 106.98G clean bases were generated from 16 libraries with an average of 6.68G clean bases for each library. Further analysis showed 569 and 2,259 differentially expressed genes (DEGs) were identified in the hypothalamus and ovary between LP and HP, respectively. The KEGG pathway enrichment analysis revealed 114 and 139 pathways in the hypothalamus and ovary, respectively which includes Calcium signaling pathway, ECM-receptor interaction, Focal adhesion, MAPK signaling pathway, Apoptosis and Apelin signaling pathways that are involved in egg production. Based on the GO terms and KEGG pathways results, 10 potential candidate genes (P2RX1, LPAR2, ADORA1, FN1, AKT3, ADCY5, ADCY8, MAP3K8, PXN, and PTTG1) were identified to be responsible for egg production. Further, protein-protein interaction was analyzed to show the relationship between these candidate genes. Therefore, this study provides useful information on transcriptome of hypothalamus and ovary of LP and HP Muscovy ducks.

INTRODUCTION

Muscovy ducks (Cairina moschata) are native to the tropical climates of Central and South America, Mexico, and southern Texas in the USA (Zeng et al., 2015). These ducks have been found in nearly 67 countries around the world. The ancestor of Muscovy duck is unique from closely all domestic duck breeds that are domesticated from Mallard (Anas platyrhynchos) (Veeramani et al., 2016). Muscovy duck has a lower egg production rate, and the number of eggs produced varies significantly among individuals (Zhu et al., 2020). The peak time of egg production of Muscovy ducks is between 35 to 53 wk, and 59 wk is the last stage of laying (Ye et al., 2017). In China, Muscovy ducks are raised on a large scale, but the low egg production rate has been a major constraint for the farmers.

Egg production performance is a polygenic trait with a low heritability value (Yu et al., 2015), which is affected by various organs involved in the reproduction (Su et al., 1999). The hypothalamic-pituitary-gonadal (HPG) axis regulates the reproductive and endocrine systems of laying birds (Padmanabhan et al., 2002). In birds, the hypothalamus controls reproduction by producing gonadotropin-releasing hormone-1 (GnRH-1) that stimulates the release of follicle stimulating hormone (FSH) and luteinizing hormone (LH) in the pituitary gland, which is responsible for the growth of ovarian follicles in the ovary during ovulation. The ovaries are pair-like of ova-producing organs (produces egg cells) that maintain the distinction of being an endocrine gland due to the secretion of hormones (estrogen and progesterone), which is vital for normal reproductive development and fertility. In a study by Liu et al. (2021), insulin-like growth factor-1 (), prolactin receptor (PRLR), follicle-stimulating hormone (FSH), very low-density lipoprotein receptor (), low-density lipoprotein receptor-related protein (), estrogen receptor 1 (), and Luteinizing hormone receptor () are essential genes for egg production traits of Pekin duck and Black Muscovy duck.

RNA sequencing (RNA-seq), a method based on next-generation sequencing (NGS), offers the opportunities to provide ample details about the transcriptional landscape of a particular organism with the usage of advanced molecular biology techniques and bioinformatics (Wang et al., 2008). Recently, different studies have used this method for identification of genes or molecular pathways responsible for egg production in different species of poultry such as in chicken (Wang and Ma, 2019), Jinghai Yellow chicken (Zhang et al., 2019), Luhua Chicken (Wu et al., 2017), Turkey (Brady et al., 2020), Black and White Muscovy duck (Bao et al., 2021), Leizhou Black duck (Zou et al., 2020), and Goose (Ouyang et al., 2020). However, no reports have used hypothalamus and ovary tissues of White Muscovy ducks to identify potential genes responsible for egg production.

In this study, we identified potential candidate genes through high-throughput RNA sequencing technology by generating a comprehensive transcriptome profile of hypothalamus and ovary tissues of selected lowest and highest producing White Muscovy ducks at 59 wk of age.

MATERIALS AND METHODS

Ethics Statement

In this study, the entire animal used for experiments were approved by the South China Agricultural University Institutional Animal Care and Use Committee (Guangzhou, People's Republic of China), and were handled strictly in compliance with the guidelines of this committee.

Experimental Animals and Sample Preparation

A total of 1,537 Muscovy laying duck raised in Muscovy duck breeding farm of Guangdong Wenshi Southern Poultry Breeding Co., Ltd Renma Farm were used in this study. All the ducks were housed in an individual cage with the same feeding and management conditions. The number of eggs was monitored and recorded from age of first egg to 59 wk of age. The ducks were classified at 300 d of age (EN300) into 2 groups; low producing and high producing Muscovy ducks. The 4 lowest (LP) and 4 highest producing (HP) Muscovy ducks within the same egg number from each group were selected as final experimental ducks. The reproduction parameters such as age at first egg (AFE), number of eggs at 300 d (NE300), and number of eggs at 59 wk (NE59) were recorded and analyzed for the LP and HP ducks. The selected ducks were slaughtered, and the hypothalamus and ovary (excluding both the white and yellow follicles) were collected (a day after withdrawing of feed) quickly. All tissues were washed with RNA-free water, wrapped in nylon polybags and frozen in liquid nitrogen, and then stored at –80°C until RNA extraction.

RNA Extraction and Quality Control

Total RNA was extracted from all samples using Trizol reagent (Ambion, Austin, TX) according to the manufacturer's instructions. RNA purity was measured using the Nanodrop2000 spectrophotometer. RNA integrity was assessed using the RNA 6000 Assay kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA).

Library Preparation and Transcriptome Sequencing

A total of 1 μg RNA per sample was used for the construction of 16 libraries. Briefly, purification of mRNA and cleavage of the mRNA into fragments were carried out. Reverse transcription and cDNA synthesis were also done, and clustering of the index-coded samples was performed according to the AMPure XP system (Beckman Coulter, Brea, CA) instructions. After cluster generation, the library preparations were sequenced on an Illumina Novaseq platform at Novogene Biotech. Co. Ltd. in Beijing, China.

Transcriptomic sequencing Analysis

The raw reads (fastq format) were firstly processed through in-house perl scripts. After this, clean data (clean reads) were obtained by removing reads containing adapter, reads containing poly-N, and low quality reads from raw data. Simultaneously, Q20, Q30, and GC content of the clean reads were calculated. Clean reads were then mapped to the Mallard duck (Anas platyrhynchos) using Hisat v2.05 with default settings. The mapped reads of each sample were assembled by StringTie (v1.3,3b).

Identification of Differentially Expressed Genes (DEGs)

The reads numbers mapped to each gene were counted using features Counts v1.5.0-p3 (Liao et al.) based on the expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced (FPKM) method. The differential expression analysis between the LP and HP hypothalamus and ovary samples was performed using the DESeq2 R package (1.20.0) (Love et al., 2014). DESeq2 provides statistical routines for determining differential expression in digital gene expression data using a model on the negative binomial distribution. The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling false discovery rate. Genes with P-value <0.05 found by DESeq2 were assigned as differentially expressed.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genome (KEGG) Pathway Enrichment Analyses of DEGs

GO enrichment analysis of DEGs was conducted by the clusterProfiler R package, in which gene length bias was corrected. Pathway enrichment analysis was performed using the KEGG database (http://www.genome.jp/kegg/). GO terms and KEGG pathways with adjusted P-value <0.05 were considered significantly enriched among DEGs.

Protein-Protein Interaction

The STRING 10 database (http://string-db.org/) was employed to identify the relationship between the candidate genes identified in this study. All the network visualization was performed using Cytoscape_v3.2.1 (Shannon et al., 2003).

Validation of RNA-seq Results Using qRT-PCR

To validate RNA-seq results’ reliability, quantitative real-time PCR (qRT-PCR) was done on 9 selected DEGs. Total RNA was extracted from the same hypothalamus, and ovary samples used for the RNA-seq using MagZol Reagent (Magen, Guangzhou, China), and 1 μg total RNA was reversely transcribed to cDNA using the Evo M-MLV reagent Kit (Hunan, China) according to manufacturers’ instructions. The MonAmp ChemoHS qPCR mix (Monad, Suzhou, China, Cat# MQ00401) was used to determine the expression level of mRNA in the 9 selected genes using GAPDH as the house-keeping gene. The qRT-PCR reaction was performed in CFX96 (Bio-Rad, Hercules, CA). The primers were designed according to the NCBI database sequence using Primer premier version 5.0 software (Applied Biosystems). 10.9 μL qPCR total mixed volume was used in each well which contained 5 μL MonAmp ChemoHS qPCR mix (Monad, Suzhou, China, Cat# MQ00401), 0.2 μL (10 μmol/L) each forward and reverse primer, 1.5 μL cDNA template and the rest was RNase-free H2O. The PCR reaction was conducted as follows: 30 s at 95°C, followed by 40 cycles each for 5 s at 95°C, and annealing temperature of each primer for 34 s and 30 s at 72°C. All reactions were performed in triplicate. The tissues of HP Muscovy ducks were used as the control. The 2−ΔΔCT method was used to calculate target gene expression (Livak and Schmittgen, 2001).

The expression value of each gene from qRT-PCR and RNA-seq were log2 transformed. A linear regression line and correlation coefficient were obtained between the results of RNA-Seq and qRT-PCR. Primer sequences of the genes are provided in Supplementary Data S9.

Statistical Data Analysis

Data were expressed as the mean ± standard deviation (SD). The significance of the difference of the egg production traits, linear regression, and correlation coefficient between the log2RNA-Seq and log2qRT-PCR was done using two-tailed Student's t test with SPSS 19.0 statistical software. P-value <0.001 was considered statistically significant.

RESULTS

Egg Production Traits of the Highest and Lowest Laying Ducks

The egg production traits of the LP and HP Muscovy ducks are shown in Supplementary Data S1. The results showed that age at first egg (days) was not significantly different (P > 0.001) between the 2 groups. However, there was a significant difference (P < 0.001) between the number of eggs at 300 d and number of eggs at 59 wk of LP and HP.

Transcriptome Data

The results of the transcriptome sequencing and mapping information are shown in Table 1 and Supplementary Data S2. Both the raw read and clean read of each library were more than 40 million. The clean reads were mapped to Mallard duck (Anas platyrhynchos) genome with a total map ranging from 50.83% to 61.06%. The average GC content of the samples was approximately 50%, while the average percentages of Q20 and Q30 bases were 97.15% and 92.90%, respectively. Supplementary Figure 1 shows the correlation coefficient was above 0.85 among the hypothalamus and ovary samples, which affirms the reliability of the mapping data for further analyses. The sequencing raw data have been deposited into the Sequence Read Archive in NCBI with accession number PRJNA704544.

Table 1

Quality analyses of transcriptome sequencing and mapping.

Sample	Raw_reads	Clean_reads	Clean base	Error rate (%)	Q20 (%)	Q30 (%)	GC content (%)	Total_mapped
HY1	45700082	44305268	6.65G	0.03	97.15	92.90	50.33	25173764 (56.82%)
HY2	41602714	40516548	6.08G	0.03	97.17	92.84	49.07	21400238 (52.82%)
HY3	43702420	42494908	6.37G	0.03	97.25	92.98	48.77	25307792 (59.55%)
HY4	45303344	44159426	6.62G	0.03	96.90	92.40	49.24	22945033 (51.96%)
HY5	46745386	45503158	6.83G	0.03	97.17	92.86	47.42	26743235 (58.77%)
HY6	44801406	43548168	6.53G	0.03	97.15	92.85	48.49	23128383 (53.11%)
HY7	47309352	46594362	6.99G	0.03	96.81	92.52	50.23	23685077 (50.83%)
HY8	46110576	43836962	6.58G	0.03	97.22	93.19	50.32	23937415 (54.61%)
OVARY_00289	47892996	46572178	6.99G	0.03	96.96	92.25	50.09	28032233 (60.19%)
OVARY_00318	45305832	43867458	6.58G	0.03	97.03	92.70	50.97	25470584 (58.06%)
OVARY_00413	55564892	53629668	8.04G	0.03	97.17	92.99	52.15	30966869 (57.74%)
OVARY_04703	41701644	40420626	6.06G	0.03	97.09	92.79	51.23	23105666 (57.16%)
OVARY_04727	45492404	43697182	6.55G	0.03	97.42	93.37	50.55	26688404 (61.08%)
OVARY_04975	46468602	44615024	6.69G	0.03	97.53	93.68	50.62	27374192 (61.36%)
OVARY_05028	47224738	45140604	6.69G	0.03	97.26	93.16	50.81	26260799 (58.18%)
OVARY_05711	45601332	44876074	6.73G	0.03	97.04	92.86	49.65	25186802 (56.13%)

Note: HY1-HY8 represents the hypothalamus tissues. Ovary_00289-Ovary_05711 represents the Ovary tissues. Raw reads- the statistics of the number of sequencing sequences generated from library construction. Clean reads-this is statistics of the number sequencing sequences filtered by sequencing data. Q20 and Q30 represent the percent of sequenced bases that have a predicted quality score of 20 and 30, respectively. QC content means the dependence between fragment counts (read coverage). Total mapped describes the statistics of the sequencing sequences that can be located on the reference genome.

Figure 1

Volcano plot of differentially expresses genes in (A) hypothalamus (B) ovary. Red dots represent the up-regulated genes, and green dots show down-regulated genes. Blue dots represent genes that did not show obvious changes between lowest-producing (LP) and highest-producing (HP) samples. (C) The Venn diagram shows the number of genes expressed in the two tissues.

Identification of Differential Expressed Genes

Based on the 2 tissues used for the sequencing, genes with significant differences (P < 0.05) at expression levels were called differentially expressed genes (DEGs). In this study, the sequencing data of the LP and HP Muscovy ducks were examined. In the hypothalamus, a total of 596 DEGs were found, including 227 up-regulated and 369 down-regulated genes, whereas 32 were novel genes (Figure 1A, Supplementary Data S3). Also, a total of 1,502 DEGs were identified in the ovary comprising 745 up-regulated and 757 down-regulated genes, of which 160 genes were novel (Figure 1B, Supplementary Data S4). Based on the identified DEGs in the hypothalamus and ovary samples, 99 genes were common in both tissues as shown in the Venn diagram (Figure 1C). The FKPM hierarchical clustering map (Figure 2) expresses the gene expression patterns between the samples and presents the reliability of the gene set.

Figure 2

Hierarchical clustering of Differentially Expressed Genes. HP_HY and LP_HY represent the hypothalamus samples of highest and lowest producing Muscovy ducks, respectively. HP_Ovary and LP_Ovary represent the ovary samples of highest and lowest producing Muscovy ducks, respectively. The red colors represent the over expressed genes, and the blue colors represent the genes with the lowest expression levels. Colored bars indicate the expression level.

Gene Ontology and KEGG Pathway Enrichment Analyses of Differentially Expressed Genes

To further express the functions of the identified DEGs between the hypothalamus and ovary of the LP and HP, analysis of GO and KEGG pathways was done. 435 GO terms comprising of 207 biological processes (BP), 46 cellular components (CC), and 182 molecular functions (MF) were identified in the hypothalamus. Only 24 BP which includes movement of cell or subcellular component, regulation of cellular metabolic process, regulation of transcription, DNA membrane, and 21 MF which includes microtubule motor activity, cytoskeletal protein binding, molecular transducer activity were enriched at adj.P < 0.05 (Figure 3A, Supplementary Data S5). Moreover, 676 terms were identified in the ovary comprising of 355 BP, 75 CC, and 246 MF. 18 enriched BP including immune system response, G-protein coupled receptor signaling pathways, regulation of cell growth, 5 CC including extracellular matrix, extracellular region, plasma membrane, and 32 MF including signaling receptor activity, cytokine receptor binding, transcription regulator activity in the ovary were enriched at adj.P < 0.05 (Figure 3B, Supplementary Data S6). The enriched DEGs in the hypothalamus are primarily affected by ion channel activity, substrate-specific channel activity process, passive-transmembrane transporter activity while signal transducer activity, peptidase activity, and metalloendopeptidase activity affected the enriched DEGs in the ovary (Figure 3C,D).

Figure 3

GO-enriched bar graph of differentially expressed genes in (A) HP_HY and LP_HY (B) HP_Ovary and LP_Ovary. The red bar indicates biological processes; green indicates cellular components; blue indicates molecular functions. The scatter plot shows 30 GO terms in (C) HP_HY and LP_HY (D) HP_Ovary and LP_Ovary. The ordinate indicates the GO terms, and the abscissa indicates the Rich factor. The size of the point indicates how many differentially expressed genes are in the GO terms, and the color of the point corresponds to the p.adj value range.

Based on the KEGG pathway enrichment analysis, 114 pathways were identified in the hypothalamus, 8 were enriched which includes Calcium signaling pathway, Neuroactive ligand-receptor interaction, RNA degradation at adj.P < 0.05 (Figure 4A).139 pathways were identified in the ovary, only 11 were enriched including ECM-receptor interaction, Focal adhesion, MAPK signaling pathway (Figure 4B). Apoptosis and Apelin signaling pathways were enriched in both hypothalamus and ovary samples (Supplementary Data S7 and S8).

Figure 4

The top 30 KEGG enrichment scatter plot of differentially expressed genes in (A) the hypothalamus (B) ovary samples. The ordinate indicates the GO terms, and the abscissa indicates the Rich factor. The size of the point indicates how many differentially expressed genes are in the KEGG pathways, and the color of the point corresponds to the p.adj value range.

Identification of Candidate Genes

Considering the significance levels of the DEGs, GO and KEGG pathway results and literature reviews suggest P2RX1, LPAR2, ADORA1, FN1, AKT3, ADCY5, ADCY8, MAP3K8, PXN, and PTTG1 as potential candidate genes for egg production of Muscovy ducks. The details of these identified genes in HP_HY and LP_HY, and HP_Ovary and LP_Ovary are listed in Table 2. Further, to generate the relationship between the potential candidate genes, protein-protein network based on String 10 database was shown using cytoscape (Figure 5).

Table 2

List of potential candidate genes involved in egg production in hypothalamus and ovary samples of the LP and HP Muscovy duck.

Tissues	Gene symbol	Expression level in LP	Expression level in HP	Gene description		P-value	Log2FC	GO terms	KEGG pathways	Expression level
Hypothalamus	P2RX1	22.81	47.87	Purinergic receptor p2 × 1	0.0383		1.0700	Signaling receptor activity	Calcium signaling pathway	Up-regulated
	LPAR2	42.93	21.98	Lysophosphatidic acid receptor 2	0.0154		-0.9657	Transmembrane signaling receptor activity	Regulation of actin cytoskeleton	Down-regulated
	ADORA1	457.73	222.51	Adenosine a1 receptor	0.0163		-1.0409	Signal transducer activity	Neuroactive ligand-receptor interaction	Down-regulated
Ovary	FN1	1520.25	3994.12	Fibronectin 1	0.0178		1.3936	Extracellular region	ECM-receptor interaction	Up-regulated
	AKT3	270.51	554.83	Akt serine/threonine kinase 3	0.0239		1.0365	Protein serine/threonine kinase activity	Apoptosis	Up-regulated
	ADCY5	409.98	1017.12	Adenylatecyclase 8	0.0239		1.0365	Plasma membrane	Apelin signaling pathway	Up-regulated
	ADCY8	85.08	237.06	Adenylatecyclase 8	0.0001		1.4818	Plasma membrane	Apelin signaling pathway	Up-regulated
	MAP3K8	41.65	95.78	Mitogen-activated protein kinase kinase kinase 8	0.0043		1.2023	Phosphotransferase activity	MAPK signaling pathway	Up-regulated
	PXN	25.19	68.08	Paxillin	0.0348		0.5535	Zinc ion binding	Focal adhesion	Up-regulated
	PTTG1	2001.78	992.06	Pituitary tumor-transforming 1	0.0457		-1.0124	Cellular component organization	Oocyte meiosis	Down-regulated

HP, highest producing ducks; LP, lowest producing ducks.

Figure 5

Protein-protein interaction of the relationship between the potential candidate genes in the hypothalamus and Ovary of the LP and HP Muscovy ducks. The yellow colors show the genes that are up-regulated, while the red colors indicate the down-regulated genes. The highlighted genes in blue are DEGs but not candidate genes predicted to be involved in the network.

Validation of RNA-Seq results using qRT-PCR

To affirm the reliability of the RNA-Seq results, 9 DEGs including 6 up regulated genes (GHRH, ZAR1L, GHSR, BRS3, RGS1, MMP9) and 3 down regulated genes (FSHR, HTR5A, CDC25A) were selected for quantitative reverse transcription polymerase chain reaction (qRT-PCR). The results of the comparison of the mRNA expression levels of these genes normalized with GAPDH, and RNA-Seq results are presented in Figure 6A. The linear regression line and Pearson correlation coefficient (r) of 0.9298 affirm a strong relationship between the mRNA expression levels from qRT-PCR and RNA-Seq as presented in Figure 6B. These results show the consistency between the qRT-PCR and RNA-Seq.

Figure 6

The graphs show (A) comparison of the mRNA expression levels (B) linear regression between RNA-Seq and qRT-PCR.

DISCUSSION

Egg Production and Analysis of the Egg Production Traits Between LP and HP Ducks

Egg production is a trait with low heritability, which is affected by multiple genes (Yu et al., 2015). This trait is also regulated by hormonal regulations in the hypothalamus-pituitary-gonadal (HPG) axis (Mellouk et al., 2018). In poultry industries, the number of eggs produced and the rate of egg-laying are major estimators of laying performance (Tao et al., 2017), which vary significantly among individuals (Zhu et al., 2020).

In this study, the egg production traits which include age at first egg (days) (AFE), number of eggs at 300 d (NW300), and number of eggs at 59 wk (NW59) of the LP and HP groups were analyzed. Based on the analysis of these traits between the 2 groups, NW300 and NW59 of the HP ducks were significantly higher (P < 0.001) than the LP counterpart. The differences between these groups might be due to the number of follicles present in the LP and HP ducks which are closely related to the secretion of important reproductive hormones (LH and FSH) as reported by a study that the higher the number of follicles, the greater the number of egg produced in laying geese (Yang et al., 2019). Moreover, a transcriptome study in Peking ducks revealed that variation in egg-laying might be due to the differences in the development of ovaries and follicles which may be closely associated with the ovarian steroid synthesis pathway (Ren et al., 2019).

Transcriptomic Sequencing and Analysis of the Differentially Expressed Genes Between LP and HP

Since the invention of RNA-sequencing technology, several studies have been carried out to provide ample information about the mechanism between LP and HP in laying birds. In this study, the transcriptomic study of the hypothalamus and ovary of LP and HP was done to unveil potentials DEGs that are important for egg production in White Muscovy duck, as this would help in the improvement of the breeding of high producing ducks to increase the availability of quality egg to the human populace and improve the economic value of the farmers (Zeng et al., 2015; Zhu et al., 2017). Based on the transcriptomic studies, 596 DEGs and 1,502 DEGs were identified between LP and HP Muscovy ducks of the hypothalamus and ovary samples, respectively. The higher number of DEGs in the ovary compared with the hypothalamus follows the result of a study where 414, 356, and 10 DEGs were identified in the pituitary gland, ovary, and hypothalamus of high egg production (HEP) and low egg production (LEP) chicken (Mishra et al., 2020). This means that the ovary of poultry birds usually have more DEGs than the hypothalamus which might be due to the multitude effects of different genes that tightly control ovarian follicle development and hormone secretion (Quan et. al., 2019; Zou et al., 2020). To confirm the reliability of the transcriptomic sequencing results, 9 DEGs were selected for quantitative real-time PCR (qRT-PCR), and the expression levels of these genes and regression coefficient proved that the 9 genes were in agreement with the RNA-seq data.

GO and KEGG Analysis of the DEGs

To further explain the biological roles of the DEGs, Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genome (KEGG) analysis were performed. The GO results shown that biological processes and molecular functions of the hypothalamus and ovary were assigned more DEGs than the cellular components. Moreover, the DEGs in the hypothalamus were more involved in GO terms related to the transcription of DNA template, regulation of the cellular biosynthetic process, nucleic acid template transcription, and regulation of nitrogen compound and primary metabolic processes. This agrees with a report that genes identified in their study were primarily related to biosynthetic processes, cellular nitrogen metabolic processes, transport, cell differentiation, cellular protein modification processes, signal transduction, and small molecule metabolic processes during prelaying and laying periods in Sichuan White geese (Gao et al., 2015). This suggests that transcription and metabolic processes are very essential in the reproductive activities of White Muscovy ducks. In the ovary, 676 GO terms were revealed with 355 terms under biological processes of which the “protein modification process” had the highest number of DEGs followed by the “phosphorus metabolic” process and “phosphate-containing compound metabolic process.” The highest number of DEGs in the biological processes in this study agrees with a report that fertilization of eggs is a complex biological activity that requires different changes in the conversion of both cellular and molecular (Wang et al. 2021). Moreover, in this study, “extracellular region” and “membrane-bounded organelle” were the cellular components with the highest DEGs in the ovary. In a study on proteomic profiles of human sperm, GO analysis revealed that most of the differentially identified sperm protein was enriched in the extracellular membrane bounded organelles, cytoplasm, and cytosol (Li et al., 2019). The molecular functions in the ovary revealed that the highest DEGs were recorded in “transferase phosphorus containing group” followed by “phosphotransferase activity alcohol group as acceptor.” It has been reported that transferase are highly essential in a wide range of physiological activities by catalyzing different biological processes that are important to a living system (Fyte et al., 2012). Moreover, phosphotransferase, also known as phosphorus transferases is a large enzyme category among all transferases (Su et al., 2013). The activities of transferase include but are not limited to essential metabolism, genetic regulations, and detoxification (Fyte et al., 2012).

KEGG pathway enrichment analysis was done to better understand the different metabolic and signaling transduction pathways the genes are involved. 114 pathways were revealed in the hypothalamus with “Neuroactive ligand-receptor interaction” had the highest DEGs followed by “calcium signaling pathway” and “apoptosis.” This result is similar to a study that confirmed that the neuroactive ligand-receptor interaction pathway is essential for egg production in geese (Ouyang et al., 2020). Also, calcium (Ca2+) is an important signaling molecule that influences a diverse range of biological activities (Park et al., 2020). Moreover, the calcium signaling pathway was found correlated with eggshell quality in chicken (Zhang et al. 2015a), and this pathway was also enriched in a study on egg production in duck ovaries (Tao et al., 2017). In the ovary, 139 KEGG pathways were involved. Although, the “neuroactive ligand-receptor interaction pathway” had more DEGs followed by “MAPK signaling pathway” and “focal adhesion”. It has been reported that MAPK (mitogen-activated protein kinase) signaling pathways regulate a diversity of biological processes through several cellular mechanisms (Yue and López, 2020). In a recent study, the MAPK signaling pathway was identified to be involved in the egg production of hen at different laying periods (Zheng et al., 2021). The high DEGs in the focal adhesion in this study agree with a report on egg production in Jinghai Yellow chicken in which 4 known reproductive pathways which include calcium signaling, WNT signaling pathway, focal adhesion, and cytokine-cytokine receptor interaction pathways were detected (Zhang et al., 2019). Also, several reports have indicated the importance of calcium signaling pathway, WNT signaling pathway, focal adhesion, and cytokine-cytokine receptor interaction pathway on egg production and egg-shell quality in chickens (Chen et al., 2007a, 2007b; Zhang et al., 2015a, 2015b, 2015c; Tao et al., 2017).

Identified Potential Candidates Genes

P2X1 is one of the 7 different subtypes of P2X receptors, specifically activated by ATP (Kawate et al., 2009; Jiang et al., 2013). It has been reported that double deletion of P2X1 and α1A-adrenergic G protein–coupled receptors in male mice could lead to infertility due to inhibition of sperm transport (White et al., 2013). The low expression of P2X1 in females might affect the ovulation of eggs. Purinergic signaling can cause an increased influx of Ca2+ from the extracellular environment and the release of Ca2+ intracellular stores (Pérez-Flores et al., 2015). Calcium signaling plays a role in the reproduction and development of mammals (Stewart and Davis, 2019). In this study, calcium signaling pathway was enriched with the P2XI as one of the genes associated with it. Therefore, this shows that P2XI might be a potential for egg production in White Muscovy duck.

LPAR2 is one of the 6 receptors of Lysophosophatidic acid (LPA) that exists in different mammalian cells and tissues (Lin et al., 2010); that affects cell migration, survival, cellular interaction, and changes in the cytoskeleton (Moolenaar et al., 2004). High expression of LPAR2 has been reported in testis, kidney, heart, lung, and brain (Hecht et al., 1996; Contos and Chun, 2001). In cattle, pigs, and sheep, LPA plays a great role in the influence of the uterine cycle during the estrous cycle and early pregnancy; majorly implantation via its role on endometrial prostaglandin (PG secretion) (Seo et al., 2008; Liszewska et al., 2009; Woclawek-Potocka et al., 2009). As reported by (Aoki et al., 2008; Smyth et al., 2008; Ye, 2008), LPA is present in different biological fluids including chicken's egg white, thereby suggesting the role of LPAR2 on egg production.

ADORA1 is one of the G-protein coupled receptor family members which includes A1, A2A, A2B, and A3 receptors (Haskó et al., 2008), which is highly distributed in the brain (Sebastião and Ribeiro, 2009). ADORA1 has been reported to play an important role in cell proliferation and hormone secretion (Wang et al., 2015). This gene plays a key role in fertilization by obstructing adenylyl cyclase (a messenger that regulates cell growth and differentiation) and has been identified as a candidate gene for sexual precocity in goats (Su et al., 2018). Based on the importance of this gene in reproductive processes, it can be affirmed that ADORA1 play a significant role in the rate of egg production in White Muscovy ducks.

FN1 gene provides information for developing 2 types of fibronectin-1 protein: soluble plasma fibronectin-1 and insoluble cellular fibronectin-1. This gene is involved in several biological processes. Its major role is the support of cell adhesion (Goossens et al., 2009). However, it is also involved in cytoskeleton organization, cell migration, and important physiological processes (Yamada, 1991; Moursi et al., 1996). FN 1 is a major component of the extracellular matrix (ECM) (Lemańska-Perek and Adamik, 2019). Also, this gene has been identified as a candidate gene in bovine preimplantation embryos (Goossens et al., 2009). Due to the importance of FN1 in cell adhesion and extracellular matrix; which are both enriched pathways identified in this study, it can be inferred that this gene can affect egg-laying performances of ducks.

AKT3 gene is one of the 3 isoforms of AKT. The AKT gene family plays a great role in the diversity of genetic processes that include cell propagation, apoptosis, development, and angiogenesis among others (Hennessy et al., 2006; Manning and Cantley, 2007). AKT plays a great role in bovine early embryonic development (Ashry et al., 2018). This gene is a fundamental regulator of many cytokines- and hormone-driven processes (Fabi and Asselin, 2014). Moreover, it is an important gene for maintaining brain homeostasis (Easton et al., 2005). In this study, apoptosis pathway was enriched in both the hypothalamus and ovary. Based on this and the importance of AKT3 in several physiological processes, AKT3 can be a potential gene that affects egg production performance in duck.

ADCY5 gene is a member of the adenylatecyclases (ACs) family, which contains 12 domain transmembrane proteins responsible for catalyzes the conversion of ATP to cyclic AMP (Defer et al., 2000). Out of the 12 domains, AC5 is prominent for controlling the extra pyramidal motor system related to several neuropsychiatric diseases, such as Parkinson's disease (Chen et al., 2012; Procopio et al., 2013). The variants of ADCY5 gene were affirmed to be related to human gestational duration (Zhang et al., 2017; Zhou et al., 2020), and fetal growth characteristics such as birth weight of the human fetus (Freathy et al., 2010). Also, ADCY5 gene has been proven to affect the development of ovarian morphological-related traits in Bovine (Li et al., 2021). Considering the importance of ADCY5 on reproductive processes, it can be affirmed that this gene might affect egg production parameters in White Muscovy ducks.

ADCY8 constitutes only a minor adenylyl cyclases (Delmeire et al., 2003; Waddleton et al., 2008), which is stimulated by Ca2+ in a calmodulin (CaM)-dependent manner (Cooper et al., 1995) that helps in support of its distinct physiological roles (Masada et al., 2009). This gene is highly expressed in the thalamus and the cerebral cortex (Nicol et al., 2005). In a study on Zebrafish (Danio rerio), mRNA expression of ADCY8 was identified as one of the genes that control the endocrine system in the HPG axis (Liang et al., 2020). Due to the ability of ADCY8 to maintain the HPG axis which is a highly essential axis in reproduction, it can be inferred that ADCY8 plays a key role in the rate of egg-laying in ducks.

MAP3K8 (also known as COT) is one of the serine/threonine protein family members. It is important in cell transformation, proliferation, migration, and invasion by activating extracellular signal-regulated kinase (ERK), Rac1, and focal adhesion kinase (FAK) (Dumitru et al., 2000; Das et al., 2005). MAP3k8 was predicted as a target gene of miRNA-509-3p that promotes the secretion of estradiol (Huang et al., 2016), which is an important hormone in egg ovulation. Besides, MAP3K8 affects several signaling pathways which include MAPK pathway in a cell type- and stimulus specific manner (Das et al., 2005). In this study, MAPK pathway is associated with MAP3K8 that can affect the egg production of Muscovy ducks.

PXN is one of the family genes of paxillin along with Hic-5 (TGFb1I1) and leupaxin (LPXN). PXN is a primary component of focal adhesions that are involved in cell movement, migration, and embryogenesis (Deakin and Turner, 2008; Hu et al., 2014). Moreover, the lack of expression of PXN during mouse embryonic development is deadly (López-Colomé et al., 2017). This was supported by the high expression of focal adhesion kinase (FAK) and PXN in the hypothalamus of female rats at birth compared with male rats (Gozin et al., 1998). It can be inferred that the low expression of PXN in LP ducks compared to HP ducks might be one of the reasons for low egg production in the LP ducks at NE300 and NE50 thereby, affirming the importance of this gene in the egg production of ducks

PTTG1 (also known as securin) is an oncogene that plays a key role in many biological processes, which includes organ development, apoptosis, DNA repair, regulation of cell cycle, and separation of sister chromatid (Panguluriet al., 2008; Moreno-Mateos et al., 2011). This gene expression has been found in germ, Leyding, and sertoil cells in testis (Pei, 1998). Moreover, there is an increase in the expression levels of PTTG1 during cell proliferation and in mitosis, which shows an important role in controlling the cell cycle (Ramos-Morales et al., 2000). PTTG1 has been proven to be immensely expressed in several tumors, including pituitary, uterine, and ovarian (Tong and Eigler, 2009; Chen et al., 2017). In this study, the oocyte meiosis pathway was related to PTTG1, this confirms the tendency of this gene to affect egg production of Muscovy ducks.

In conclusion, the hypothalamus and ovary are important tissues along the HPG axis that play a major role in egg production. In this study, transcriptome sequencing of these 2 tissues was done. Based on the percentage of Q30 of the samples, we assumed that the low mapping rate in this study was due to the unavailability of the reference genome of Muscovy ducks (Cairina moschata). A total of 596 DEGs and 1,502 DEGs were identified in the hypothalamus and ovary, respectively, between lowest and highest producing White Muscovy ducks. Ten potential candidate genes (P2RX1, LPAR2, ADORA1, FN1, AKT3, ADCY5, ADCY8, MAP3K8, PXN, and PTTG) based on their corresponding GO terms and pathways were identified. These potential candidate genes identified in this study can be used as selection markers to improve egg production in White Muscovy ducks.