Transcriptome Profiling Identifies Differentially Expressed Genes in Huoyan Goose Ovaries between the Laying Period and Ceased Period

The Huoyan goose is famous for its high egg-laying performance and is listed as a nationally protected domestic animal by the Chinese government. To elucidate the key regulatory genes involved in Huoyan goose egg laying, RNA from ovarian tissue during the ceased and laying periods was sequenced using the Illumina HiSeq 2000 sequencing platform. More than 12 million reads were produced in ceased and laying libraries that included 11,896,423 and 12,534,799 clean reads, respectively. More than 20% of the reads were matched to the reference genome, and 23% of the reads were matched to reference genes. Genes with a false discovery rate (FDR) ≤0.001 and log2ratio ≧1 or ≤−1 were characterized as differentially expressed, and 344 up-regulated and 344 down-regulated genes were classified into functional categories. Twelve genes that are mainly involved in pathways for reproduction regulation, such as steroid hormone biosynthesis, GnRH signaling pathways, oocyte meiosis, progesterone-mediated oocyte maturation, steroid biosynthesis, calcium signaling pathways, and G-protein coupled receptor signaling pathway were selected for validation by a quantitative real-time polymerase chain reaction (qRT-PCR) analysis, the qRT-PCR results are consistent with the general expression patterns of those genes from the Illumina sequencing. These data provide comprehensive gene expression information at the transcriptional level that might increase our understanding of the Huoyan goose's reproductive biology.


Introduction
The Huoyan goose is a native Chinese breed that tolerates coarse feed and cold temperatures, is extremely adaptable and weakly broody, and has high egg-laying rates and outstanding reproductive performance. It is considered a national treasure by the Chinese goose industry and was listed as one of the nationally protected domestic animals by the Chinese government in 2000 [1,2]. But recently, problems of variety degeneration, especially a decrease in the number of eggs laid, have become prominent and are hindering the conservation and utilization of this local goose breed. In order to elucidate the molecular mechanism that controls egg-laying in geese, the reproductive biology of the Huoyan goose need to be investigated.
The number of eggs laid by a bird is determined by the number of follicles destined for ovulation and the capacity of the oviduct to transform the ova into a hard-shelled egg. Beside environment and metabolism factors, ovarian follicle growth and development are also controlled by a myriad of endocrine, paracrine and autocrine factors, including gonadotropins, sex steroid hormones and growth factors. In poultry, the reproductive endocrine system and reproductive activity are strictly controlled by the hypothalamic-pituitary-gonadal axis [3]. The hypothalamus regulates reproduction by releasing neurohormones (gonadotropin-releasing hormones, GnRH) to the pituitary gland. The pituitary gland then synthesizes and releases gonadotropins (luteinizing hormone, LH; follicle-stimulating hormone, FSH), which in turn act on the gonads to stimulate gametogenesis (spermatogenesis, oogenesis) and sex steroid hormone secretion (gestagens, androgens, and estrogens). The ovary performs numerous roles critical for oocyte development and ovulation. Within the ovary, granulose cells are an important site for estrogen production for local use as well as providing endocrine signaling to other tissues [4]. The somatic cells of the ovary act to support the growth and development of the oocyte in preparation for the surge of luteinizing hormone (LH) which elicits a physiological response leading to ovulation. This response includes promoting meiosis, steroidogenesis, follicular development, cumulus cell expansion, luteinization, and progesterone production, ultimately promoting oocyte maturation [5,6]. The dynamic and highly regulated process of follicle development requires the coordinated actions of a great number of genes, which is orderly orchestrated at the transcriptional and posttranscriptional levels. Recently, genomic and transcriptomic studies have identified genes that are related to high egg production and expressed specifically in the hypothalamic-pituitary-gonadal axis (HPG). Shiue (2006) identified nine transcripts related to high egg production in the chicken hypothalamus/pituitary gland, thus providing a valuable resource for the identification of markers of high egg production in chickens [7]. Kang et al. (2009) identified 18 known and 8 unknown differentially expressed genes in the ovaries from the pre-laying to egg-laying stage [8]. Yen (2006) found that at least 19 genes were more highly expressed in the pituitary glands of laying geese compared with pre-laying geese [9]. Our team used the suppression subtractive hybridization (SSH) method to investigate gene expression profiles in the hypothalamus and pituitary gland of Huoyan geese during the laying period and ceased period [2,10]. However, these limited genetic data for geese do not provide global transcriptome information and are insufficient for elucidating the molecular mechanisms of productivity in the laying goose.
In recent years, the development of next-generation sequencing (NGS), which is also known as high-throughput or deep sequencing technology, has provided a powerful, highly reproducible and cost-efficient tool for transcriptomic research [11]. An advantage of the technology is that it may be used for gene discovery and expression profiling of organisms without a reference genome by the de novo assembly of generated short reads [12,13,14]. There are several advanced and alternative NGS platforms, such as Roche's 454 GS FLX, Illumina HiSeq 2000 and Applied Biosystems' SOLiD [15].
In this study, high-throughput sequencing technology was applied to generate comprehensive transcriptomic profiles of ovarian tissue from Huoyan geese during the ceased and laying periods using the Illumina HiSeq 2000 sequencing platform. A comparative analysis of transcriptomic data was performed, and a large number of genes were found to be differentially expressed between the ceased and laying periods. Our analysis found that certain genes involved in pathways for reproduction regulation, such as steroid hormone biosynthesis, GnRH signaling pathways, oocyte meiosis, progesterone-mediated oocyte maturation, steroid biosynthesis, calcium signaling pathways, G-protein coupled receptor signaling pathway, dopaminergic synapse, and MAPK signaling pathway are differentially regulated. Additionally, expression profiling of these differentially regulated genes were validated by quantitative real-time polymerase chain reaction (qRT-PCR).

