Transcriptome analysis of follicles reveals the importance of autophagy and hormones in regulating broodiness of Zhedong white goose

Broodiness, a maternal behavior and instinct for natural breeding in poultry, inhibits egg production and affects the poultry industry. Phenotypic and physiological factors influencing broodiness in poultry have been extensively studied, but the molecular regulation mechanism of broodiness remains unclear. Effective research strategies focusing on broodiness are hindered by limited understanding of goose developmental biology. Here we established the transcriptomes of goose follicles at egg-laying and broody stages by Illumina HiSeq platform and compared the sequenced transcriptomes of three types of follicles (small white, large white and small yellow). It was found that there were 92 up-regulated and 84 down-regulated transcription factors and 101 up-regulated and 51 down-regulated hormone-related genes. Many of these genes code for proteins involved in hormone response, follicular development, autophagy, and oxidation. Moreover, the contents of progesterone and estradiol in follicles were altered, and the autophagy levels of follicles were enhanced during the broody stage. These results suggest that hormone- and autophagy-signaling pathways are critical for controlling broodiness in the goose. We demonstrated that transcriptome analysis of egg-laying and broody Zhedong white goose follicles provided novel insights into broodiness in birds.

Scientific RepoRts | 6:36877 | DOI: 10.1038/srep36877 occurrence of broodiness 8,9 . Interestingly, in avian species, the granulosa cells only produce progesterone, and the theca folliculi produces testosterone or estradiol 13 . In female canaries and hens, progesterone dropped significantly at the onset of incubation 14,15 . Meanwhile, before incubation, the concentration of estradiol reduced and was maintained at low levels in the Chabo hen 16 . GnRH, synthesized and secreted by the hypothalamus, affects the ovarian activity controlling egg laying 17 .
Further studies showed that many other factors, including autophagy, homeostasis and ambient temperature, are associated with broodiness. An elevated ambient temperature accelerated the broody process in turkeys by increasing the level of PRL 10 . Our previous study indicated that the broodiness was associated with autophagy 7 . During the broody stage, the level of autophagy, identified by scanning electron microscope and the ratio of LC3-II to LC3-I, was elevated, and the homeostasis of follicles was impaired in goose follicles 7 . Also, many autophagy-related genes and miRNAs altered their expression in broody follicles tested by qPCR and miRNA transcriptome 7,18 . However, the mechanism of the occurrence of broodiness is not clear.
RNA-Seq, a next-generation sequencing technology, is a powerful tool to uncover transcriptional profiles for gene expression analysis to predict the possible molecular mechanism. RNA-Seq data are sufficient for gene model prediction, identification of novel transcripts, and measurement of transcript expression 19,20 . Moreover, RNA-Seq technology is much more sensitive and efficient for comparing gene expression profiles 20 . The successful sequencing of the goose genome opened a prelude to more detailed study of geese 21 . To date, RNA-Seq has been applied in geese, mainly based on the ovarian 22,23 , spleen 24 , exocrine 25 , hepatic 26 , follicles 27 , and mixed tissues 2 . Identified genes are responsible for reproductive biology, development and metabolism processes, and laying performance.
The genetic mechanism behind broody behavior has been widely studied for its potential effect on egg production. This mechanism of broodiness was considered to be controlled by a set of hormone-related genes of the ovarian tissues, including PRL and PRLR (prolactin receptor), and autophagy-related genes, including p53, Beclin1, Atg12, Atg9, and Caspase3 in previous studies. However, based on our previous study, we speculated that their different expressions in follicles might be the cause inducing broodiness 7 . To deeply explore the mechanism of broodiness occurring in geese, we sampled the SWF, LWF, and SYF during egg-laying and broody stages, and analyzed the whole transcriptomic data. It was revealed that certain genes involved in pathways for reproduction regulation, such as hormone biosynthesis, phagosome, pathway in cancer, focal adhesion, and regulation of actin cytoskeleton, were differentially regulated.

