Transcriptomic Analyses of the Hypothalamic-Pituitary-Gonadal Axis Identify Candidate Genes Related to Egg Production in Xinjiang Yili Geese

Simple Summary Transcriptome is the study of gene expression at the RNA level. In this study, high-throughput transcriptome sequencing technology was used to compare and analyze the gene expression differences in the hypothalamus, pituitary, and ovary tissues of Xinjiang Yili geese with high and low egg production. We preliminarily screened 30 candidate genes related to egg production regulation of Xinjiang Yili geese, including the calcium signaling pathway. Therefore, this study laid a theoretical foundation for further elucidating the molecular mechanism of egg production of Yili geese in Xinjiang. Abstract The study was conducted to investigate the transcriptomic differences of the hypothalamic-pituitary-gonadal axis between Xinjiang Yili geese with high and low egg production and to find candidate genes regulating the egg production of Xinjiang Yili geese. The 8 selected Xinjiang Yili Geese with high or low egg production (4 for each group) were 3 years old, with good health, and under the same feeding condition. High-throughput sequencing technology was used to sequence cDNA libraries of the hypothalami, pituitary glands, and ovaries. The sequencing data were compared and analyzed, and the transcripts with significant differences were identified and analyzed with bioinformatics. The study showed that the transcriptome sequencing data of the 24 samples contained a total of 1,176,496,146 valid reads and 176.47 gigabase data. Differential expression analyses identified 135, 56, and 331 genes in the hypothalami, pituitary glands, and ovaries of Xinjiang Yili geese with high and low egg production. Further annotation of these differentially expressed genes in the non-redundant protein sequence database (Nr) revealed that 98, 52, and 309 genes were annotated, respectively. Through the annotations of GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) databases, 30 candidate genes related to the egg production of Xinjiang Yili geese were preliminarily selected. The gap junction, focal adhesion, and ECM-receptor interaction signaling pathways were enriched with the hypothalamic, pituitary, and ovarian differentially expressed genes, and the calcium signaling pathway was enriched with the pituitary and ovarian differentially expressed genes. Thus, these pathways in the hypothalamic-pituitary-gonadal axis may play an important role in regulating egg production of Xinjiang Yili geese. The results provided the transcriptomic information of the hypothalamic-pituitary-gonadal axis of Xinjiang Yili geese and laid the theoretical basis for revealing the molecular mechanisms regulating the egg-laying traits of Xinjiang Yili geese.


Total RNA Extraction
The total RNA was extracted using Trizol (Invitrogen, Carlsbad, CA, USA) by following the manufacturer's protocol, and 4 measures were used to ensure the quality of the samples was good enough for transcriptome sequencing. First, agarose gels were used to analyze the extent of RNA degradation and whether the samples were contaminated. Then, RNA was assessed for quantity and quality using a Nanodrop ND-1000 spectrophotometer (Implen, Weslake Village, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) according to the manufacturer's instructions.

Construction of cDNA Libraries and Illumina Sequencing
After the samples passed the quality control, the cDNA libraries were constructed according to the instruction of the Illumina kit (Illumina, San Diego, CA, USA). Next, we used Qubit 2.0 for the preliminary quantification of the constructed libraries and diluted the libraries to 1 ng/uL, and then used Agilent 2100 to determine the insert sizes of the libraries. If the insert sizes met the expectation, we used qRT-PCR to determine the effective concentrations of the libraries. To ensure the quality of the library, the effective concentration of a library should be higher than 2 nM.
After the libraries were qualified, they were pooled and sequenced using an Illumina HiSeq 2500 platform according to the effective concentrations and the requirement of target data quantity. Paired-end sequencing with the 150 bp sequencing read length was performed. The sequencing was done by Novogene (Beijing, China) Bioinformatics Technology Co.Ltd.