Ethics statement
Experimental procedures were approved by the animal welfare committee of the College of Animal Science and Veterinary Medicine of Shenyang Agricultural University (No. 2011036) and performed in accordance with the Regulations for the Administration of Affairs Concerning Experimental Animals (China, 1988) and EU Directive 2010/63/EU for animal experiments. All of the surgery was performed according to recommendations proposed by the European Commission (1997), and all efforts were made to minimize the suffering of the animals.

Experimental animal and sample preparation
The Huoyan geese were selected from the Liaoning Huoyan Goose Stock Breeding Farm and raised according to the program used at the farm. During the experiment, geese were fed ad libitum with rice grain and supplemented with green grass or water plants whenever possible. Feed was provided during the daytime when the geese were released into an open area outside of the building. Nine peak-laying period geese were sampled at 12 months of age (average BW53.5¡0.6 kg). Nine ceased period geese were sampled at 15 months of age (average BW53.5¡0.4 kg). Geese were killed by exsanguination to obtain ovary samples, which comprised the whole ovary including the small and large yellow follicles. For those ovarian follicles at different development stages, the walls of ovarian follicles were isolated for RNA extraction. All of the samples were quickly dissected, frozen in liquid nitrogen, and stored at 280˚C. Numbers of eggs to sampling day were recorded for each goose in laying period group and ceased period group from the day of first egg to sampling day. Their average eggs to sampling day (mean ¡ S.E.M.) were 83.7¡2.3 and 115.3¡6.6, respectively.