Results
Transcriptome sequencing and quality control. To obtain a global view of the Zhedong white geese follicular transcriptome and identify the genes involved in the regulation of broodiness, cDNA libraries from the different grade follicles of broody and egg-laying geese were constructed and sequenced. Totals of 30,417,008 to 100,673,900 raw reads were generated in each library, and the valid data were 29,846,346 to 99,189,858 reads, with the valid ratios being more than 96% for 18 libraries. The data qualities were good, with Q20 (base sequencing error probability < 1%) > 98% and Q30 (base sequencing error probability < 0.1%) > 82% in each library (Table 1).
To determine the precision of estimated FPKMs at the current sequencing depth, an FPKM saturation analysis was performed. With an increasing number of reads, the number of detected genes increases as well. Moreover, we found the sample correlations were more than 0.8 ( reasonable. Furthermore, the evenly distributed reads in every position of the genes indicate that the randomness of the breaking of these 18 samples was good. 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 goose genome. The mapping results show that 58.6% to 71.5% of the reads from each library were perfectly matched to the reference genome, whereas in unique mapped reads, 3.12% (371,746) to 3.45% (432,179), and in multi mapped reads, 3.12% (371,746) to 3.45% (432,179) were matched (Table S2). Gene expression of each sample showed a normal distribution ( Figure S1), suggesting that the gene expression of three biological replicates had the consistent trend.
Most importantly, 90% of the reads from each library were perfectly matched to the reference exons ( Fig. 1). Gene coverage was calculated as the absolute depth of read coverage across the genes from each of the broody and egg-laying goose follicles libraries. 0~1 of the absolute depth of the samples had coverage between 51~59%, and > 30 of the absolute depth of the samples had coverage between 12~25% ( Table 2). The results show that the sequencing quality complied with the requirements of further analyses.

Analysis of differential expression genes (DEGs).
We focused on discovering genes differentially expressed between egg-laying and broody hierarchical follicles, including SWF, LWF, and SYF. Among the 6 libraries, 20,868 genes were detected. Between the LSWF and BSWF libraries, 4,836 significantly differentially expressed genes (P < 0.05 and |log 2 ratio | ≥ 1) composed of 1,755 up-regulated and 3,081 down-regulated genes were identified (Fig. 2a,b). To determine the function of DEGs, we mapped them according to terms of the GO database. A total of 2,918 genes were categorized into the 3 main categories of GO classification, including biological processes, cellular components, and molecular functions ( Fig. 2c-e). Among them, the most important biological processes included "transcription, DNA-dependent", "regulation of transcription, DNA-dependent", protein transport, apoptosis, and multicellular organismal development. The important cellular components included the membrane, nucleus, and cytoplasm. For molecular function, genes were involved in ATP, protein, zinc ion, and DNA binding. For the LWF libraries, the differentially expressed genes contained 2,921 upregulated genes and 2,306 down-regulated genes in the broody follicles compared to those in the egg-laying follicles. Meanwhile, for the SYF libraries, the differentially expressed genes contained 2,670 upregulated genes and 2,169 down-regulated genes in the broody follicles compared to those in the egg-laying follicles. GO enrichment results of the two libraries show that the cellular components, molecular function, and the biological processes were the same with those of SWF libraries.
Expression analysis of hormone-related genes. Endogenous hormones, including PRL, P2, and E4, are important for the development of follicles and the regulation of broodiness in the goose. In this study, we identified 22, 39, and 40 upregulated hormone-related genes and 23, 14, and 14 downregulated hormone-related genes in SWF, LWF, and SYF of broody geese, compared to these of egg-laying geese ( Fig. 3a,b; Table S3). The differentially expressed genes were divided into GnRH-contained, progesterone-contained, and steroid-related groups.

Expression analysis of transcription factors.
Considering the functional importance of transcription factors, we identified a total of 634 transcription factors expressed in follicles. In the LSWF and BSWF libraries, 176 genes were significantly differentially expressed, including 92 up-regulated and 84 down-regulated (Fig. 4a,b; Table S7). In the LLWF and BLWF libraries, 216 genes were significantly differentially expressed, including 141 up-regulated and 75 down-regulated. In the LSYF and BSYF libraries, 188 genes were significantly differentially expressed, 126 up-regulated, and 62 down-regulated. Together, 58 transcription factors were upregulated and 31 transcription factors downregulated in broody follicles compared with those in egg-laying follicles.
Of the 58 upregulated transcription factors (Fig. 4c, Table S7)  Expression analysis of autophagy-related genes. Based on the GO database, we identified 33 autophagy-related genes, including BECN1, MAP1LC3B, ATGs, PIK3CB, and others (Fig. 4d, Table S8). The expression of PIK3CB, MAP1LC3B, and BECN1 was upregulated in the SYF of the broody goose, which was confirmed by qRT-PCR. These results suggest that the autophagy-related genes might play roles in regulating goose broodiness.

Analysis of the Endogenous hormones level.
Hormones, including P4 and E2, play regulatory role in goose broodiness. We found that the levels of E2 increased significantly in the SYF of the broody goose, but there was no difference observed in the SWF and LWF during egg-laying and broody stages (Fig. 5a). P4 levels were elevated in the LWF and SYF during the broody stage compared to those of the egg-laying stage (Fig. 5b). The results indicate that P4 and E2 were associated with the broodiness.

Analysis of autophagy levels.
Our previous study showed that autophagy was associated with broodiness in geese 7 . Combined with the differential expression results from this study, we postulated that the autophagy played a critical role in regulating broodiness in geese. We confirmed the expression of autophagy-related genes, including TP63 and BECN1, by qRT-PCR (Fig. 6a,b). We also detected the level of LC3 protein, a biomarker of autophagy, in the follicles of different sizes during both egg-laying and broody stages (Fig. 6c). Since LC3-II levels increased in the broody follicles, it suggests that the autophagy levels were raised at the broody stage in geese.

Discussion
Broody behavior of birds is significant for their reproduction, but it affects the egg production with the degeneration of follicles. During the broody stage, LYF, SYF, and LWF are significantly reduced in geese and rapidly increase after broodiness 6 . Some studies have shown that the degeneration of follicles is associated with autophagy, apoptosis and homeostasis imbalances, including the altered hormones and ROS level 7 . Furthermore, the transcriptomic results indicate that the development or degeneration of follicles requires the coordinated actions of abundant genes controlled at the transcriptional and posttranscriptional levels 27 . In addition, miRNAs also function in the regulation of follicular development to regulate the bird broodiness 18 . However, previous studies mainly focused on the gene expressions of the ovary or the mixed follicles at the egg-laying and broody stages; they did not extensively analyze the regulatory mechanism of goose broodiness in follicles of different grades.
In the current study, we provided new insight into goose follicular transcriptomes at egg-laying and broody stages by RNA-seq transcriptome analysis. We performed transcriptome profiling of follicles of different grades (SWF, LWF, and SYF) during egg-laying and broody stages in geese to identify genes that are differentially expressed. Many hormone-related genes, autophagy-and redox-related transcription factors were identified.
Also, the hormones, such as progesterone and estradiol function in regulating the development of follicles, ovulation, and sexual behaviors in birds 28 and the altered levels of reproductive hormones are associated with broodiness 7 . In this study, we demonstrated that progesterone and estradiol increased and the expressions of hormone-related genes were altered in the broody follicles. These results indicate that hormones played an important role in regulating broodiness in geese. The cyclins are involved in regulating the follicular size. CCNB2, as an important factor for cellular divisions, rose sharply to increase follicle sizes 29 . CPEB1 is involved in Cyclin B translation and meiotic resumption to regulate the development of oocytes 30 . In this study, we found that CCNB2, CCNA1, CPEB1, CDC27, and CDC16 decreased significantly in the broody follicles, suggesting that the development of follicles was inhibited during the broody stage.
Insulin and insulin-like growth factor 1 (IGF1) could stimulate granulosa cells producing estradiol 31 . Correspondingly, an elevated level of IGF1 and increased contents of estradiol were detected in the broody follicles in geese. MAPKs play important roles in follicle development 32 . The expression of MAPK14 was reduced during the broody stage in geese. Polo-like kinase 1 (Plk1), involved in meiotic arrest of oocytes 33 , also decreased in the broody follicles in geese. An increase in KSR1 expression leads to an increase in proliferative capacity of cells 34 . Conversely, the expression of KSR1 was downregulated in the broody follicles, resulting in the bad follicular development. Our results suggest that broodiness affected the development of follicles.    Autophagy is involved in both cell survival and cell death 35 . Our previous study showed that autophagy was associated with broodiness in geese 7 . In the present study, we detected an elevated level of autophagy in the broody follicles and many autophagy-related genes, such as BECN1, TP63, and ATGs, alerted their expression in the broody follicles. ATF4 regulates a set of autophagy genes, including Atg16l1, Map1lc3b, Atg12, Atg3, Beclin1, and Gabarapl2, implicated in the formation, elongation, and function of the autophagosome to control autophagy and respond to severe hypoxia 36 . MYC also functions in induction of autophagy 37 .
Increased TNFSF11 production is associated with increased oxidative stress in osteoblasts 38 . Akt/PKB activation blocks FoxO3 activation and autophagy, and this effect is not prevented by rapamycin. FoxO3, blocked by Akt/PKB, controls the expression of autophagy-related genes, including LC3 and Bnip3, to regulate autophagy 39 .  The overproduction of chemerin (CMKLR1) in mice reduced the generation of mitochondrial reactive oxygen species (ROS) and increased mitochondrial autophagy, as determined by the increased conversion of LC3-I to LC3-II and higher expression levels of Beclin1 and autophagy-related protein-5 and 7 40,41 . TRP63 and TRP73, the members of TRP53 family, can bind and regulate target genes through the same consensus site as TRP53 to activate a global autophagy program 42,43 . TFEB, having a much broader activity, controls a number of autophagy genes involved in multiple crucial steps of the autophagic pathway such as autophagosome formation, autophagosome lysosome fusion and lysosome-mediated degradation of the autophagosomal content. The regulation of TFEB activity, not affected by rapamycin, is similar to that of FOXO3a, and TFEB and FOXO3a act in parallel mechanisms to regulate autophagy 39,44 .
In this study, many homeobox domain genes, involved in the regulation of patterns of anatomical development, were differentially expressed in broody goose follicles. Newborn ovary homeobox (NOBOX), an oocyte-specific transcription factor, is preferentially expressed in the germ cells throughout folliculogenesis 45,46 . Female mice lacking NOBOX are infertile due to postnatal oocyte loss and a disrupted transition in follicular development from primordial to primary follicle 46 . In the broody goose, the NOBOX genes were downregulated in LWF and SYF compared to those in egg-laying goose. Combined with the sharply reduced numbers of LWF and SYF in geese during broodiness, we suggest that the altered expression of homeobox genes controlled the follicular development during broodiness.
In the present study, a comparative RNA-Seq transcriptomics approach led to the identification of hundreds of genes differentially expressed in goose follicles during egg-laying and broody stages. Many of them functioned predictably on the pathways of hormone response, follicular development, autophagy, and oxidation. Furthermore, the rise of broodiness in geese was accompanied by the altered levels of progesterone and estradiol and the elevated autophagy levels of follicles marked by LC3-II to LC3-I. Overall, this study provided novel insights into broodiness in birds and demonstrated that RNA-Seq is a powerful tool for examining the molecular mechanism in regulating broodiness.

Materials and Methods
Ethics statement. All experimental protocols employed in the present study that relate to animal experimentation were performed in accordance with Measures for the Administration of Affairs Concerning Experimental Animal of Zhejiang Province, China (Approved by the Zhejiang Provincial Government in 2009 and promulgated by Decree No. 263). All animal experiments of this study were approved by the Animal Care and Use Committee of the Zhejiang Agriculture and Forestry University (Lin'an, Zhejiang, China), in order to ensure compliance with international guidelines for animal welfare.
Animals, feeding, and follicle collection. A total of 6 Zhedong white geese, including 3 broody and Construction of mRNA library and sequencing. Total RNAs of goose follicles were extracted using Trizol reagent (Invitrogen, CA, USA) following the manufacturer's procedure. Their quantity and purity were analyzed by the Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RNA integrity > 7.0. Approximately 5-μ g total RNA of each sample was subjected to isolating Poly (A) mRNA with poly-T oligo attached magnetic beads (Invitrogen). Following purification, the mRNA is fragmented into small pieces using divalent cations under elevated temperature. Then the cleaved RNA fragments were constructed into the cDNA library according to the protocol for the Illumina RNA ligation based method (Illumina, San Diego, USA). We performed the single end sequencing on an Illumina Hiseq 2500 at the LC Sciences, China, following the vendor's recommended protocol.
Sequencing and primary analysis. A cDNA library was sequenced run with Illumina 2500 sequence platform. Using the Illumina paired-end RNA-seq approach, we sequenced the transcriptome of SWF, LWF, and SYF sampled during egg-laying and broody stages, respectively. Prior to assembly, the low quality reads (containing sequencing adaptors, sequencing primers, or nucleotide with q quality scored lower than 20) were removed. The raw sequence data were submitted to the NCBI Short Read Archive with accession numbers.

RNA-seq reads mapping.
We aligned the reads of 18 samples to the UCSC (http://genome.ucsc.edu/) goose reference genome using Tophat package, which initially removed a portion of the reads based on quality information accompanying each read and then mapped the reads to the reference genome of goose. Tophat allows multiple alignments per read (up to 20 by default) and a maximum of two mismatches when mapping the reads to the reference.
Transcript abundance estimation and differentially expressed testing. The aligned reads were processed by Cufflinks, which use the normalized RNA-seq fragment counts to measure the relative abundances of the transcripts. The unit of measurement is fragment per kilobase of exon per million fragments mapped (FPKM). The reference GTF annotation used in Cuffinks was downloaded from the UCSC database. Cufflink was used to de novo assemble the transcriptome. The second, Cuffmerge was used to co-merge all transcripts of samples to generate unique transcripts. The downloaded UCSC GTF file was passed to Cuffdiff along with the original alignment (SAM) files produced by Tophat. Cuffdiff re-estimates the abundance of the transcripts listed in the GTF file using alignments from the SAM file and concurrently tests for different expression. Only the comparisons with "p value" less than 0.05 and status marked as "OK" in the Cuffdiff output were regarded as showing differential expression.
The differentially expressed genes, including transcription factor, hormone related genes, and autophagyrelated genes, were analyzed using MeV software (http://www.tm4.org/mev.html), and confirmed partially by qRT-PCR methods. Gene Ontology (GO) of differential genes was carried out using GeneOntology database (http://www.geneontology.org/). KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis was analyzed by KEGG database.
Analysis of P4 and E2. P4 and E2 in the different grade follicles (SWF, LWF, and SYF of egg-laying and broody goose) were analyzed by using a goose progesterone assay kit (H089), and goose estradiol assay kit (H102), respectively (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). The samples were processed according to the previous report 7 . All the analyses were carried out according to the kit protocols in three biological repeats. Statistical analysis was carried out by Microsoft Excel. Variance was compared based on t-test and the significance was determined when P < 0.05.