Sequencing Data Analyses
The original image data files obtained from the Illumina HiSeq 2500 platform were transformed into raw sequenced reads by Consensus Assessment of Sequence and Variation (CASAVA pipeline v2.0, Illumina Inc., San Diego, CA, USA) base calling analyses. In order to ensure the quality of bioinformatic analyses, the raw reads were filtered to obtain clean reads. The steps for data processing were as follows: First, the adaptor sequences in the raw reads were trimmed; second, reads with more than 10% unknown nucleotides or with more than 50% low-quality bases (bases quality ≤20) were removed. The follow-up analyses were based on the obtained clean reads. Sequence alignment software HISAT 2.0.4 (Johns Hopkins University, Baltimore, MD, USA) was used to align the sequences of each sample with the reference genome (GRCg6a-galGal6-Genome-Assembly-NCBI.https://www.ncbi.nlm.nih.gov/ assembly/GCF_000002315.6) to obtain the position of the clean reads on the reference genome and the barcode sequence information unique to the sequenced samples.
The gene expression levels of each sample were analyzed by HTSeq (version 0.6.1) software. The Union model was chosen, and FPKM (Fragments Per Kilobase of exon per Million fragments mapped) values were calculated. Next, we used the DESeq R package (version 1.10.1) [15] software to standardize the gene expression of different samples and then calculated the p values according to the model and finally corrected the multiple hypothesis tests (FDR). After the correction, for pituitary glands and ovaries tissue, the Padj of 0.05 and log2 (Fold_change) with no limitations served as the threshold of significance for differential expression. Only one candidate gene was screened when Padj <0.05 gene was selected as the differential expression gene in the hypothalamus. Therefore, the adjusted threshold for differentially expressed genes in the hypothalamus was p-value < 0.005. The differentially expressed genes were analyzed by GOseq software (Illumina, San Diego, California, USA) [16], and the pathway analyses were done by KOBAS (KO Based Annotation System) [17].

Real-Time PCR Validation of Sequencing Results
Randomly selected from the transcriptome sequencing results, 15 differentially expressed genes related to HEP and LEP of Xinjiang Yili geese were used for fluorescence-based quantitative validation. Oligo 7.0 software (Molecular Biology Insights, CO, USA) was used to design primers (Table S1). SYBR GREEN reagent (TaKaRa) was used to amplify the target gene and internal reference gene (beta-actin) mRNAs on a ROCHE 480 quantitative PCR instrument (Eppendorf, Hamburg, Germany). The PCR reaction system (20 µL) included the 2 µL total RNA as the template, 1 µL upstream primer (10 µmol/L), 1 µL downstream primer (10 µmol/L), 10 µL 2× master mix, and 6 µL ddH2O. The reaction condition was as follows: 95 • C for 15 min, followed by 40 cycles of 95 • C for 10 s, 58 • C or 60 • C for 20 s, 72 • C for 20 s. Quantitative expression results were calculated according to the cross point (CP) values, and the relative expression levels were calculated according to the 2 −∆∆Ct method [18]. Primer sequences are given in Table S1.

Evaluation of RNA and Sequencing Data Quality
The Nanodrop (Implen, Weslake Village, USA), Qubit 2.0 (Life Technologies, California, USA), and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) were used to determine the purity, concentrations, and integrity of 24 RNA samples. Table S2 shows that the RNA quality of the samples was good, and the purity and integrity of the samples met the requirements of sequencing library construction.
The quality assessment of the output data of 24 samples is shown in Table S3. According to the  table, 579,563,136 and 596,933,010 reads were obtained by sequencing Xinjiang Yili geese with HEP and LEP. We obtained a total of 176.47 gigabases of data, with at least 5.98 gigabases for each sample. The Q20 and Q30 of each sample were at least 95.10% and 88.57%, respectively, and the GC contents (calculate the percentage of the total number of bases g and C in the total number of bases) were between 48.28% and 50.57%. These results indicated that the quality of transcriptome sequencing results met the needs of subsequent analyses.
The statistics of the sequencing read alignments of the 24 samples to the reference genome are shown in Table 1. According to the statistic results, 72.51% to 77.98% of the clean reads were mapped to the reference genome, among which 71.08% to 77.22% of the clean reads were uniquely mapped. The results showed that there were enough reads of each sample mapped to the reference genome, and the selected reference genome was suitable.

Identification of Differentially Expressed Genes
Gene expression displays temporal and spatial patterns. Both the external stimuli and internal environment can affect gene expression. Under two different conditions, genes with significant differences at expression levels were called differentially expressed genes (DEGs). Padj < 0.05 was used as the criteria for identifying differentially expressed genes in the pituitary and ovary, and p-value < 0.005 was used as the criteria for determining differentially expressed genes in the hypothalamus. In this study, the sequencing data of the HEP and LEP Xinjiang Yili geese were compared. A total of 135 differentially expressed genes were identified in the hypothalamus, including 79 up-regulated genes and 56 down-regulated genes. 56 differentially expressed genes were identified in the pituitary, including 25 up-regulated genes and 31 down-regulated genes. A total of 331 differentially expressed genes were identified in the ovary, including 312 up-regulated genes and 19 down-regulated genes (Figures 1-3, Table 2).