RNA extraction, library preparation and sequencing
The total RNA was extracted using TRIzol reagent (Invitrogen, USA) following the manufacturer's protocol. The RNA quality was characterized initially on an agarose gel and NanoDrop 8000 spectrophotometer (NanoDrop, Thermo Scientific), and the integrity of the RNA samples was evaluated using a 2100 Bioanalyzer (Agilent Technologies, USA). To minimize the effect of transcriptome unevenness among individuals, the RNA samples from nine individuals were pooled within each group in equal amounts to generate one mixed sample per group. These mixed RNA samples were subsequently used to construct a complementary DNA (cDNA) library and perform Illumina deep sequencing.
The cDNA libraries of the ceased period and laying period were constructed according to the Illumina kit's manufacturer's instructions (Illumina, San Diego, CA). Briefly, a fragmentation buffer was mixed with magnetic beads and Oligo(dT) was used to isolate the messenger RNA (mRNA), and then the mRNA was fragmented into shorter fragments. The first strand of cDNA was then synthesized with random hexamer-primer using the mRNA fragments as templates. Buffer, deoxynucleotide triphosphates (dNTPs), Ribonuclease H (RNase H) and DNA polymerase I were added to synthesize the second strand. Double-stranded cDNAs were purified with the QiaQuick PCR extraction kit (Qiagen, Germany) and eluted with EB buffer for end repair and poly(A) addition. Finally, sequencing adapters were ligated to the 59 and 39 ends of the fragments, and the fragments were purified by agarose gel electrophoresis and enriched by PCR amplification to create a cDNA library. During the quality control (QC) steps, the Agilent 2100 Bioanalyzer and the ABI StepOnePlus Real-Time PCR System were used for quantification and qualification of the sample library. Finally, the library was sequenced using an Illumina HiSeq 2000 at the Beijing Genomics Institute (BGI) in Shenzhen, China. All of the technical steps were performed in duplicate.

