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

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.


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 (Etches et al., 1990;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 (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 (Lovell et al., 2003;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 (Woods and Johnson, 2005;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 (Padmanabhan et al., 2002;Johnson and Woods, 2009). Previous studies reported various follicular sizes exhibited different effects on granulosa cell apoptosis, autophagy, and follicular atresia (Lin and Rui, 2010;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 (Woods and Johnson, 2005;Onagbesan et al., 2009;Dupont et al., 2012;Johnson, 2015b;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 . 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 (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 (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 (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 (Kang et al., 2009). 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. 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 (Tao et al., 2017). 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 (Ren et al., 2019). 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 (Wu et al., 2020). These results provided transcriptome data of the Xinjiang Yili goose hypothalamicpituitary-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 (Lohmann Tierzucht GmBH, 2005). In contrast, the Jilin Black (JB) chicken, an indigenous Chinese breed (dualpurpose type for egg and meat production), is characterized by a low rate of egg production (LR), low laying peak, and strong broodiness (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.

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 (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 (Johnson and Woods, 2009;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, Supplementary  Figure 1). After the surrounding vascular and connective tissues of the follicles were removed with fine forceps and a scalpel (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 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 2 by using sequence alignment program HISAT 2.1.0 (Pertea et al., 2016;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) 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 (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 4 .

GO and KEGG Enrichment Analysis for DEGs
An online biological tool, Gene Ontology (GO 5 ), was utilized to annotate and analyze the molecular and functional characteristics of the DEGs (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 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 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 8 , KEGG pathway 9 , PANTHER 10 , and STRING 11 . Results with p < 0.05 were considered significantly enriched by DEG. Additionally, the UniProt (UniProt Knowledgebase 12 ) database, EC (Enzyme Commission 13 ) database, and eggNOG (Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups 14 )

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 (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 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.

Statistical Analysis
Statistical analysis was fulfilled using the SPSS12.0 software package (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.

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% (Supplementary Table 1). By quality filtering to strain off lowquality 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 (Supplementary Table 2  their biological replicates were evaluated higher than 0.9179, indicating high repeatability and reliability of the sequencing data ( 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 ( Figure 1B). The genes that led to the maximum amount of variance (PC1) were selected and GO terms obtained by using the GO Consortium.

Functional Annotation of Unigenes
After sequences of unigenes were aligned to the reference genome 15 of the Ensembl database, unigenes were analyzed and annotated by using the NCBI and five other databases ( Figure 1C) 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 15 http://asia.ensembl.org/Gallus_gallus/Info/Index unigenes were found in LWF follicles, of which significantly up-and down-regulated unigenes were 137 and 99, respectively (| log2FoldChange| > 1, P < 0.05;  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 (Figures 3-5). This finding indicated functional annotation data of DEG from SYF follicles clarifies their biological effects on hen egg-laying production. As shown in 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, Supplementary Table 9). Moreover, hypotaxis of the significant top 10 GO terms of MF was displayed in the gene ontology DAG structure (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 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).

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 ( Figure 1D); among them, 30 known genes are listed in Table  S14 including P2RX1, CAB39L, BLK, CSMD3, GPR65, and ADRB2 (Table 3). Sixty-eight co-expressed DEGs of SYF and LYF follicles were enriched ( Figure 1D) and 42 of known genes were focused on (Supplementary Table S15) including CSMD1, PLPP4, ATF3, PRLL, STMN3, and RORB (Table 3). In 54 co-expressed DEGs in LWF and LYF follicles, 35 of known genes were concentrated on (Figure 4 and Supplementary Table 16) and the most relevant genes included PIK3R1, PERP1, ACSBG1, MRTO4, CDKN1A, and EDA2R (Table 3). Moreover, there were 13 co-expressed DEGs obtained in LWF, SYF, and LYF follicles of HR and LR hens (Figure 4 and 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 (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 Figure 6. Additionally, the uncharacterized genes in this study comprised of ENSGALG00000033498, ENSGALG00000036119, ENSGALG00000037683, and ENSGALG00000046112 (Supplementary Table 17).

Validation of the Selected DGEs by RT-qPCR
After the aforementioned analyses, 22 of most relevant DGEs in the Table 3 were selected for validation by RT-qPCR. As shown in 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 upregulated, two genes (NDUFAB1 and NCAM2) were downregulated, 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 (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.

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 (Komatsu and Masubuchi, 2016). The follicle is the major ovarian compartment which enables the ovary to fulfill its dual function of gametogenesis and steroidogenesis (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 (Johnson and Woods, 2009;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 (Padmanabhan et al., 2002;Woods and Johnson, 2005;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 (Channing et al., 1980;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 (Lovell et al., 2003;Xu et al., 2018;Zhang et al., 2019;Mishra et al., 2020). Particularly, genetic regulation of 6.0-8.0 mm diameter follicles are generally involved in follicle selection (Johnson and Woods, 2009;Son et al., 2011;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 egglaying 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 egglaying 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 (Morrow et al., 1990), mainly comprising two α subunits, two β subunits, and either one γ or δ subunit (Barrera et al., 2008;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 (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 (Chang et al., 2016). Moreover, the AC/cAMP/PKA pathway was involved in granulosa cell proliferation and differentiation of ovarian follicles in hen (Dupont et al., 2012;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 (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 (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 (Feng et al., 2009). Moreover, NDUFAB1 is linked directly to mitochondrial ATP production, indirectly to potential reactive oxygen species (ROS) generation (Pecorelli et al., 2013), and reduced mitochondrial ATP production directly resulted in increased cell apoptosis (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 (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 (Eppig et al., 1997;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 (Bird, 2002;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.
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 (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 (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 (Choi et al., 2005;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 (Simoni et al., 1997;Johnson and Woods, 2009;Johnson, 2015b). cAMP signaling provokes a variety of processes required for growth, selection, differentiation, and maturation of ovarian follicles (Dupont et al., 2012;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 (Chang et al., 2016;Nemer et al., 2018;Richards and Ascoli, 2018;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 PRJNA669967 and 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