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, egg production, differentially expressed gene, ovarian follicle, transcriptome Introduction Egg-laying production efficiency is an important economic trait in the world-wide poultry industry. Efficiency is determined by the number of ovarian follicles destined for ovulation or atresia and the capacity of the oviduct to transform the ova into a hard-shelled egg ([37]Etches et al., 1990; [38]Luan et al., 2014). Follicle development is an intricate and highly coordinated cellular process, which is characterized by a well-organized follicular hierarchy in high production egg-laying hens. In chickens with low egg-laying rates, such as the broiler breeder hen, follicular development is not well-organized, which results in reduced reproductivity ([39]Johnson, 2012). Growing follicles are mainly classified into primordial follicles (smaller than 0.08 mm in diameter), primary follicles (from 0.08 mm up to 1 mm in diameter), undifferentiated prehierarchical follicles (from 1 mm up to 8 mm in diameter), and preovulatory/hierarchical (larger than 9 mm in diameter) follicles ([40]Lovell et al., 2003; [41]Onagbesan et al., 2009). During growth and development of the hen ovary, a majority of the follicles undergo atresia, resulting in a limited number of follicles further developing into preovulatory follicles. In this process, follicle selection (cyclic recruitment) plays a critical role in determining whether prehierarchical follicles enter the follicular hierarchy or undergo atresia in the ovary. During follicle selection, only one dominant follicle is generally recruited into the preovulatory hierarchy on an approximate daily basis, which occurs from a small cohort of prehierarchical follicles measuring approximately 6.0- to 8.0-mm in diameter ([42]Woods and Johnson, 2005; [43]Johnson, 2012). In different ovarian follicle developmental stages, many divergent biological processes affect oocyte growth, and proliferation and differentiation of granulosa and theca cells within the various-sized follicles. Moreover, at the different follicular stages, studies reported ovaries were strictly controlled by multiple distinct reproductive endocrine hormones, paracrine and autocrine regulatory factors, as well as regulatory activities of hypothalamic-pituitary-gonadal axis ([44]Padmanabhan et al., 2002; [45]Johnson and Woods, 2009). Previous studies reported various follicular sizes exhibited different effects on granulosa cell apoptosis, autophagy, and follicular atresia ([46]Lin and Rui, 2010; [47]Zheng et al., 2019). Furthermore, numerous studies established ovarian follicular function is regulated by a plethora of endocrine, paracrine, and autocrine factors besides the pituitary gonadotropins (follicle stimulating hormone and luteinizing hormone) including steroids, neuropeptides, growth factors, adipokines, cytokines, transcription factors, and so on ([48]Woods and Johnson, 2005; [49]Onagbesan et al., 2009; [50]Dupont et al., 2012; [51]Johnson, 2015b; [52]Xu et al., 2018). Obviously, follicular size has a close relationship with its developmental biology, physiological function, and molecular regulation. In recent years, transcriptome analyses revealed and identified many genes relevant to high egg production and which are expressed specifically in ovaries and the glands of the hypothalamus and pituitary of chickens. There are five differentially expressed genes (DEGs) that were identified in the ovaries of Jinghai Yellow chickens related to egg production including ZP2, WNT4, AMH, IGF1, and CYP17A1, which serve as potential candidate genes for improving egg production ([53]Zhang et al., 2019). It was reported that 414, 356, and 10 DEGs from the hypothalamus, pituitary gland, and ovary, respectively, of the Chinese domestic Luhua chicken were potentially associated with greater rates of egg production ([54]Mishra et al., 2020). These results provided valuable information for understanding the reproductive biology of chicken and identifying effective molecular markers for genetic improvement of the domestic chicken breed. Another report has demonstrated high- and low-yielding Chinese Dagu hens possessed 7 and 39 DEGs in the hypothalamus and pituitary, respectively. The main DEGs identified by the authors consisted of OVCH2, SLC7A10, VCAN, CATHL2, GFRA4, ESR1, LAMA1, CFD, and MANBAL1 ([55]Wang and Ma, 2019). Additionally, an earlier study identified transcripts in the hypothalamus and pituitary gland of Taiwan Country chickens including PLAG1, NCAM1, PRL, SAR1A, BDH, PCDHA@, PGDS, SCG2, and STMN2, which could serve as potential molecular markers of high egg production in chicken breeding ([56]Shiue et al., 2006). Several studies in the literature reported similar findings in geese and ducks. One study identified 26 genes (18 known and 8 unknown DEGs) in the ovaries of pre-laying and laying geese, with FSHR, ESR1, ESR2 and PRLR, and FTH serving as representative candidate genes ([57]Kang et al., 2009). [58]Luan et al. (2014) identified 12 candidate genes involved in reproduction regulation that differed during the laying and termination period of egg production of the Huoyan goose breed. This study provided transcriptome profiling data that increased our understanding of geese reproductive biology. [59]Tao et al. (2017) verified 30 genes were relevant to the reproductive process, with 25 associated with egg production in the ovaries of laying ducks. These findings provided useful information that may contribute to future functional studies of genes involved in bird reproduction ([60]Tao et al., 2017). [61]Ren et al. (2019) found seven key DEGs, including EPOX, StAR, CYP17, 3β-HSD, CYP1B1, CYP19A1, and SR-B1, were related to Peking duck ovarian follicle development ([62]Ren et al., 2019). [63]Wu et al. (2020) revealed 135, 56, and 331 genes in Xinjiang Yili geese hypothalamus, pituitary gland, and ovaries, respectively, were related to high and low egg production ([64]Wu et al., 2020). These results provided transcriptome data of the Xinjiang Yili goose hypothalamic-pituitary-gonadal axis, which laid a foundation for further exploration of the molecular mechanisms underlying egg-laying performance of birds. However, the role and molecular regulation DEG play in the formation of various-sized ovarian follicles and egg-laying traits are still poorly understood. The Lohmann Brown (LB) layer, a commercial egg-laying breed, is characterized as possessing a high rate of egg production (HR), high laying peak, and great consistency in laying performance ([65]Lohmann Tierzucht GmBH, 2005). In contrast, the Jilin Black (JB) chicken, an indigenous Chinese breed (dual-purpose type for egg and meat production), is characterized by a low rate of egg production (LR), low laying peak, and strong broodiness ([66]Chen et al., 2004). The divergences in the egg-laying performance between LB and JB, therefore, may serve as a proper model for discovering potential key genes involved in egg production and for further investigating the molecular mechanism underlying formation of egg-laying traits in chicken. To reveal the potential candidate genes associated with egg production, high-throughput transcriptome technology was performed to enrich and analyze differentially expressed genes of follicles sized 3.5–5.5, 6.0–8.0, and 8.5–10.5 mm in diameter between the LB and JB breeds. The objective of this study was to identify key functional genes implicated in ovarian follicular growth and development (including follicle selection) that directly influences egg-laying productivity in hens, and to explore molecular functions, biological processes, and signaling pathways by which the candidate genes contribute to egg-laying production. The present study has laid a foundation for further work ascertaining the exact roles and regulatory mechanisms of the promising genes in ovary development. Finally, these results can be utilized for discovering potential molecular markers of egg-laying traits that can facilitate layer breeding and genetic improvement programs. Materials and Methods Ethics Statement All procedures executed in chickens were approved by the Institutional Animal Care and Use Committee (IACUC) of Jilin Agricultural University (Changchun, China). All birds were sacrificed prior to the removal of organs based upon the IACUC Guidelines for the euthanasia of experimental animals [Permission No. GR (J) 19-030]. Euthanasia of the chickens is conducted fully compliant with Chinese applicable laws and regulations concerning care and use of laboratory animals, which was issued on the basis of the Regulations for the Administration of Affairs Concerning Experimental Animals of the State Council of the People’s Republic of China (1988). All of efforts have been made to minimize the suffering of the animals. Animals and Samples Preparation Lohmann Brown and Jilin Black hens were reared in laying battery cadges according to the husbandry practices as previously reported ([67]Qin et al., 2015). Under the feeding conditions, average body weights of LB and JB hens at the age of 21 weeks were 1.71 ± 0.13 kg and 1.78 ± 0.19 kg, average egg production rate at the age were 53 and 11%, respectively. Ten birds of each breed were obtained from the population, euthanized at 21 weeks of age, and their ovarian follicles were harvested. Growing follicles were categorized by size (3–5 or 6–8 mm) or according to color (large white follicles or small yellow follicles) types ([68]Johnson and Woods, 2009; [69]Johnson, 2012). In this study, white follicles 3.5–5.5 mm in diameter were categorized as large white follicles (LWF), yellow follicles 6.0–8.0 mm were categorized as small yellow follicles (SYF), and yellow follicles 8.5–10.5 mm were categorized as large yellow follicles (LYF, [70]Supplementary Figure 1). After the surrounding vascular and connective tissues of the follicles were removed with fine forceps and a scalpel ([71]Politis et al., 1990), follicles were immediately snapped frozen in liquid nitrogen and preserved at –80°C for RNA extraction. Total RNA was isolated from follicles of each hen using Trizol Reagent (Invitrogen Life Technologies, Carlsbad, CA, United States) according to the recommended manufacturer’s protocol, and concentration, quality, and integrity were determined by a Thermo NanoDrop 8000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, United States). Three micrograms of RNA were used as input material for RNA sample preparations. The other parts of the samples were stored at –80°C for validation experiments. cDNA Library Construction and Illumina Sequencing cDNA libraries were generated using the TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, United States). Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. First strand cDNA was synthesized using random oligonucleotides and SuperScript II after fragmentation was performed utilizing divalent cations under elevated temperature in an Illumina proprietary fragmentation buffer. Second strand cDNA synthesis was subsequently carried out by using DNA polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities and enzymes were removed. Following adenylation of the DNA fragment 3′ends, Illumina PE adapter oligonucleotides were immediately ligated to prepare for hybridization. To screen cDNA fragments 200 bp in length, library fragments were purified by the AMPure XP system (Beckman Coulter, Beverly, CA, United States). DNA fragments with ligated adaptor molecules on both ends were selectively enriched by applying Illumina PCR Primer Cocktail in a 15 cycle PCR reaction. Finally, products were purified (AMPure XP system) and quantified utilizing the Agilent high sensitivity DNA assay on a Bioanalyzer 2100 system (Agilent Biotechnologies, Palo Alto, CA, United States). Effective concentration of libraries were greater than 2 nM and libraries were qualified, sequenced, and paired-end sequencing with the 150 bp sequencing read length was performed. The cDNA libraries were sequenced on a Hiseq platform (Illumina, Inc., San Diego, CA, United States) by Shanghai Personal Biotechnology Cp. Ltd. RNA-Seq Data Analysis To obtain clean reads (clean data), raw data (raw reads) of FASTQ format^[72]1 were processed by filtering out low-quality reads and trimming the adaptor sequences using in-house perl scripts. Meanwhile, Q20, Q30, and N (fuzzy base) content of the raw data were evaluated. Follow-up analyses were based on clean data with high quality. All of the clean reads were aligned with the reference genome^[73]2 by using sequence alignment program HISAT 2.1.0 ([74]Pertea et al., 2016; [75]Kim et al., 2019). Identification of Differentially Expressed Genes HTSeq (ver. 0.6.1) software was used to count the reads mapped to each gene. The gene expression level was normalized based on its FPKM (Fragments per Kilobase of transcript per Million mapped reads) values using the Cufflinks package (v2.1.1). The expression pattern assessment for differentially expressed genes (DEGs) was performed by the Multi-Experiment Viewer (MeV) software (version 4.9.0)^[76]3 based upon the normalized FPKM + 1 value between case group and control group. The resulting p-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate and an absolute value of the |log2FoldChange| (Log2FC) was used as the threshold for judging significance of the gene expression. Candidate genes with an adjusted p-value (padj) < 0.05 and Log2FC > 1 were identified as differentially expressed. Volcano plots to visualize up- and down-regulated genes of each sample were generated by using the R package “ggplot2” (version 3.5.0), based on the values of padj and Log2FC ([77]Wang et al., 2010). To explore the similarity of these samples based on their mRNA expression pattern, all statistical analyses were performed in the R statistical program, and Principle component analysis (PCA) was finished using companion to applied regression (CAR) package in R^[78]4. GO and KEGG Enrichment Analysis for DEGs An online biological tool, Gene Ontology (GO^[79]5), was utilized to annotate and analyze the molecular and functional characteristics of the DEGs ([80]Ashburner et al., 2000). The calculated p-value was Bonferroni corrected, taking the corrected p-value < 0.05 as a threshold for GO annotation, by which the DEGs were enriched in three GO terms including biological process (BP), cellular component (CC), and molecular function (MF) of the GO database. A topGO version 3.0 package^[81]6 provided the tools for testing GO terms while accounting for the topology of the GO graph. Another online biological tool, KEGG (Kyoto Encyclopedia of Genes and Genomes^[82]7), provided the comprehensive database resources for the KEGG pathway enrichment analysis of the DEGs. In this step, four databases were utilized to reveal high-level functions and biological systems of the DEGs, including Reactome^[83]8, KEGG pathway^[84]9, PANTHER^[85]10, and STRING^[86]11. Results with p < 0.05 were considered significantly enriched by DEG. Additionally, the UniProt (UniProt Knowledgebase^[87]12) database, EC (Enzyme Commission^[88]13) database, and eggNOG (Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups^[89]14) database were used for functional and biological analysis of the enriched DEGs. Data Validation by Quantitative Real-Time RT-PCR To verify the accuracy and repeatability of the RNA-Seq results of EDGs, transcription levels of 22 selected genes in the follicles were assessed by using quantitative real-time reverse transcriptase PCR (RT-qPCR) according to our previously described method ([90]Xu et al., 2018). The primers utilized for amplification of the candidate genes including P2RX1, CAB39L, BLK, CSMD3, GPR65, and ADRB2, et al., were listed in [91]Table 1. Using the 2^–ΔΔCt method, mRNA expression results were normalized against18S rRNA as internal control. To quantify mRNA expression levels by RT-qPCR analysis, four amplified products in independent reactions per bird were utilized. All the experiments were carried out in triplicate using different batches of sampled follicles. TABLE 1. Primer pairs designed for quantitative real-time PCR analysis. Gene Forward primer (5′ – 3′) Reverse primer (5′ – 3′) Accession No. Size P2RX1 TTTCCTCTATGAGAAGGGCTAC GGGCAGGTTCCTTGTTTC [92]NM_204519 210 bp CAB39L GATGTCCCAGTGATGTGC AATTCAGCCTCTTGCTCT [93]NM_001006272 150 bp BLK GATGCTGGCGAAACAAAC GGCTCTGCCTCATACTGC [94]XM_004935895 103 bp CSMD3 AGGTGGACAACTCAACAG GGTGCCAACTAAGTAATCT [95]XM_015283035 205 bp GPR65 TGCCTCTGTGGATTGATT TCTTGTGCGTAAGTGCTG [96]XM_421305 185 bp ADRB2 CGCTGTGGACCGATACTTT AGATGATGGAGGAGGCAAT [97]MK139003 239 bp CSMD1 AGGAGGATTCTTTGTTCT GATGGCATTTGTACTTGA [98]XM_015284920 137 bp PLPP4 TTTGGATGGCAGATTGTC CAGTCTAAAGCACCGTAA [99]XM_426544 172 bp ATF3 GGTGACAGTGTCCGAAAG TATCAAATGCTGCTTCTC [100]XM_015283874 238 bp PRLL GCAGTAGATGAAGCGATGT GCCATACCCAACTCCCAC [101]NM_001165912 172 bp STMN3 GAAGTCCGCAGGAACAAG AACGCTGGCATTTAGAACA [102]NM_001004393 144 bp RORB GGTGCTGAGACGGGATTT GGAGTAGGAGGCATTGTTC [103]NM_205093 215 bp PIK3R1 TGCCATTTGTTCCTGTCG GCCCTAAGCATTTATCCC [104]XM_424759 150 bp PERP1 AGACCTTGCCCTATGTGC GAAGTTGAACCGAAGTGTAT [105]NM_001277703 206 bp ACSBG1 TTACATAAGGGAGCATTTC TAACTCATTGCCACCACT [106]XM_413747 103 bp MRTO4 TACACCGTCAGCCTGGAC CTCTGCCATCTCGTAGCC [107]XM_025142548 201 bp CDKN1A GACCACGGAAGGGACTGA AGACGGTGACGGACCACA [108]NM_204396 129 bp EDA2R CCCTGGTCATCATCGTGC GCTTCCACTGGTCCTTCTGT [109]NM_001083360 241 bp GABRA1 TGTGTTTTCTGCCCTCATC ATCCTTCACCTTCTTTGGC [110]NM_204318 109 bp NDUFAB1 AGGACGAGTTCGGCTTTG TGACATCAGGGCCATCC [111]XM_004945263 157 bp NCAM2 CGGCTACAAACAGAATAGGAA ATGAATGAAACCTTGGCAGT [112]XM_425540 124 bp LOC424014 TGAGGATGGCTCGGTTGA AATGGTGTAGGTGCTGATGG [113]XM_001232693 130 bp 18S rRNA ATTGGAGGGCAAGTCTGGTG CCAGCTCGATCCCAAGATCC [114]AF173612 109 bp [115]Open in a new tab Statistical Analysis Statistical analysis was fulfilled using the SPSS12.0 software package ([116]Lin et al., 2002). Data were analyzed by executing a Student’s t-test for comparisons between the RNA-Sequencing and RT-qPCR determination after confirmation of normal distributions for non-parametric analysis. Values were presented as mean ± SEM and bars with superscript symbols that indicate the significant difference compared with control groups at ^∗∗p < 0.01, ^∗p < 0.05. Results Overview of RNA Sequencing Data In this study, a total of 18 cDNA libraries were constructed from HR and LR ovarian follicles ranging in size from 3.5 to 5.5 mm (LWF), 6.0 to 8.0 mm (SYF), and 8.5 to 10.5 mm (LYF) in diameter. Raw reads of each library ranged from 41.61 million up to 53.22 million reads. Among the raw reads, more than 97.35 and 93.87% had quality scores at the Q20 and Q30 levels, respectively. The highest N (base) content of the raw data was 0.0071% ([117]Supplementary Table 1). By quality filtering to strain off low-quality and adaptor sequences, more than 41 million clean reads of each library were procured. In the clean data, clean reads of high quality accounted for 98.91–99.27% of the sequences of the raw reads, and bases (bp) of the clean data accounted for 99.17–99.56% of the raw data ([118]Supplementary Table 2). Total mapped sequences (reads) possessed 85.11–86.49% of the clean reads. In the mapped reads, multiple mapped reads accounted for 3.54–5.25% of the clean reads, and uniquely mapped reads (URs) accounted for 94.75–96.46% of the clean reads ([119]Supplementary Table 3). Among the RNASeq mapped events, total reads mapped to genes (MGs) shared 81.69–83.29% of the URs and mapped to genomic interval regions (interGene) shared 16.71–18.31%. Sequences mapped to exons covered 87.21–91.37% of the MGs ([120]Supplementary Table 4). To quantify gene expression levels, Pairwise Pearson’s correlation coefficients of the samples with their biological replicates were evaluated higher than 0.9179, indicating high repeatability and reliability of the sequencing data ([121]Figure 1A). To obtain an overview of the transcriptomic variation, principal component analysis (PCA) was performed and the values of PC1 and PC2 were 97.95 and 1.37%, respectively ([122]Figure 1B). The genes that led to the maximum amount of variance (PC1) were selected and GO terms obtained by using the GO Consortium. FIGURE 1. [123]FIGURE 1 [124]Open in a new tab 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. Functional Annotation of Unigenes After sequences of unigenes were aligned to the reference genome^[125]15 of the Ensembl database, unigenes were analyzed and annotated by using the NCBI and five other databases ([126]Figure 1C). Total, 18,346 unigenes were annotated in the 18 cDNA libraries. Among them, there were 12,312 unigenes matched successfully in the NCBI databases, accounting for 67.10% of the total unigenes. The current data demonstrated that 15,747 unigenes (85.83%) had successfully matched in the GO database, 9,014 unigenes (49.13%) in KEGG database, 8,761 unigenes (47.75%) in UniProt ID database, and 2,675 unigenes (14.58%) in EC database, and 11,58 unigenes (6.31%) in egg NOG database ([127]Figure 1C). Identification and Analysis of Differentially Expressed Genes (DGEs) In the present study, DEGs in the follicles from HR and LR hens were identified. A total of 236 DEGs out of 16,310 unigenes were found in LWF follicles, of which significantly up- and down-regulated unigenes were 137 and 99, respectively (| log2FoldChange| > 1, P < 0.05; [128]Table 2, [129]Figures 1D, [130]2A, and [131]Supplementary Table 5). A total of 544 (447 up- and 97 down-regulated) and 386 unigenes (229 up- and 157 down-regulated) showed significantly differential expression in SYF and LYF follicles, respectively (| log2FoldChange| > 1, P < 0.05; [132]Figures 2B,C and [133]Supplementary Tables 6, [134]7). The current data indicated the number of total DEGs of SYF follicles between HR and LR hens were much greater than LWF and LYF follicles. TABLE 2. Analysis of differentially expressed genes in various-sized ovarian follicles of HR and LR hens. Follicles Total unigenes Nodiff. Up- Down- DEGs LWF (3.5–5.5 mm) 16310 16074 137 99 236 SYF (6–8 mm) 16366 15822 447 97 544 LYF (8.5–10.5 mm) 16340 15954 229 157 386 [135]Open in a new tab No diff., non-differentially expressed unigenes (|log2FoldChange| < 1 or P > 0.05); Up-, significantly up-regulated unigenes; Down-, significantly down-regulated unigenes (|log2FoldChange| > 1, P < 0.05). FIGURE 2. [136]FIGURE 2 [137]Open in a new tab 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 [MATH: 12 :MATH] 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. GO Functional Classification of DEGs For GO analysis, DEGs were annotated in three ontologies: biological process (BP), cellular component (CC), and molecular function (MF) of the GO database. Based upon comparison of LWF follicles between HR and LR hens, 236 DEGs were assigned to 203 different GO terms covering 162 BP, 179 CC,159 MF, and the top 10 significant GO terms for BP, CC and MF are presented in [138]Figure 3 (P < 0.05, [139]Supplementary Table 8), respectively. When comparing prehierarchical follicles (6.0–8.0 mm in diameter) of HR and LR hens, 544 DEGs were enriched to 486 different GO terms (including 404 BP, 444 CC, 375 MF) and the top 10 significant categorized GO terms are displayed in [140]Figure 4 (P < 0.05, [141]Supplementary Table 9). According to comparison of LYF follicles between HR and LR hens, 386 DEGs were grouped to 336 different GO terms (comprising 246 BP, 301 CC, 241 MF), and top 10 of the significant categorized GO terms are exhibited in [142]Figure 5 (P < 0.05, [143]Supplementary Table 10). FIGURE 3. [144]FIGURE 3 [145]Open in a new tab 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. [146]FIGURE 4 [147]Open in a new tab 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. [148]FIGURE 5 [149]Open in a new tab 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. Comparison of significant GO terms of DEGs in the three follicle groups demonstrated the average top-10 term frequency in SYF follicles from HR and LR hens was much greater than undifferentiated LWF and preovulatory LYF ([150]Figures 3[151]–5). This finding indicated functional annotation data of DEG from SYF follicles clarifies their biological effects on hen egg-laying production. As shown in [152]Figure 4, the main GO terms for SYF molecular functions included, transmembrane signaling receptor activity (GO:0004888), signaling receptor activity (GO:0038023), molecular transducer activity (GO:0060089), receptor activity (GO:0004872), signal transducer activity (GO:0004871), and transmembrane receptor activity (GO:0099600) (P < 0.05, [153]Supplementary Table 9). Moreover, hypotaxis of the significant top 10 GO terms of MF was displayed in the gene ontology DAG structure ([154]Supplementary Figure 2) and revealed that in 486 DEGs out the 12,502 GO termed unigenes by MF annotation, 375 DEGs were involved in the molecular function of GO:003674. Enrichment of Signaling Pathways of DEGs Annotated by GO Terms To reveal what signaling pathways were significantly associated with egg production, and to further identify which DEGs annotated by the GO ontology were implicated in those transducing activities, enrichment analysis of the KEGG pathway was accomplished on the follicles between HR and LR hens. In LWF follicles, most relevant signaling pathways remarkably enriched consisted of cell adhesion molecules (CAM; ko 04514), PPAR signaling pathway (ko 03320), neuroactive ligand-receptor interaction (ko 04080), fatty acid biosynthesis (ko 00061), AMPK signaling pathway (ko 04152), cAMP signaling pathway (ko 04024), adipocytokine signaling pathway (ko 04920), ECM-receptor interaction (ko 04512), p53 signaling pathway (ko 04115), and insulin secretion (ko 04911). Collectively, there were 44 DEGs (24 up- and 20 down-regulated) to be significantly enriched in these 10 KEGG pathways ([155]Supplementary Table 11). In SYF follicles, the prominently enriched KEGG pathways included phagosome (ko 04145), CAMs (ko 04514), the NF-kappa B signaling pathway (ko 04064), the calcium signaling pathway (ko 04020), the Toll-like receptor signaling pathway (ko 04620), neuroactive ligand-receptor interaction (ko 04080), and the GABAergic synapse (ko 04727). In total, 87 DEGs (80 up- and 7 down-regulated) were implicated in the seven signaling pathways ([156]Supplementary Table 12). In LYF follicles, there were 11 relevant KEGG pathways encompassing the neuroactive ligand-receptor interaction (ko 04080), CAMs (ko 04514), circadian entrainment (ko 04713), glutamatergic synapse (ko 04724), mucin type O-glycan biosynthesis (ko 00512), the PPAR signaling pathway (ko 03320), D-arginine and D-ornithine metabolism (ko 00472), carbohydrate digestion and absorption (ko 04973), the JAK-STAT signaling pathway (ko 04630), the cAMP signaling pathway (ko 04024), and cytokine-cytokine receptor interactions (ko 04060). Totally, 71 DEGs (57 up- and 14 down-regulated) were encircled in the 11 signaling pathways ([157]Supplementary Table 13). Obviously, there were two co-enriched pathways (CAMs and neuroactive ligand-receptor interaction) in LWF, SYF, and LYF follicles and two co-enriched pathways (PPAR signaling pathway and cAMP signaling pathway) in LWF and LYF follicles of HR and LR hen ovaries. The relevant up-regulated DEGs included ACSBG1, CNTN1, DMB2, FABP7, GRIK1, GABRA1, GABRB2, GABRG1, GABRQ, GALR2, GRIA2, GRIA4, GRIN2A, HRH1, HRH2, HTR1A, HTR2A, ITGB2, LEPR, LRRC4C, NCAM1, NLGN1, NRCAM, PIK3R1, PRLR, PTPRC, P2RX7, RYR2, VCAM1, VCAN, VIP, and VIPR2, et al. ([158]Supplementary Tables 11–[159]13). The related down-regulated DEGs comprised ACOX2, CCKBR, CDKN1A, CHRNA4, CREB3L3, FABP3, F2RL1, GRPR, MC2R, NCAM2, NRXN1, OPRD1, P2RX1, RBPJL, and SSTR2, et al. ([160]Supplementary Tables 11–[161]13). Characterization of the Potential, Critical DEGs Involved in Egg-Laying Production To further explore critical candidate genes potentially involved in high and low egg-laying chickens, 47 co-expressed DEGs in LWF and SYF follicles of HR and LR hens were identified ([162]Figure 1D); among them, 30 known genes are listed in Table S14 including P2RX1, CAB39L, BLK, CSMD3, GPR65, and ADRB2 ([163]Table 3). Sixty-eight co-expressed DEGs of SYF and LYF follicles were enriched ([164]Figure 1D) and 42 of known genes were focused on ([165]Supplementary Table S15) including CSMD1, PLPP4, ATF3, PRLL, STMN3, and RORB ([166]Table 3). In 54 co-expressed DEGs in LWF and LYF follicles, 35 of known genes were concentrated on ([167]Figure 4 and [168]Supplementary Table 16) and the most relevant genes included PIK3R1, PERP1, ACSBG1, MRTO4, CDKN1A, and EDA2R ([169]Table 3). Moreover, there were 13 co-expressed DEGs obtained in LWF, SYF, and LYF follicles of HR and LR hens ([170]Figure 4 and [171]Supplementary Table 17). Among them, three known genes, one novel gene, and four uncharacterized genes related to egg production were selected from the co-expressed DEGs, the characterized genes included GABRA1, NDUFAB1, NCAM2, LOC424014 ([172]Table 3). In order to lay a useful molecular foundation for further elucidating the genetic mechanism underlying ovarian follicle development associated with egg production in chicken, a potential physical interaction network of proteins encoded by four genes as shown in [173]Figure 6. Additionally, the uncharacterized genes in this study comprised of ENSGALG00000033498, ENSGALG00000036119, ENSGALG00000037683, and ENSGALG00000046112 ([174]Supplementary Table 17). TABLE 3. Characterization of key co-expressed DEGs in LWF, SYF, and LYF follicles of HR and LR hens. Samples Genes RP (H1 vs. L1) RP (H2 vs. L2) RP (H3 vs. L3) KEGG Pathway Britehierarchy (KO) UniProtID PFAM Domains eggNOG LWF&SYF P2RX1 Up Up – Neuroactive ligand-receptor interaction gga04080 [175]Q8AWC0 PF00864 ENOG410IFJF CAB39L Up Up – mTOR signaling pathway gga04150 F1NV80 PF08569 KOG1566 BLK Up Up – Non-receptor tyrosine kinases gga422035 E1BYL7 PF07714 KOG0197 CSMD3 Up Up – Phosphatase activity gga01009 – PF00084 – GPR65 Up Up – GPCR ligand binding gga04030 F1NJ83 PF00001 KOG3656 ADRB2 Up Up – Neuroactive ligand-receptor interaction gga04080 A0A1D5P7A9 PF00001 – SYF&LYF CSMD1 – Up Up Phosphatase activity gga01009 – PF00084 – PLPP4 – Up Up Glycerolipid metabolism gga00561 R4GJ86 PF01569 KOG3030 ATF3 – Up Up cyclic AMP-dependent transcription gga03000 A0A3Q3B2P0 PF00170 – PRLL – Down Down Prolactin signaling pathway gga04917 C6ZDB7 PF00103 ENOG410II74 STMN3 – Down Down Regulation of Rac signal transduction gga04131 [176]O93388 PF 00836 ENOG410IEAZ RORB – Down Up Circadian rhythm gga04710 [177]Q98934 PF00104 KOG4216 LWF&LYF PIK3R1 Up – Up PI3K-Akt signaling pathway gga427171 A0A1D5PUK7 PF16454 – PERP1 Up – Up p53 signaling pathway gga421683 R4GHP1 PF00822 KOG4671 ACSBG1 Up – Up Fatty acid biosynthesis gga00061 F1NLD6 PF00501 KOG1256 MRTO4 Down – Down Ribosome biogenesis gga03009 F1P353 PF00466 KOG0816 CDKN1A Down – Down ErbB signaling pathway gga04012 [178]Q8AYE7 PF02234 – EDA2R Down – Down Cytokine-cytokine receptor interaction gga04060 A2SWM7 PF00020 – LWF&SYF &LYF GABRA1 Up Up Up Neuroactive ligand-receptor interaction gga04080 [179]P19150 PF02931 KOG3642 NDUFAB1 Down Down Down Oxidative phosphorylation gga00190 E1C4C0 PF00550 KOG1748 NCAM2 Down Down Down Cell adhesion molecules gga04514 – PF13895 – LOC424014 Down Up Down putative methyltransferase gga424014 A0A1D5P5L5 PF08241 – [180]Open in a new tab H1 vs. L1, indicating DEGs expressed in LWF follicles of HR and LR hen ovaries; H2 vs. L2, DEGs in SYF follicles of HR and LR; H3 vs. L3, DEGs in LYF follicles of HR and LR; Genes, co-expressed DEGs in the samples; RP, regulation pattern; KO, KEGG Orthology. FIGURE 6. [181]FIGURE 6 [182]Open in a new tab 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. Validation of the Selected DGEs by RT-qPCR After the aforementioned analyses, 22 of most relevant DGEs in the [183]Table 3 were selected for validation by RT-qPCR. As shown in [184]Figure 7, all of the six co-expressed DGEs (P2RX1, CAB39L, BLK, CSMD3, GPR65, and ADRB2) were up-regulated in LWF and SYF follicles sampled from HR and LR hens. For the candidate genes from SYF and LYF follicles, three genes (CSMD1, PLPP4, and ATF3) were up-regulated and two genes (PRLL and STMN3) were down-regulated at the different development stages examined, but candidate RORB was down-regulated at the SYF stage, while up-regulated at the LYF stage. For the candidate genes co-expressed in LWF and LYF follicles, three genes (PIK3R1, PERP1, and ACSBG1) were up-regulated; however, three genes (MRTO4, CDKN1A and EDA2R) were down-regulated. For the DGEs co-expressed in LWF, SYF and LYF follicles, one gene (GABRA1) was up-regulated, two genes (NDUFAB1 and NCAM2) were down-regulated, and a novel gene (LOC424014) was down-regulated in the LWF and LYF follicles, while it was up-regulated in the SYF group. Our results show the significantly differential expression levels of the most transcripts detected by RT-qPCR analysis were consistent with the observations determined by RNA-seq ([185]Figure 7). While there are significant differences in the expression levels of the other genes, including P2RX1, CAB39L, CSMD1, CSMD3, PRLL, NDUFAB1, NCAM2, et al. the observed expression trends coordinate with the RNA-seq results, confirming the transcriptome sequencing results are reliable, and can be further studied for elucidating the genetic mechanism underlying ovarian follicle development associated with egg production in chicken. FIGURE 7. [186]FIGURE 7 [187]Open in a new tab 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. Discussion The ovary is known as one of the most important reproductive organs in adult poultry, which contains many follicles at various developmental stages. Therefore, various regulatory mechanisms in addition to physical contact mediate follicular interactions and directly contribute to egg-laying production ([188]Komatsu and Masubuchi, 2016). The follicle is the major ovarian compartment which enables the ovary to fulfill its dual function of gametogenesis and steroidogenesis ([189]Qin et al., 2015). The various-size growing follicles are mainly composed of oocytes, one or more than two rows of granulosa cells (GCs) surrounding the oocyte, and a theca cell layers at the exterior of the GCs, marking the basement membrane ([190]Johnson and Woods, 2009; [191]Qin et al., 2015). The theca cell layer not only serves as the target apparatus for action of gonadotrophins (FSH and LH) of the hypothalamic-pituitary-gonadal axis that controls the whole reproductive system ([192]Padmanabhan et al., 2002; [193]Woods and Johnson, 2005; [194]Johnson, 2015b), but also acts as a key reproductive moderator and effector to dominate development of oocytes, proliferation and differentiation of granulosa and theca cells, and steroid synthesis (progesterone, estradiol, and testosterone) and secretion ([195]Channing et al., 1980; [196]Richards and Ascoli, 2018). Moreover, a myriad of local autocrine factors and relevant proteins (especially receptors, such as FSHR, ESRs, AMHRs, and VIPR) are synthesized in response to endocrine hormones, paracrine factors, and extra- and intra-ovarian signaling stimuli ([197]Lovell et al., 2003; [198]Xu et al., 2018; [199]Zhang et al., 2019; [200]Mishra et al., 2020). Particularly, genetic regulation of 6.0–8.0 mm diameter follicles are generally involved in follicle selection ([201]Johnson and Woods, 2009; [202]Son et al., 2011; [203]Johnson, 2015b) and may possess an exclusive influence on hierarchy of undifferentiated prehierarchical follicles. Obviously, follicles at various developmental stages have their own distinct molecular biological characteristics and play different roles in contributing to ovary growth and development. Therefore, identification of DEGs of various-sized follicles that differ between HR and LR hens will identify potential candidate genes that control egg production traits. Moreover, our data can provide evidence for highlighting the molecular mechanisms underlying their involvement in formation of egg production traits. Additionally, intensively mining the valuable molecular markers related to egg-laying traits can facilitate breeding programs aiming to improve egg productivity in chickens. In this study, high-throughput transcriptome technology was performed to analyze and identify DEGs of the follicles sized 3.5–5.5, 6.0–8.0, and 8.5–10.5 mm in diameter between the LB and JB chicken breeds. To our knowledge, this is the first study to highlight potential pivotal genes associated with high and low rates of egg production in hen ovarian follicles at various developmental stages via RNA sequencing analysis. Our data showed 47 co-expressed DEGs in LWF and SYF follicles, 68 co-expressed DEGs in SYF and LYF, and 54 co-expressed DEGs in LWF and LYF. These DEGs were significantly related with various physiological functions of follicle growth and development, which have potential to be involved in hen egg-laying production. Moreover, by further enrichment and analyses of the DEGs, 13 of the co-expressed DEGs were obtained in LWF, SYF, and LYF follicles of HR and LR hens. Among them, three known genes (GABRA1, NDUFAB1, and NCAM2) and one novel gene (LOC424014) were characterized. The GABRA1 gene encodes γ-aminobutyric acid type A receptor alpha1 subunit (GABRA1), which is one member of the gamma-aminobutyric acid signaling pathway of the PANTHER subfamily (PTHR18945:SF514). Previously, studies demonstrated the γ-aminobutyric acid type A receptor (GABA-A receptors, GABAR) is a ligand-gated Cl^– ion channel ([204]Morrow et al., 1990), mainly comprising two α subunits, two β subunits, and either one γ or δ subunit ([205]Barrera et al., 2008; [206]Chuang and Reddy, 2018). This configuration shows the receptor’s vital role in apoptosis regulations, mitochondrial function, inhibitory effects of different neurotransmitters, and synaptic plasticity by interaction with neuroactive steroids such as androstane and progesterone derivatives ([207]Tuem and Atey, 2017). Although the exact roles of GABRA1 and the GABAR in ovarian follicular function remain unclear in chicken, a study reported the inhibitory effect of xanthohumol on evoked glutamate release was antagonized by suppressing the GABAR/Ca^2+-calmodulin/AC/cAMP/PKA cascade in rat hippocampus ([208]Chang et al., 2016). Moreover, the AC/cAMP/PKA pathway was involved in granulosa cell proliferation and differentiation of ovarian follicles in hen ([209]Dupont et al., 2012; [210]Nemer et al., 2018). In this study, our results revealed that transcript levels of GABRA1 were remarkably up-regulated in all sizes of ovarian follicles from LR hens compared to HR hens. Therefore, up-regulated GABRA1 gene may be closely related with low hen egg-laying production. This was also supported by the result that expression levels of GABRA1 mRNA were down-regulated in goose ovaries during the laying period when compared with the non-laying period ([211]Luan et al., 2014). The NDUFAB1 gene encodes the acyl carrier protein (ACP), NDUFAB1 subunit of mitochondrial complex I (NADH: ubiquinone oxidoreductase), which is involved in cytosolic and mitochondrial fatty acid synthesis (FAS) pathways. The active form, holo-ACP, serves as a FAS platform, using its 4′-phosphopantetheine group to present covalently attached FAS intermediates to enzymes responsible for the acyl chain elongation process ([212]Masud et al., 2019). Down-regulation of mitochondrial ACP expression in HEK293T cells has brought about a slowing of cell growth rate or a high level of cell death ([213]Feng et al., 2009). Moreover, NDUFAB1 is linked directly to mitochondrial ATP production, indirectly to potential reactive oxygen species (ROS) generation ([214]Pecorelli et al., 2013), and reduced mitochondrial ATP production directly resulted in increased cell apoptosis ([215]Yu et al., 2013). Neural cell adhesion molecule 2 (NCAM2), encoded by the NCAM2 gene, was originally reported to be a key regulator of cell adhesion, axon growth, and fasciculation ([216]Maness and Schachner, 2007). It was also shown to play a critical role in the regulation of granulosa cell differentiation and promoted the development of the cumulus cell phenotype ([217]Eppig et al., 1997; [218]Kaese et al., 2015). In this study, transcript levels of NDUFAB1 and NCAM2 were significantly down-regulated in the ovarian follicles of all size categories from LR hens compared to HR hens. This result suggested down-regulated expression of NDUFAB1 and NCAM2 transcripts may be strongly associated with hen low egg-laying production. The LOC424014 gene encodes a putative methyltransferase that interacts with histone modifications to catalyze DNA methylation to inhibit gene transcription and induce gene silencing ([219]Bird, 2002; [220]Court et al., 2019). This is an uncharacterized novel gene in chicken, which was identified to be significantly down-regulated in the ovarian follicles of 3.5–5.5 and 8.5–10.5 mm in diameters, but up-regulated in 6.0–8.0 mm diameter ovarian follicles of LR hens compared to HR hens. Undoubtedly, these four candidate genes are most potentially associated with egg production traits in chicken, although further studies need to be conducted to confirm this hypothesis. Additionally, the other 18 selected DEGs validated by RT-qPCR, including P2RX1, CAB39L, BLK, CSMD3, GPR65, ADRB2, CSMD1, PLPP4, ATF3, PRLL, STMN3, RORB, PIK3R1, PERP1, ACSBG1, MRTO4, CDKN1A, and EDA2R were significantly differentially expressed in ovarian follicles at the various development stages examined, which indicates they may also be promising molecular biomarkers of chicken egg production. In this study, four signaling pathways including CAMs, the neuroactive ligand-receptor interaction, the PPAR signaling pathway and the cAMP signaling pathway were found to be the most important pathways associated with high and low egg production. Among them, the CAMs pathway was previously reported to mediate many biologic actions such as cell to cell recognition, cell to matrix adhesion, and in early vertebrae embryos ([221]Li et al., 2007). In the current pathway, 15 candidate DEGs including CDKN1A, CNTN1, CREB3L3, DMB2, ITGB2, LRRC4C, NCAM1, NCAM2, NLGN1, NRCAM, NRXN1, PTPRC, RBPJL, VCAM1, and VCAN were initially identified to be involved in ovarian follicle development. Previous studies demonstrated the neuroactive ligand-receptor interaction may be the most important pathway contributing to differential egg production rates between HR and LR hens ([222]Tao et al., 2017; [223]Zhang et al., 2019). In the present study, this pathway was enriched to be involved in chicken ovarian follicle development, of which 47 DEGs including CCKBR, CHRNA4, F2RL1, GABRA1, GABRB2, GABRG1, GABRQ, GALR2, GRIA2, GRIA4, GRIK1, GRIN2A, GRPR, HRH1, HRH2, HTR1A, HTR2A, LEPR, MC2R, OPRD1, PRLR, P2RX1, P2RX7, SSTR2, VIP, and VIPR2 were implicated. Previous studies proved the PPAR signaling pathway plays an important role in ovarian follicle development and normal ovarian function by being indirectly involved in oocyte maturation and ovulation via regulation of steroid hormone synthesis in granulosa cell ([224]Dupont et al., 2012). In the hypothalamus/pituitary gland of laying hens, peroxisome proliferator-activated receptor A (PPARA) and prostaglandin-D synthetase (PGDS) mRNA levels were highly correlated to high egg production in L2 Taiwan country chickens ([225]Chen et al., 2010). In this study, seven DEGs consisting of ACOX2, ACSBG1, ACSBG2, FABP3, FABP7, RXRG, and SCD of the PPAR pathway were found to be co-expressed in LWF and LYF follicles of HR and LR hen ovaries. Moreover, it is well-known that FSH acts with its receptor FSHR in an endocrine dependent manner and plays an essential role in the reproductive system including steroidogenesis, folliculogenesis, and follicular maturation ([226]Choi et al., 2005; [227]Johnson, 2015a). FSHR, a member of the superfamily of G-protein-coupled receptors, is exclusively expressed on granulosa cells of ovarian follicles and mediates FSH signal transduction through cAMP signaling pathway ([228]Simoni et al., 1997; [229]Johnson and Woods, 2009; [230]Johnson, 2015b). cAMP signaling provokes a variety of processes required for growth, selection, differentiation, and maturation of ovarian follicles ([231]Dupont et al., 2012; [232]Richards and Ascoli, 2018). In the current study, 10 DEGs consisting of CREB3L3, GRIA2, GRIA4, GRIN2A, HTR1A, MC2R, PIK3R1, RYR2, SSTR2, and VIPR2 were discovered to be involved in this signaling pathway to control of granulosa cell proliferation and apoptosis, oocyte meiosis and maturation, follicular differentiation and atresia, and secretion of gonadotrophin-release hormone via crosstalk or intracellular interactions with several pathways. These pathways include PPAR signaling, neuroactive ligand-receptor interaction, the calcium signaling pathway, the GnRH signaling pathway, and the PI3K/AKT signaling pathway ([233]Chang et al., 2016; [234]Nemer et al., 2018; [235]Richards and Ascoli, 2018; [236]Liu et al., 2019). Taken together, our present data provide additional evidence for supporting these four signaling pathways were related to chicken egg production. In summary, this study represents the first transcriptome analysis of hen ovarian follicles at different developmental stages between chickens that possess high and low rates of egg production. It established a foundation for further investigation of regulatory mechanisms of the candidate genes involved in ovary follicular development, for understanding the molecular regulation of egg-laying traits, and discovering potential molecular markers that can facilitate chicken breeding programs to enhance egg-laying productivity. Data Availability Statement Our raw RNA sequencing data is available at BioProject under accession numbers [237]PRJNA669967 and [238]PRJNA670553. Ethics Statement All procedures executed in chickens were approved by the Institutional Animal Care and Use Committee (IACUC) of Jilin Agricultural University (Changchun, China). All birds were sacrificed prior to the removal of organs based upon the IACUC Guidelines for the euthanasia of experimental animals [Permission No. GR (J) 19-030]. Author Contributions XC, XS, and RX conceived and designed the experiments. XC, XS, IC, NQ, XX, SL, and RX performed the experiments. XC, RX, XX, and JG analyzed the data. XC and RX wrote the manuscript. JG reviewed and edited the manuscript. All the authors contributed to the article and approved the submitted version. Conflict of Interest 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. Acknowledgments