Edited by: Marco Milanesi, University of Tuscia, Italy
Reviewed by: Genxi Zhang, Yangzhou University, China; Yan Wang, Sichuan Agricultural University, China
†These authors have contributed equally to this work
This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Ovarian follicular development is an extremely complex and precise process in which the hypothalamic-pituitary-ovarian (HPO) axis plays a crucial role. However, research on the regulatory factors of the HPO axis is sparse. In this study, transcriptomes of the tissues in the entire HPO axis at 15, 20, 30, and 68 w of age were analyzed. In total, 381, 622, and 1090 differentially expressed genes (DEGs) were found among the hypothalamus, pituitary, and ovary, respectively. In particular, the greatest number of DEGs (867) was identified from the comparison of ovary at 30 and 15 w, which might be related to ovarian development and function at high ovulation capacity. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses indicated that most of these DEGs in the significantly enriched biological process (BP) terms and pathways were primarily involved in tissue development and the regulation of reproductive hormone biosynthesis and secretion. The latter is highly related to the HPO axis. Therefore, a number of hub candidate genes strongly associated with the HPO axis in each tissue were filtered by analyzing the Protein-protein interaction (PPI) network and seven known reproductive hormone-associated key genes were obtained:
Egg production is an essential trait for laying hens and for breeding in poultry. The growth and development of chicken ovarian follicles has been found to be directly related to the egg-laying performance (
In laying hens, ovarian follicular development is an extremely complex and precise process involving pre-hierarchical follicle recruitment, selection, and dominance, pre-ovulatory follicles, and ovulation (
In recent years, the rapid advancement in molecular biotechnology has been utilized to carry out research into the molecular mechanisms. Hence, to explore the novel genes and regulatory mechanisms of the HPO axis during different stages of ovarian development in hens, the transcriptomes of whole HPO tissues (hypothalamus, pituitary, and ovary) were analyzed by RNA-seq during four stages: initial ovarian development, sexual maturation, the peak and late laying periods. This study will provide new insights into the function of the HPO axis at different developmental stages of the ovary during reproduction.
All animal experiments and animal care in this study were performed in accordance with the regulations for the administration of affairs concerning experimental animals (Revised Edition, 2017). The protocols were approved by the Henan Agricultural University Institutional Animal Care and Use Committee (Permit Number: 19-0068).
A total of 36 healthy Hy-line brown laying hens at 15, 20, 30, and 68 weeks of age were provided by the poultry germplasm resource farm of Henan Agricultural University and were housed in separate cages, in which each age phase consisted of 9 hens. Commercial feed and clean drinking water were supplied
Total RNA was isolated from 36 tissues samples (12 hypothalamus, 12 pituitary, and 12 ovary samples) using Trizol RNA extraction reagent (Invitrogen, Carlsbad, CA, United States), with each sample consisting of a homogenous mixture from three hens to reduce individual variation. The purity, concentration, and integrity of the RNA were measured and checked using the Nano Photometer spectrophotometer (IMPLEN, CA, United States), Qubit RNA Assay Kit in Qubit 2.0 Fluorometer (Life Technologies, CA, United States), and an RNA Nano 6000 assay kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, United States), respectively. All cDNA libraries were sequenced using an Illumina Hiseq 2500 platform and 125 bp paired-end reads were generated.
Raw data (raw reads) in the fastq format were first processed through in-house Perl scripts. Clean data (clean reads) were screened by removing reads containing adapters, reads containing ploy-N, and low-quality reads from the raw data. At the same time, Q20, Q30, and GC content of the clean data were calculated. All the downstream analyses were based on the clean data of high quality.
FPKMs (Fragments per kilo-base of exon per million fragments mapped) of gene expression levels in each sample were calculated using StringTie (v2.1.1) and the obtained read counts were normalized by TMM (
Hierarchical cluster analysis of the DEGs filtered from the three tissues at different weeks of age was performed using g-plots in the R package (
The DEGs associated with the regulation of the HPO axis on ovarian function were extracted according to the combinations of the GO and KEGG analyses. All protein interaction predictions for these DEGs were performed using the STRING database (
RNA expression was detected using a Primer Script RT reagent Kit (TaKaRa, Dalian, China), according to the manufacturer’s protocol. Real-time PCR was performed using iTaqTM Universal SYBR Green Supermix Kit (Bio-Rad Laboratories Inc., Waltham, MA, United States) using a Light Cycler 96 instrument (Roche, Basel, Switzerland) at 95°C for 3 min; 40 cycles of 95°C for 10 s, annealing at 60°C for 30 s, 72°C for 30 s, and 72°C for 1 min. The sequences of all primers are listed in
To evaluate the quality of the Illumina sequencing, a total of 3,635,530,764 raw reads were generated from the 36 RNA-Seq datasets. A total of 3,455,715,304 high-quality clean reads were obtained after screening the raw reads and were used for mapping and transcriptome assembly. The average Q20, Q30, GC contents, and error rates were 96.21, 90.86, 49.52, and 0.02%, respectively, which was high enough for further analysis (
Raw data, clean data, quality and GC content of 36 transcriptomes from the tissues of hypothalamus-pituitary-ovarian (HPO) axis at four ages of hens.
Sample name | Raw reads | Clean reads | Q20 (Clean reads) | Q30 (Clean reads) | GC Content (clean reads) | Error rate |
H15_1 | 107070970 | 103730874 | 94.47% | 89.22% | 49.01% | 0.03% |
H15_2 | 110926506 | 107375016 | 94.87% | 89.02% | 48.45% | 0.03% |
H15_3 | 108520724 | 104843110 | 94.60% | 88.60% | 49.07% | 0.03% |
H20_1 | 102591682 | 91978252 | 95.40% | 88.84% | 50.35% | 0.02% |
H20_2 | 93401116 | 89429338 | 97.01% | 92.38% | 50.04% | 0.02% |
H20_3 | 110509498 | 106746240 | 97.18% | 93.01% | 51.93% | 0.02% |
H30_1 | 99591894 | 96672562 | 96.45% | 91.56% | 48.00% | 0.02% |
H30_2 | 108078742 | 102427084 | 96.13% | 90.73% | 48.91% | 0.02% |
H30_3 | 112978466 | 110575642 | 96.82% | 92.37% | 48.24% | 0.02% |
H68_1 | 110166072 | 106727232 | 97.47% | 93.72% | 50.05% | 0.01% |
H68_2 | 90082030 | 86223358 | 97.03% | 92.45% | 49.94% | 0.02% |
H68_3 | 87295316 | 83762850 | 96.87% | 92.12% | 48.79% | 0.02% |
P15_1 | 102387630 | 98132488 | 96.81% | 92.06% | 46.33% | 0.02% |
P15_2 | 91929454 | 87835908 | 96.88% | 92.17% | 47.09% | 0.02% |
P15_3 | 102458208 | 97687892 | 96.92% | 92.19% | 48.88% | 0.02% |
P20_1 | 90529502 | 86091934 | 96.92% | 92.24% | 48.17% | 0.02% |
P20_2 | 92065624 | 88161182 | 97.06% | 92.57% | 46.17% | 0.02% |
P20_3 | 82165648 | 79923374 | 96.96% | 92.40% | 50.17% | 0.02% |
P30_1 | 94008258 | 89616974 | 96.92% | 92.33% | 50.56% | 0.02% |
P30_2 | 86961510 | 80863012 | 94.58% | 87.05% | 46.55% | 0.03% |
P30_3 | 87202180 | 80953636 | 94.31% | 86.50% | 48.45% | 0.03% |
P68_1 | 111188934 | 105174500 | 95.02% | 88.63% | 48.79% | 0.02% |
P68_2 | 110463626 | 106352586 | 97.09% | 93.53% | 44.79% | 0.01% |
P68_3 | 101323228 | 97473420 | 97.10% | 93.55% | 44.68% | 0.01% |
O15_1 | 93111000 | 86789686 | 96.15% | 90.54% | 49.53% | 0.02% |
O15_2 | 93903244 | 88062558 | 95.92% | 90.01% | 50.09% | 0.02% |
O15_3 | 83885806 | 78689920 | 95.97% | 90.19% | 49.53% | 0.02% |
O20_1 | 98291872 | 92073930 | 95.84% | 89.91% | 53.08% | 0.02% |
O20_2 | 89376294 | 83980312 | 95.87% | 89.95% | 51.09% | 0.02% |
O20_3 | 102757886 | 94330090 | 95.60% | 89.30% | 52.32% | 0.02% |
O30_1 | 95300758 | 88280870 | 95.36% | 88.81% | 52.24% | 0.02% |
O30_2 | 99835164 | 96702800 | 96.54% | 91.44% | 52.70% | 0.02% |
O30_3 | 143326170 | 131267236 | 95.36% | 89.29% | 54.16% | 0.02% |
O68_1 | 91523808 | 85509102 | 95.30% | 88.72% | 51.79% | 0.02% |
O68_2 | 131756526 | 127234318 | 97.37% | 93.42% | 51.42% | 0.01% |
O68_3 | 118565418 | 114036018 | 97.58% | 93.88% | 51.50% | 0.01% |
Total | 3635530764 | 3455715304 | ||||
Average | 96.21% | 90.96% | 49.52% | 0.02% |
To investigate the key genes of the HPO axis that were involved in follicular development and ovarian function, the differentially expressed genes (DEGs) of the hypothalamus, pituitary, and ovary at the four stages of development of laying hens were analyzed (
Number of DEGs in the three tissues from the HPO axis at four different ages. The name of the abscissa consists of letters and numerals, the capital letters correspond to the tissue (H: hypothalamus; P: pituitary; and O: ovary), while the numerals indicates the weeks of age (i.e., 15, 20, 30, and 68 w). This naming convention is the same in the following figures. The blue bars represent the number of significantly down-regulated DEGs while red bars represent significantly up-regulated DEGs.
Venn diagrams of the DEGs in the three tissues of HPO axis at four ages. The capital letters in the upper left correspond to the analysis of different tissues [
Hierarchical clustering heatmaps of all the DEGs with different combinations of ages in the three tissues of the HPO axis were produced (
Hierarchical clustering heatmap of DEGs in the three tissues of HPO axis at four ages [
Upregulated and downregulated DEGs in the six combinations of each tissue were used for GO and KEGG analysis. In the GO enrichment analysis, there were 112 significantly enriched biological process (BP) items in the hypothalamus, and 178 and 709 significant BP terms were obtained in the pituitary and ovary, respectively (
Top 20 BP terms of gene ontology of differentially expressed genes in hypothalamus, pituitary and ovary [
KEGG pathway analysis of the HPO axis tissues is presented in
The significant KEGG classification of differentially expressed genes in hypothalamus, pituitary and ovary (
To identify the candidate DEGs that were probably related to the regulation of reproductive hormone synthesis and secretion according to the top 20 enriched BP terms and significant KEGG pathways, PPI network analysis was conducted, and 43, 42, and 11 DEGs in the hypothalamus, pituitary, and ovary, respectively, were identified. In this network, only
The protein-protein interaction (PPI) network of the DEGs that are associated with HPO axis function according to the GO and KEGG analyses [
Six genes with high expression from the tissues of the HPO axis were selected to validate the accuracy of the sequencing results:
The validation of differentially expressed mRNAs between qRT-PCR and transcriptome analysis [
The present study identified 381, 622, and 1090 differentially expressed genes from the hypothalamus, pituitary, and ovary in multiple comparisons among the four developmental stages of hens. A larger number of DEGs were found in the ovary compared to the other tissues, and most of them were upregulated at 30 and 68 w and were downregulated at 15 and 20 w, indicating that these genes may be associated with ovarian function and ovulation process. Some researchers have reported that the ovary of hens grows slowly before sexual maturation (18–20 weeks post hatch), and then rapidly develops until reaching the largest size with numerous ovarian follicles during the peak period of egg production (28–32 w post hatch). This could explain why the greatest number of DEGs was found in the comparison between the15 and 30 w (
The significantly enriched BP terms from the GO analysis in the three tissues of the HPO axis elucidated that most of the DEGs were involved in tissue development and the regulation of hormone secretion. The BP terms of catecholamine uptake, regulation of GTPase activity, regulation of phospholipase activity, neurotransmitter reuptake, and regulation of the phospholipase C-activating G-protein coupled receptor (GPCR) signaling pathway in the hypothalamus are well known to be involved in the release and function of GnRH/GnIH, which indicated that GnRH/GnIH conducted its effects via the related biology process and stimulated/inhibited the secretion of FSH and LH in the pituitary gland and ultimately affecting the reproductive system (
In the analysis of the KEGG pathway, 7, 4, and 9 significantly enriched pathways were obtained in the hypothalamic, pituitary, and ovary, respectively (
A total of 36 genes strongly associated with the regulation of hormone synthesis and secretion were observed in the PPI network analysis of the significantly enriched GO and KEGG pathways (
In the pituitary, there are seven highly related genes involved in intracellular cholesterol biosynthesis; in particular, 24-dehydrocholesterol reductase (
Overall, ovarian follicle development, maturation, and its functions were mainly controlled by the HPO axis and age-related alterations, and these novel potential genes derived from the significantly enriched BP terms of the GO and KEGG pathways associated with the regulation of reproductive hormone secretion and biosynthesis might be used to further study the regulatory mechanisms of the HPO axis on the reproductive system.
In conclusion, our study presented RNA-seq analysis of the HPO axis at four ages of hens and revealed that the HPO axis could be highly associated with the regulatory mechanisms of follicular development and ovarian function through the control of sex hormone biosynthesis and secretion in the ovaries of hens. It not only provides novel candidate DEGs and potential pathways for further research, but also provides a better understanding of the HPO axis function in the reproduction system and ovulatory cycle of poultry.
The RNA-seq datasets generated and/or analyzed during the present study are available at NCBI project PRJNA634976 with accession number
The animal study was reviewed and approved by the Administration of affairs concerning experimental animals (revised edition 2017) and Henan Agricultural University Institutional Animal Care and Use Committee (Permit Number: 19-0068).
YT and XK acquired the funding for the study. JL and CL performed research, analyzed the data, and wrote the manuscript. QL and GL analyzed the data and was involved in the study design. HL and WL performed the statistical analysis. YT and XK were involved in the design of the study. YT conceived the study and involved in its design and coordination. All the authors read and approved the final manuscript.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The author would like to thank all the students of Henan Innovative Engineering Research Center of Poultry Germplasm Resource for their contribution in sample collection and analysis.
The Supplementary Material for this article can be found online at:
The list of top 20 significantly enriched BP terms of all the DEGs in hypothalamus (
The list of top 20 significantly enriched BP terms of all the DEGs in pituitary.
The list of top 20 significantly enriched BP terms of all the DEGs in ovary.
The list of the significant KEGG classification of differentially expressed genes.
The list of gene primers sequencing used for qRT-PCR validation.