Identification of Differentially Expressed Genes
Gene expression displays temporal and spatial patterns. Both the external stimuli and internal environment can affect gene expression. Under two different conditions, genes with significant differences at expression levels were called differentially expressed genes (DEGs). Padj < 0.05 was used as the criteria for identifying differentially expressed genes in the pituitary and ovary, and p-value < 0.005 was used as the criteria for determining differentially expressed genes in the hypothalamus. In this study, the sequencing data of the HEP and LEP Xinjiang Yili geese were compared. A total of 135 differentially expressed genes were identified in the hypothalamus, including 79 up-regulated genes and 56 down-regulated genes. 56 differentially expressed genes were identified in the pituitary, including 25 up-regulated genes and 31 down-regulated genes. A total of 331 differentially expressed genes were identified in the ovary, including 312 up-regulated genes and 19 down-regulated genes. (Figures 1-3, Table 2) Figure 1. Volcano plot of differentially expressed genes in the hypothalamus. DEGs were filtered using a p-value < 0.005 as a threshold. Red spots represent up-regulated genes, and green spots indicate down-regulated genes. Blue spots represent genes that did not show obvious changes between the HEP and LEP samples.  Differentially expressed genes (DEGs) were filtered using Padj < 0.05 as a threshold. Red spots represent up-regulated genes, and green spots indicate down-regulated genes. Blue spots represent genes that did not show obvious changes between high-egg production (HEP) and low-egg production (LEP) samples.  Differentially expressed genes (DEGs) were filtered using Padj < 0.05 as a threshold. Red spots represent up-regulated genes, and green spots indicate down-regulated genes. Blue spots represent genes that did not show obvious changes between high-egg production (HEP) and low-egg production (LEP) samples.

GO and KEGG Pathway Enrichment Analyses of Differentially Expressed Genes
Blast2 GO software was used to compare the differentially expressed genes and annotate GO functions. The GO annotation includes the enrichment of molecular functions, biological processes, and cellular components [19] (Figures 4-6).
passive transmembrane transporter activity, and calcium ion binding.
Besides, the GO terms of steroid biosynthetic process, steroid hormone-mediated signaling pathway, reproduction, developmental process, and G-protein coupled receptor signaling pathway were found to be associated with the hypothalamic DEGs; the terms of developmental and G-protein coupled receptor signaling pathway and calcium ion binding were found to be associated with the pituitary DEGs; the terms of developmental process, reproduction, and reproduction process were found to be associated with the ovarian DEGs. The differentially expressed genes on the above terms may be involved in the regulation of the reproductive traits of Xinjiang Yili geese. From the above GO terms, 19 differentially expressed genes were selected, including 14 up-regulated genes and 5 down-regulated genes, which could be used as candidate genes associated with the egg production of Xinjiang Yili geese (Table 3).      The hypothalamic differentially expressed genes were mainly associated with the biological process terms of organophosphate biosynthesis, single-organism development, and development; the cell component term of BBSome and molybdopterin synthase complexes; the molecular functions of anion transmembrane transporter activity, organic anion transmembrane transporter activity, and nickel cation binding. The pituitary differentially expressed genes were mainly associated with the biological process terms of protein metabolism, protein oligomerization, and anion transport; the cellular component term of Nem1-Spo7 phosphatase complexes; the molecular functions of catalytic activity, nucleotide binding, nucleotide phosphate binding, etc. The ovarian differentially expressed genes were mainly associated with the biological process terms of regulation of cellular process, cellular response to stimulus, cell communication, signal transduction, and single organism signaling; the cell component terms of extracellular matrix, proteinaceous extracellular matrix, etc.; the molecular functions of ion channel activity, substrate-specific channel activity, passive transmembrane transporter activity, and calcium ion binding.
Besides, the GO terms of steroid biosynthetic process, steroid hormone-mediated signaling pathway, reproduction, developmental process, and G-protein coupled receptor signaling pathway were found to be associated with the hypothalamic DEGs; the terms of developmental and G-protein coupled receptor signaling pathway and calcium ion binding were found to be associated with the pituitary DEGs; the terms of developmental process, reproduction, and reproduction process were found to be associated with the ovarian DEGs. The differentially expressed genes on the above terms may be involved in the regulation of the reproductive traits of Xinjiang Yili geese. From the above GO terms, 19 differentially expressed genes were selected, including 14 up-regulated genes and 5 down-regulated genes, which could be used as candidate genes associated with the egg production of Xinjiang Yili geese (Table 3).  The pathway enrichment analysis was based on the KEGG database, and a hypergeometric test was used to identify pathways in which the differentially expressed genes were significantly enriched by comparing them with the reference genes of the whole genome. The enrichment analysis results are shown in Figures 7-9. The hypothalamic differentially expressed genes were mainly enriched in terms of metabolic pathways, phagosome, sphingolipid metabolism, tight junction, gap junction, focal adhesion, and ECM-receptor interaction signaling pathway; the pituitary differentially expressed genes were mainly enriched in terms of neuroactive ligand-receptor interaction, sphingolipid metabolism, regulation of actin cytoskeleton, metabolic pathways, gap junction, focal adhesion, ECM-receptor interaction, and calcium signaling pathway, etc. The ovarian differentially expressed genes were mainly enriched in terms of focal adhesion, ECM-receptor interaction, vascular smooth muscle contraction, phagosome, tight junction, gap junction, regulation of actin cytoskeleton, calcium signaling pathway, and MAPK signaling Pathway, etc. The differentially expressed genes of the hypothalamus, pituitary gland, and ovary were all enriched in terms of gap junction, focal adhesion, ECM-receptor interaction signaling pathway, and the differentially expressed genes of both the pituitary and ovary were enriched in the calcium signaling pathway, indicating that these pathways played an important role in regulating egg production of Xinjiang Yili geese. Eleven differentially expressed genes were selected from the associated pathways, all of which were up-regulated genes. These genes are candidate genes related to the egg production of Xinjiang Yili geese (Table 3). calcium signaling pathway, indicating that these pathways played an important role in regulating egg production of Xinjiang Yili geese. Eleven differentially expressed genes were selected from the associated pathways, all of which were up-regulated genes. These genes are candidate genes related to the egg production of Xinjiang Yili geese (Table 3).   The vertical axis represents the pathway name, the horizontal axis represents rich factor, the size of the point represents the number of differentially expressed genes in this pathway, and the color of the point corresponds to different Q-value ranges (the same below).

