Transcriptome Analysis of Ovarian Follicles Reveals Potential Pivotal Genes Associated With Increased and Decreased Rates of Chicken Egg Production

Abstract

Egg production is an important economic trait in the commercial poultry industry. Ovarian follicle development plays a pivotal role in regulation of laying hen performance and reproductive physiology. However, the key genes and signaling pathways involved in the various-stages of laying hen follicular development remain poorly understood. In this study, transcriptomes of ovarian follicles at three developmental stages, the large white follicle (LWF), small yellow follicle (SYF), and large yellow follicle (LYF), were comparatively analyzed in hens with high (HR) and low (LR) egg-laying rates by RNA-sequencing. Eighteen cDNA libraries were constructed and a total of 236, 544, and 386 unigenes were significantly differentially expressed in the LWF, SYF, and LYF follicles of HR and LR hens, respectively. Among them, 47 co-transcribed differentially expressed genes (DEGs) in LWF and SYF, 68 co-expressed DEGs in SYF and LYF, and 54 co-expressed DEGs in LWF and LYF were mined. Thirteen co-expressed DEGs were found in LWF, SYF, and LYF follicles. Eighteen candidate genes, including P2RX1, CAB39L, BLK, CSMD3, GPR65, ADRB2, CSMD1, PLPP4, ATF3, PRLL, STMN3, RORB, PIK3R1, PERP1, ACSBG1, MRTO4, CDKN1A, and EDA2R were identified to be potentially related to egg production. Furthermore, Kyoto Encyclopedia of Genes and Genomes analysis indicated neuroactive ligand-receptor interaction, cell adhesion molecules, peroxisome proliferator-activated receptor pathway, and cAMP signaling pathway might elicit an important role in formation of egg-laying traits by influencing ovarian follicle development. This study represents the first transcriptome analysis of various-sized follicles between HR and LR hens. These results provide useful molecular evidence for elucidating the genetic mechanism underlying ovarian follicle development associated with egg production in chicken.

Keywords: chicken; differentially expressed gene; egg production; ovarian follicle; transcriptome.

FIGURE 1

Basic information of sequencing data. (A) Pairwise Pearson’s correlation coefficients of sequencing data from six samples with three replicates; (B) Principal component analysis (PCA) of transcriptomic variation. Shapes denote different groups of the follicles sampled, colored dots indicate different samples; (C) annotation of gene numbers based on different databases; (D) Venn diagram analysis of distribution of annotated genes of different follicle groups of HR and LR hens. A1, LWF of HR hens; B1, LWF of LR hens; A2, SYF of HR hens; B2, SYF of LR hens; A3, LYF of HR hens; B3, LYF of LR hens.

FIGURE 2

Volcano Plot of differently expressed unigenes of follicles in HR and LR chickens. The two vertical lines of dashes present twice of the difference threshold, and horizontal line of dashes indicates p-value $\frac{1}{2}$ of 0.05. Red dots represent significantly up-regulated genes, blue dots indicate down-regulated genes, and gray dots imply non-differentially expressed genes (p-value adjusted for multiple testing < 0.05). (A) A1 vs. B1; (B) A2 vs. B2; (C) A3 vs. B3; A1, LWF of HR hens; B1, LWF of LR hens; A2, SYF of HR hens; B2, SYF of LR hens; A3, LYF of HR hens; B3, LYF of LR hens.

FIGURE 3

GO functional classification of DEGs in the prehierarchical follicles sized 3.5–5.5 mm in diameter of HR and LR hens. Y-axis: three major functional annotations of GO terms including biological process, cellular component, and molecular function. X-axis: the number of DEGs annotated in the group.

FIGURE 4

GO functional classification of DEGs in the prehierarchical follicles sized 6.0–8.0 mm in diameter of HR and LR hens. Y-axis: three major functional annotations of GO terms including biological process, cellular component, and molecular function. X-axis: the number of DEGs annotated in the group.

FIGURE 5

GO functional classification of DEGs in the preovulatory follicles sized 8.5–10.5 mm in diameter of HR and LR hens. Y-axis: three major functional annotations of GO terms including biological process, cellular component, and molecular function. X-axis: the number of DEGs annotated in the group.

FIGURE 6

Predicted protein to protein interaction network of the four co-expressed DEGs in LWF, SYF and LYF follicles of HR and LR hens. (A) Presumed regulatory network of GABRA1 in neuroactive ligand-receptor interaction pathway by using the STRING database; (B) Inferred interaction network of NDUFAB1 in oxidative phosphorylation pathway; (C) Deduced regulatory network of NCAM2 in cell adhesion molecules interaction pathway; (D) Predicted interaction network of LOC424014 in putative methyltransferase.

FIGURE 7

Evaluation and comparison of mRNA expression levels (Log2FC) of 22 candidate genes when analyzed by RNA Seq and quantitative RT-qPCR. (A) Relative expression levels of the DEG mRNA co-expressed in LWF and SYF; (B) Relative expression levels of the DEG mRNA co-expressed in SYF and LYF; (C) Relative expression levels of the DEG mRNA co-expressed in LWF and LYF; (D) Relative expression levels of the DEG mRNA co-expressed in LWF and SYF and LYF. The relative expression level of each gene was quantified by the 2^−ΔΔCT method and 18S rRNA was used as the internal control for normalization of the results.