Sequence data analysis and assembly
The original image data were translated into sequence data as raw data or raw reads via base calling. Because the algorithms used in de novo transcriptome construction of the short reads provided by the Illumina platform may be severely inhibited by sequencing errors, a stringent filtering process of the cDNA sequence was employed to select only clean reads. This process including the removal of reads with adaptors, reads with more than 5% unknown nucleotides and lowquality reads in which the percentage of low-quality bases (base quality #5) was higher than 50%. An analysis was then performed of the clean read data, which included the sequence quality assessment, saturation analysis of sequencing, experimental repeatability analysis, distribution of clean reads copy number, and alignment statistics of clean reads.

Reads mapping to reference genome and genes
Because of a lack of genomic resources for geese, an essential dataset containing reference genes expressed in the chicken genome was prepared to identify the genes corresponding to the reads in each library. The chicken reference genome and gene sets were downloaded from the official website of the National Center of Genome Research (ftp://ftp.ncbi.nih.gov/genomes/Gallus_gallus, ftp://ftp.ncbi. nih.gov/genomes/Gallus_gallus/Assembled_chromosomes/). All of the clean reads were mapped to the chicken genome using the Short Oligonucleotide Analysis Package alignment tool SOAPaligner/SOAP2 [16], which does not allow mismatches of more than two nucleotides.

Identification of differentially expressed genes
To identify the differentially expressed genes between the ceased period and laying period, expression levels were calculated using the RPKM (Reads Per Kb per Million mapped reads) method as described by Mortazavi et al. [17]. Differentially expressed genes (DEGs) were identified using a rigorous algorithm developed by BGI based on the method described by Audic et al. [18]. The false discovery rate (FDR) was used to determine the threshold P-Values in the multiple tests and analyses through manipulation of the FDR values. We used P#0.01, FDR #0.001 and |log 2 ratio| §1 as the threshold to judge the significance of gene expression differences [19].
For gene ontology (GO) annotations and GO functional classifications of DEGs, the Blast2GO program was used [20]. A gene ontology enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed to identify which DEGs were significantly enriched in GO terms and metabolic pathways at the Bonferroni-corrected P-value of #0.05 compared with the whole-transcriptome background using the formula described by Ting Shi et al. [21].

Quantitative Real-time PCR (qRT-PCR) Validation
To validate the reliability of the Illumina analysis, twelve putative genes involved in the regulation of reproduction and reproductive processes and sex steroid hormone biosynthesis were selected and detected by qRT-PCR. The qRT-PCR primers were designed using Primer 3.0 (http://frodo.wi.mit.edu/primer3), and those used for qRT-PCR analysis are listed in Table S1. Tissues from every three birds were pooled in each group, and then there were three sample pools in each group. The total RNA was extracted as described for the transcriptome library preparation and sequencing.
Two micrograms of total RNA was reverse-transcribed using the PrimerScript RT Reagent Kit (TaKaRa, Dalian, China), and the RT-PCR analysis was conducted on the Bio-Rad iQ5 Real-time PCR Detection System (BIO-RAD, California, USA). Each 25 ml reaction volume contained 1 ml 10 mM (each) forward and reverse primers, 12.5 ml 26 SYBR Premix Ex Taq II (Takara, Dalian, China), and 2 ml cDNA products, and the final volume was adjusted using PCR-water. The following PCR program was used for amplification: 15 min at 95˚C, 40 cycles of denaturation at 95˚C for 10 s and annealing and extension at 60˚C for 30 s. Beta-actin (ACTB) was selected as an internal reference gene, and the expression level of ACTB was used to normalize the qRT-PCR results for each gene. Negative controls without the cDNA template were included in this experiment. The standard curve testing was performed using a series of 10-fold diluted samples for the different genes. The slopes of standard curves and PCR efficiency of these genes were calculated to confirm if the qRT-PCR data were precise and reliable. The melting curves were analyzed to make sure that a single PCR product was amplified for each pair of primers, and the expression levels of these genes were verified from three independent biological replicates along with the internal reference gene. Fold changes were calculated using the 2 2DDCt method [22]. All of the data were expressed as the mean ¡ standard deviation, and a statistical analysis using Student's t tests was performed with SPSS 16.0 for Windows (SPSS Inc.) to evaluate whether the means were significantly different (P,0.05) or highly significantly different (P,0.01).

Accession code
The Illumina HiSeq 2000 sequencing data for the goose ovary transcriptome have been deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra) with the accession number SRR1119186.

Illumina transcriptome sequencing and reads assembly
To obtain a global view of the goose ovary transcriptome and identify the genes involved in the regulation of egg-laying, cDNA libraries from the ovarian tissues of ceased period geese and laying period tissues were constructed and sequenced using the Illumina HiSeq 2000 sequencing platform. More than 11.9 million and 12.6 million raw tags were generated in each library. Approximately 99% of the raw reads passed the filter, with 11,896,423 and 12,534,799 clean reads acquired in the two libraries ( Figure 1). The total produced base pairs were 582,924,727 and 614,205,151 in each library ( Table 1).

Assessment of cDNA Libraries
To confirm that the number of detected genes increases proportionally to the sequencing amount, a saturation analysis was performed. A sequence saturation analysis is used to measure the sequencing data of a sample. With an increasing number of reads, the number of detected genes increases as well. However, when the number of reads reaches a certain value, the growth rate of the detected genes flattens and saturates. Figure 2 shows a trend of saturation where the number of detected genes almost ceases to increase when the number of reads reaches approximately 5 million.  During preparation of the cDNA sequencing libraries, the mRNA was first fragmented into short segments by chemical methods and then sequenced. If the randomness of breaking is poor, the read preference from a specific gene region will directly affect the subsequent analysis. We used the distribution of read location on the genes to evaluate the randomness of breaking. Because genes have different lengths, the read location on a gene was standardized to a relative position (which was calculated as the ratio between the reads location on the gene and gene length), and then the number of reads in each position was counted [23]. In this study, the evenly distributed reads in every position of the genes indicated that the randomness of breaking of these two samples is good (Figure 3).

Reads mapping to the reference genome dataset
To identify the genes corresponding to clean reads in each library, the clean reads were mapped to the reference genes expressed in the chick genome. The mapping results showed that 3.19% (379,968) and 3.32% (416,772) of the reads from each library were perfectly matched to the reference genome, whereas approximately 3.12% (371,746) and 3.45% (432,179) of the reads from each library were perfectly matched to the reference genes. The percentages of unique reads matched to the genome were 17.81% (2,118,775)   Gene coverage was calculated as the percentage of a gene covered by reads from each of the ceased period and laying period libraries. Figure 4 shows the distribution gene coverage of the two samples. Approximately 1% of the total genes of the ceased and laying period samples had coverage between 90-100% (105 and 129 genes, respectively), and 3% of the total genes of the ceased and laying period samples had coverage between 80-90% (291 and 302, respectively) respectively.

Analysis of gene differential expression
To discover the differentially expressed genes between the ceased period and laying period, the gene expression levels were calculated using the RPKM method (Reads per Kb per Million mapped reads). Between the two libraries, 12120 entities were detected, and according to the applied criteria (FDR #0.001 and |log 2 ratio | §1), 688 significantly differentially expressed genes composed of 344 up-regulated and 344 down-regulated genes were identified (Table S2, Figure 5), which included genes related to reproduction and reproductive processes, sex steroid hormone biosynthesis, and hormone secretion, such as D (2)
A KEGG pathway analysis allowed us to determine functional classes of the differentially expressed genes during different egg-laying stages. In this study, a total of 585 genes mapped onto 214 pathways were selected following this process. A summary of the genes involved in these pathways has been included in Table S5. The largest category was metabolic pathways (Ko01100, 12.14%), which had 71 annotated genes. Notably, a specific enrichment of genes was observed for pathways involved in reproduction regulation, such as steroid hormone biosynthesis, GnRH signaling pathways, oocyte meiosis, progesterone-mediated oocyte maturation, steroid biosynthesis, calcium signaling pathways, G-protein coupled receptor signaling pathway, dopaminergic synapse, and MAPK signaling pathway (Table 2). Specifically, several genes are involved in multiple pathways, such as ADCY7 and ADCY2 (progesterone-mediated oocyte maturation, GnRH signaling pathways, oocyte meiosis, calcium signaling pathways, G-protein coupled receptor signaling pathway), PRKX (progesterone-mediated oocyte maturation, GnRH signaling pathways, oocyte meiosis, calcium signaling pathways, dopaminergic synapse, MAPK signaling pathway), PRKCB (GnRH signaling pathways, calcium signaling pathways, dopaminergic synapse, MAPK signaling pathway), PGR (progesterone-mediated oocyte maturation, oocyte meiosis), and LHCGR (calcium signaling pathways, G-protein coupled receptor signaling pathway).

Quantitative real-time PCR validation
As shown in Figure 6 and Table 3, the general expression patterns of twelve genes from the Illumina sequencing are consistent with the qRT-PCR results, which further support the reliability of the Illumina sequencing data. The discrepancies with respect to ratio may be attributed to the essentially different algorithms and sensitivities between the two techniques.

Discussion
During formation of the egg in the bird egg-laying cycle, the ovary undergoes dynamic hormonal, biochemical and cellular changes. The dynamic and highly regulated process of follicle development requires the coordinated actions of a great number of genes orchestrated at the transcriptional and posttranscriptional levels [24]. A number of researchers have proven that transcriptome sequencing of RNA provides a rapid and cost-effective method of generating wholetranscriptome sequences from numerous types of biological samples and identifying differentially expressed genes. In the present study, the Illumina sequencing method was used to analyze the transcriptome of goose ovarian tissue at two stages of the egg-laying cycle and identified 344 up-regulated and 344 down-regulated genes in the laying period compared with ceased period. Some of those involved in steroid hormone biosynthesis, GnRH signaling pathways, oocyte meiosis, progesterone-mediated oocyte maturation, dopaminergic synapse, steroid biosynthesis, calcium signaling pathways, G-protein coupled receptor signaling pathway, and MAPK signaling pathway. Poultry reproduction is tightly regulated by sex steroid hormones. Sex steroid hormones, such as progesterone and estradiol, are involved in the regulation of ovulation, gonadal differentiation, and sexual and nesting behaviors in birds through interactions with their intracellular receptors [25]. The hormone progesterone exerts a regulatory function on follicular development. The follicular and ovulation stages are dependent on gonadotropin release. In poultry production, progesterone is secreted by follicle granulose cells, and the target sites  Figure 6. Validation of the gene expression profile by qRT-PCR. Total RNA extracted from the ovarian tissues that was measured by RT-PCR analysis; relative expression levels were calculated according to the 2 2DDCT method using beta-actin as an internal reference gene and the ceased period group as the calibrator (relative expression 51). Error bar shows the standard deviation. Gene symbols indicating different genes refer to Table 2. *P,0.05, **P,0.01.
doi:10.1371/journal.pone.0113211.g006 Table 3. Evaluation of the expression profile variation between RNA-Seq and qRT-PCR for selected genes. Transcriptome of Different Egg-Laying Stages Goose Ovaries of progesterone action are the oviduct and hypothalamic-pituitary-ovarian axis. The secretion of progesterone is tightly correlated with the growth and status of the ovary. In the ovulation process, progesterone acts as a positive feedback signal to trigger the GnRH/LH pulse surge prior to ovulation [26]. Similarly to all steroid hormones, progesterone induces its biological effects by binding to an intracellular protein, which is the progesterone receptor (PGR) in this case. The PGR is localized in target tissues such as the brain, ovary and oviduct [27,28,29] , and its expression is regulated by several hormones, such as estradiol, progesterone and gonadotropins [30]. The expression of PGR is up-regulated after FSH treatment and down-regulated after LH treatment in the interstitial cells of chick ovarian medulla [31]. In the ovaries of laying hens after estrogen treatment, there is an increased number of PGR-positive stromal, thecal and granulosa cells [32]. Furthermore, a significant decrease in PGR concentration is observed in all of the cell types of chick oviducts after injection of progesterone [33]. In the present study, the expression of PGR is up-regulated in the laying period, which is consistent with progesterone acting via its receptor to affect the regulation of follicular maturation, ovulation and oviposition in poultry. Prolactin (PRL) is a single-chain peptide hormone that is primarily synthesized by the anterior pituitary and belongs to the prolactin/growth hormone (GH) family [34]. In birds, PRL has been reported to be involved in the regulation of many important physiological processes, including egg laying, induction and maintenance of incubation behavior, osmoregulation, immune-modulation, and gonadal development and functions [35,36]. Notably, PRL plays an inhibitory role in avian reproductive activities at all levels of the hypothalamic-pituitary-gonadal axis via inhibition of gonadotropin secretion, stimulation of incubation behavior, and development of atresia in ovarian follicles [37,38]. The role of PRL is to specifically combine with prolactin receptor (PRLR) in the membrane of effector cells. PRLR belongs to the class I cytokine receptor superfamily [39]. In adult chickens, turkeys and geese, the PRLR is widely distributed in various tissues [40,41], and several studies have indicated that PRLR may play an important role in ovarian development and the egg-laying process in geese [39,42,43]. The highest level of expression of PRLR in the ovaries of Eastern Zhejiang White Geese was found at the incubation stage; a lower level of expression occurred during the egg-laying stage, and the lowest level of expression was found in the out-of-lay stage [44]. Our transcriptome analysis showed that the expression of PRLR in the laying period was significantly decreased compared to that of the ceased period. Therefore, the down-regulation of PRLR might be reducing the inhibitory effect of PRL on reproductive activities.

Gene
The production of steroid hormones by ovarian cells is essential for follicular recruitment, oocyte maturation, ovulation, and hypothalamo-pituitary-gonadal (HPG) axis regulation. Steroidogenesis is therefore integral to successful reproduction. Sex steroid hormones are derived from cholesterol through a series of enzyme-catalyzed reactions. First, cholesterol is transported from the outer to the inner mitochondrial membrane by the steroid acute regulatory protein (STAR) [45,46,47]. At the cytochrome P450 cholesterol side-chain cleavage (P450scc) enzyme site (encoded by CYP11A1), CYP11A1 cleaves the side chain of cholesterol and converts it to pregnenolone [48]. Pregnenolone is further catalyzed by other steroidogenic enzymes, which HSD3B2 catalyze the conversion of pregnenolone to progesterone, 17-hydroxypregnenolone to 17-hydroxyprogesterone (17OHP), dehydroepiandrosterone (DHEA) to androstenedione and androstenediol to testosterone. CYP17A1 is responsible for the conversion of pregnenolone to 17-hydroxypregnenolone, progesterone to 17-hydroxyprogesterone (17OHP), 17OHP to androstenedione and 17-hydroxypregnenolone to DHEA. HSD17B1 catalyzes the conversion of DHEA to androstenedione and estrone to estradiol, HSD17B3 catalyzes the conversion of androstenedione into testosterone, and CYP19A1 catalyzes the conversion of testosterone to estradiol and androstenedione to estrone [49]. Among these enzymes, STAR and CYP11A1 are the first and rate-limiting steps in the steroidogenic pathway. In avian species, the ovary is organized into a stroma with primordial follicles (,1 mm in diameter), small white follicles (1 to 4 mm), yellowish follicles (4 to 8 mm), and a hierarchy of yellow follicles (F n to F 1 ; 9 to 35 mm) that ovulate on successive days. The process of growth and maturation of the ovarian follicle is associated with the differentiation of granulosa cells, which occurs before selection of the yellowish follicle into the preovulatory hierarchy [50,51]. During the transition of the prehierarchical follicles to a preovulatory hierarchy, the cells of the granulosa layer, which were initially stimulated by FSH and subsequently by LH, begin to express the STAR protein and CYP11A1 and produce progesterone (P 4 ), which is catalyzed by HSD3Bs [46,52]. Numerous investigations have suggested that the elevated P 4 production in the granulosa layers of the F 3 RF 1 follicles is predominantly associated with an increasing expression of STAR, CYP11A1 and HSD3B mRNAs [26,47,53,54] and enhanced HSD3B activity [55]. CYP17A1, which has both 17 alpha-hydroxylase and 17, 20-lyase activities, is a key regulatory enzyme in the steroidogenic pathway. In the gonads, CYP17 deficiency prevents gonadal sex steroid production and leads to female external genitalia development [56]. CYP19A1 and HSD17B1 are responsible for estrogen biosynthesis. After prenatal and neonatal flutamide exposure in porcine ovarian follicles, the elevated level of 17-estradiol is associated with overexpression of the follicle stimulating hormone receptor (FSHR), CYP19A1, and CTNNB1 genes [57]. However, despite the role of the key enzymes of the steroid synthesis pathway, steroid hormones have a negative feedback effect on steroidogenic enzyme expression and regulate their own production in ovarian tissue [58,59]. Notably, the genes encoding steroidogenic enzymes in our study were down-regulated in the laying period compared with the ceased period, which may be related to the negative feedback regulation of sex hormones on steroidogenic enzyme expression during the egg-laying period.
In addition to the above-mentioned genes, several G-protein-coupled receptors, including DRD2, VIPR2, NPY1R, and MEL1C, were identified. These proteins play key roles in the regulation of steroid hormone biosynthesis and reproductive functions. Dopamine belongs to a group of neurotransmitters called catecholamines and exerts its action by binding to specific membrane receptors.
Dopamine receptors are members of the G-protein-coupled receptors with seven transmembrane domains [60]. In avian species, two distinct dopamine receptor subtypes, DRD1 and DRD2, have been reported. Dopamine can inhibit PRL secretion via DRD2 at the pituitary level by operating through vasoactive intestinal peptide (VIP) [61]. The DRD2 gene is widely distributed in the hypothalamic-pituitary-ovarian axis of chickens, turkeys and geese, and it plays an important role in the regulation of hormone secretion and reproductive behaviors [62,63,64]. The elevation of DRD2 expression during the Huoyan goose egglaying period might be capable of down-regulating the expression of PRL. VIP belongs to a family of regulatory peptides, and in birds, a major role of VIP involves regulating the synthesis and secretion of the anterior pituitary hormone prolactin (PRL). VIP can increase the transcription of the gene encoding PRL, enhance PRL mRNA stability, and elevate hormone secretion [65]. The actions of VIP are mediated through an interaction with two receptor types (VIPR1 and VIPR2) that belong to the class II subfamily of the 7-transmembrane G-proteincoupled receptor superfamily. VIP receptors are present on the surface membranes of the anterior pituitary, hypothalamus, small intestine and granulosa cells [66]. The VIPR1 gene in the pituitary is involved in the regulation of reproductive activity in avian species [67]. The expression of VIP receptor mRNA in the pituitary is markedly reduced in nonphotostimulated and photorefractory turkey hens compared with laying and incubating hens [68]. In the ovaries of Huoyan geese during the egg-laying period (April-August), the rise of VIPR2 expression might be a result of the long photoperiod. Neuropeptide Y (NPY) is known to participate in the regulation of several physiological functions, such as gonadotropin release, sexual behavior, food intake, energy metabolism, and stress responses. In birds, NPY has multiple functions related to the HPG system and GnRH secretion [69,70]. The great variety of physiological functions of NPY is mediated through G-protein-coupled receptors (NPY1R, NPY5R) [71,72,73]. Melatonin plays an important role in the regulation of animal reproduction via binding to G-protein-coupled receptors. Three melatonin receptor subtypes (Mel-1a, Mel-1b and Mel-1c) have been shown to be differentially expressed in chicken ovarian granulosa cells [74]. In the gonads of birds, melatonin and its receptors play an important role in regulating the viability and maturation of ovarian granular cells, gonadal steroidogenesis, folliculogenesis and maturation of oocytes [75]. Changes in the expression levels of the melatonin receptors and estrogen concentrations in the hierarchical follicles of the Sichuan white goose are consistent with the changes that occur in the Huoyan goose [76]. In our study, the increase of Mel-1C expression during the egg-laying period suggests that melatonin may regulate ovarian function through its receptor.
In both avian and mammalian species, the process of ovarian follicle maturation is coupled with a functional differentiation of the granulosa cell layer. granulosa cells from mature, preovulatory follicles and granulosa-luteal cells produce appreciable amounts of steroid hormones [77]. granulosa cell differentiation and steroidogenesis during follicle maturation are mediated by the coordinately timed expression of endocrine, paracrine, and autocrine factors and the interaction of consequent signaling pathways [78]. For instance, in ovarian steroidogenic cells, gonadotropins bind to specific membrane G-protein coupled receptors (GPCRs) and lead to the activation of multiple signal transduction pathways, including the adenylate cyclase-/cAMP-dependent protein kinase A (PKA) signaling pathway and calcium-/calmodulindependent pathways [79]. Studies indicate that activation of the mitogen-activated protein kinase (MAPK) and protein kinase C (PKC) signaling pathways modulates steroid biosynthesis during follicle development [80]. In our study, list of differentially expressed genes are involved in pathways calcium signaling pathways, G-protein coupled receptor signaling pathway, and MAPK signaling pathway. Like PRKCB, ADCY7, PRKX, MAP2K3, and TGFBR2, their roles in follicle development and productivity of goose need to be further investigated.

Conclusions
In conclusion, this study presents the transcriptomic profiles of ovarian tissue from Huoyan geese during the ceased and laying periods using Illumina RNA-Seq technology. Differentially expressed genes which involved in pathways for reproduction regulation, such as steroid hormone biosynthesis, GnRH signaling pathways, oocyte meiosis, progesterone-mediated oocyte maturation, steroid biosynthesis, calcium signaling pathways, G-protein coupled receptor signaling pathway, dopaminergic synapse, and MAPK signaling pathway were identified, furthermore, the expression profiles were validated by qRT-PCR. These data provide comprehensive gene expression information at the transcriptional level that might increase our understanding of the Huoyan goose's reproductive biology.