Fluorescence Quantitative Polymerase Chain Reaction
In order to validate the differentially expressed genes identified by the transcriptome sequencing, 15 genes were randomly selected and confirmed by qRT-PCR using beta-actin as the internal reference gene. The results showed that the expression trends of the 15 genes were consistent with those of transcriptome sequencing results. Therefore, the transcriptome sequencing results are reliable and can be further studied and analyzed ( Figure 10).

Fluorescence Quantitative Polymerase Chain Reaction
In order to validate the differentially expressed genes identified by the transcriptome sequencing, 15 genes were randomly selected and confirmed by qRT-PCR using beta-actin as the internal reference gene. The results showed that the expression trends of the 15 genes were consistent with those of transcriptome sequencing results. Therefore, the transcriptome sequencing results are reliable and can be further studied and analyzed ( Figure 10).
In order to validate the differentially expressed genes identified by the transcriptome sequencing, 15 genes were randomly selected and confirmed by qRT-PCR using beta-actin as the internal reference gene. The results showed that the expression trends of the 15 genes were consistent with those of transcriptome sequencing results. Therefore, the transcriptome sequencing results are reliable and can be further studied and analyzed ( Figure 10). Figure 10. Comparison of qRT-PCR and RNA-seq Results. Total RNA extracted from the hypothalami, pituitary glands, and ovaries tissues that were measured by qRT-PCR analysis; relative expression levels were calculated according to the 2 −∆∆Ct method using beta-actin as an internal reference gene.

Discussion
The egg production trait of the goose is a trait with low heritability. It is time-consuming to improve the egg production performance of geese by traditional breeding methods. The application of molecular biology technologies provides a new way to improve the egg production performance Figure 10. Comparison of qRT-PCR and RNA-seq Results. Total RNA extracted from the hypothalami, pituitary glands, and ovaries tissues that were measured by qRT-PCR analysis; relative expression levels were calculated according to the 2 −∆∆Ct method using beta-actin as an internal reference gene.

Discussion
The egg production trait of the goose is a trait with low heritability. It is time-consuming to improve the egg production performance of geese by traditional breeding methods. The application of molecular biology technologies provides a new way to improve the egg production performance of geese. RNA-Seq can reveal the expression, structural characteristics, and regulatory profiles of genes by analyzing the transcriptomes of different cells or organs [4], thus identifying candidate genes for the egg production traits. In this study, a high throughput sequencing technique was used to sequence the transcriptomes of the hypothalamic-pituitary-gonadal axis of Yili geese with high and low egg production. Of the differentially expressed genes, 135, 56, and 331 were identified in the hypothalamus, pituitary, and ovary, respectively. Some of the differentially expressed genes (DRD1, IGF1, ADCY3, SLC4A4, MRAP) are involved in signaling pathways related to reproductive performance, such as gap junction, adhesion, ECM-receptor interaction, calcium signaling pathway, MAPK signaling pathway, GnRH signaling pathway, oocyte meiosis signaling pathway, etc. Dopamine (DA) is a major neurotransmitter in the central nervous system, which plays an important regulatory role in the nervous, endocrine, and reproductive systems, and its biological effects are mediated through its receptors. Dopamine receptors belong to the G protein-coupled receptor family, which consists of 7 transmembrane regions. Five dopamine receptors were found until today: D1, D2, D3, D4, D5. D1 and D5 are D1-like receptors, which induced the increase of intracellular cAMPs after activation. D2, D3, and D4 are D2-like receptors, which decrease the level of intracellular cAMP [20]. Niu Shuling [21] found that subcutaneous injection of cAMP significantly improved the laying rate of Hailan hens, a commercial egg-laying hen. Young Ren [22] showed that dopamine could stimulate the secretion of vasoactive intestinal peptide and prolactin by stimulating DRD1 to regulate the reproductive performance of birds. Xu [23] reported that there were 7 SNPs in the CDS region of the DRD1 gene in chickens. The haplotypes of G123A and T198C were significantly correlated with the egg production traits and nesting ability of chickens, and the C1107T and G123A loci were correlated with the hatching rate of chickens. Wang [24] found that the mutation of the DRD1 gene, C681T, was significantly correlated with the pre-laying weight of ducks, and the mutation of A765T was correlated with egg production traits of ducks. Wang Cui [25] scanned the SNPs of the CDS sequence of the duck DRD1 gene. Five SNPs were detected: T189C, C507T, C681T, A765T, and A1044G. The variation loci were significantly correlated with the reproductive traits, egg production, and fat metabolism of ducks. Our study also indicates that DRD1 is up-regulated in the ovaries of HEP geese and may be involved in the regulation of the calcium signaling pathway. Neurotransmitters activate DRD1, stimulate the activation of adenylate cyclase, and increase intracellular cAMP levels, thus possibly improving the egg production rate.
Insulin-like growth factor-I (IGF-I) is a metabolic hormone secreted by granulosa cells. It can stimulate mitosis of follicular granulosa cells, promote the secretion of oxytocin and progesterone from the ovary, stimulate hormone synthesis and secretion from membrane cells, and Leydig cells of the testis, and synergistically work with FSH and estrogen to promote the synthesis of progesterone and estradiol. IGF-I is highly homologous to insulin and is one of the important proteins in the growth axis of animals. It plays an important role in regulating cell growth and differentiation [26]. The binding of IGF-I to its receptor can promote the proliferation of granulosa cells, maintain the function of aromatase, increase the synthesis of estrogen, and promote the further development of follicles [27]. The results of follicular culture in vitro showed that IGF-I could significantly increase the number and thickness of SYF granular and membrane cells in laying hens [28]. Zhu Wenqi [29] analyzed the effects of IGF-1 on the reproductive performance of Wenchang chickens. The results showed that IGF-1 was involved in calcium metabolism in vivo and had a significant effect on the numbers of eggs produced by Wenchang chicken at 300 and 400 days, and it was also related to the continuous egg production traits of Wenchang chickens. Yu Mingyue [30] found that there was a correlation between IGF-1 gene expression and egg production traits in Shiqiza Chicken. The results of our experiments indicated that IGF-1 was up-regulated in the ovaries of HEP geese, and its role in focal adhesion, oocyte meiosis, and progesterone-mediated oocyte maturation may regulate the reproductive performance of Xinjiang Yili geese by participating in oocyte capacitation and maturation.
Adenylate cyclase (ADCY) is a key signaling molecule downstream of G protein-coupled receptors, which are widely distributed in mammalian cells. This enzyme catalyzes adenosine triphosphate (ATP) to produce cyclic adenosine phosphate (cAMP) and release pyrophosphate. Among all the ADCY subtypes, 9 subtypes are closely related to ADCY3 [31]. According to Thomas [32], ADCY is localized on the ovarian cell membrane and is a marker enzyme for the separation of ovarian cell membranes. Lee [33] found this enzyme in the ovary of rats. Animal studies have shown that ADCY3 plays an important role in the regulation of the hypothalamus of rats [34]. In addition, ADCY3 is highly expressed in olfactory epithelial cells [35] and is primarily responsible for the detection of odors. It was found that the absence of ADCY3 could lead to the deficiency of odor-induced signals and the weakening of maternal behaviors in mice. In addition, ADCY3 is also associated with normal sperm cell development, sperm functions, and maintenance of male fertility [36]. The results of our experiments show that the ADCY3 gene is up-regulated in the ovaries of HEP geese and may participate in a variety of signaling pathways: Calcium signaling pathway, GnRH signaling pathway, oocyte meiosis, progesterone-mediated oocyte maturation, tight junction, gap junction, purine metabolism and so on.
Animal reproduction is an extremely complex physiological process that is constrained by exogenous and endogenous factors [37]. In the female reproductive tract, epithelial cells of the uterus and fallopian tube are capable of secreting high concentrations of HCO 3− , which is important for sperm capacitation and the subsequent fertilization processes [38]. The SLC4 family are very important HCO3-transmembrane transporters. SLC4A4, also known as electrosodium bicarbonate cotransporter isoform 1 (NBCe1), is a bioelectric Na + /HCO 3 − cotransporter [39] and is widely expressed in a variety of tissues, including the central nervous system [40], the male and female reproductive tracts (the epididymis, and vas deferens in the male reproductive tract, and the endometrial and fallopian tube epithelium in the female genital tract) [41,42]. Liu [41] detected mRNA expression of 5 NBCe1 (SLC4A4) isoforms, NBCe1-A to -E, in the reproductive tract of mice, among which NBCe1-D and NBCe1-E are novel isoforms. Studies by Gholami K [43] have shown that estradiol (E2) can up-regulate the expression of SLC4A4, which is translocated on the apical and basolateral membranes of the uterine and glandular epithelium under the influence of E2. The basolateral and apical NBCe1 is involved in mediating E2-induced uterine HCO 3 -secretion. By RT-PCR and other techniques, Wang et al. [44] found that nbce1 was expressed in the female genital tract, where nbce1 may participate in HCO 3 − Animals 2020, 10, 90 15 of 18 secretion of endometrial epithelial cells. In this study, the SLC4A4 gene was up-regulated in both the pituitary glands and ovaries of HEP Yili geese and may be involved in the regulation of the above pathways to affect the egg production performance of Xinjiang Yili geese. Melanocortin receptor 2 (MC2R) belongs to the class A of 7 transmembrane alpha-helix G protein receptors. It regulates the circadian rhythm of steroid hormone secretion and stress-induced changes through CA/cAMP/PKA signal transduction pathway [45]. MC2R functions are strictly dependent on its accessory proteins (MRAP), and the functional expression of MC2R can only be achieved in the presence of MRAP [46]. According to the recent research results, there are two types of MRAP: MRAP1 and MRAP2 [47,48]. Rouault [49] demonstrated that MRAP2 regulates the activity of Orexin Receptor 1 (OX1R) and interacts with various GPCRs to control food intake and energy expenditure. The data of our experiment showed that MRAP was down-regulated in the pituitary glands of HEP Yili geese. This is consistent with our observation that the Xinjiang Yili geese with LEP generally weighed heavier. Therefore, MRAP may play a role in egg production by regulating the goose appetite.

Conclusions
The results of our study showed that we obtained 1,176,496,146 valid reads and 176.47 gigabases of data from the transcriptomes of 24 samples after quality control. Differential gene analyses showed that 135 genes were differentially expressed in the hypothalami of the high-and low-egg yields Yili geese, including 79 up-regulated genes and 56 down-regulated genes; and 56 genes were differentially expressed in the pituitary glands of the high-and low-egg yields Yili geese, including 25 up-regulated genes and 31 down-regulated genes. A total of 331 differentially expressed genes were identified in the ovaries of the high-and low-egg yields Yili geese, of which 312 were up-regulated and 19 were down-regulated. According to the annotations of GO and KEGG analyses, 30 candidate genes related to the egg production of Xinjiang Yili geese were preliminarily selected. It was found that the adhesion spot, ECM receptor interaction, gap junction, and Ca 2+ signaling pathway may play an important role in regulating the egg production of Xinjiang Yili geese by the hypothalamic-pituitary-gonadal axis.