{"paper_id":"d03bbfbb-8eea-429a-a75b-86d745bf9c52","body_text":"Egg production is an important indicator of the reproductive performance of poultry and an economically important trait. Eggs are an economical source of high-quality protein and important vitamins and minerals that play a vital role in the human diet [ 1 ]. Eggs are consumed worldwide, and their consumption will continue to grow in the future due to their low price and nutritional benefits [ 2 ]. The ovary is a vital component of the female reproductive system. Its main roles are to produce eggs and regulate oocyte secretion by releasing sex hormones for ovulation, and the secretion of hormones essential for regulating oocyte development, follicle maturation, and the ovulation process [ 3 ]. Eggs develop from follicles, and the process of egg formation in the ovary is a complicated process that is influenced by environmental and nutritional factors [ 4 ], as well as by the hypothalamic-pituitary-gonadal (HPG) axis, which produces specific neuropeptides or hormones that stimulate oocyte maturation and ovulation [ 5 ]. Furthermore, ovarian follicular development also involves various endocrine, autocrine, and paracrine factors, which control the proliferation and differentiation of oocytes, granulosa cells, and thecal cells [ 6 ]. RNA sequencing (RNA-seq) is a transcriptome analysis approach that provides information at the single nucleotide level regarding the complete biological transcript [ 7 ]. It utilizes deep sequencing techniques to address biological questions and allows the detection of differences in gene expression profiles between individuals of different developmental stages at the transcriptional level of whole-genome sequencing, facilitating the study of gene function. Transcriptome sequencing has been conducted in numerous studies to identify ovarian genes in livestock such as pigs [ 8 ], cattle [ 9 ], sheep [ 10 ] and goats [ 11 ]. Zou et al. identified 1027 differentially expressed genes in the ovaries of high- and low-yielding Leizhou black ducks, of which 495 genes were up-regulated and 532 genes were down-regulated, and found fifty genes that were related to reproduction and reproductive processes [ 12 ]. Hu et al., identified 1683 DEGs (1136 up-regulated and 1045 down-regulated) in the ovaries of black Muscovy ducks at the early, peak and late egg-laying stages; they predicted that the  HOXA10, HtrA3, StAR, ZP2  and  TAT  genes were involved in regulating ovarian development at different laying stages, and several important pathways were found, including steroid hormone biosynthesis and ovarian steroidogenesis [ 13 ].\nNinghai indigenous chickens are commonly raised in southeastern Zhejiang Province, being a native chicken with high economic value for both meat and eggs. However, the poor egg-laying performance and short peak egg-laying period restrict the economic benefits of enterprises relating to this type of chicken. In this study, RNA-seq was used to identify differentially expressed genes in different egg-laying stages of Ninghai indigenous chickens. Through comprehensive analysis of the differential genes and pathways, in order to identify key pathways and candidate genes involved in the control of egg laying in Ninghai indigenous chickens, we aim to provide a basis for subsequent improvement of egg-laying performance, and also to help understand the molecular regulatory mechanisms of the egg-laying characteristics of chickens.\n\nTwelve Ninghai chickens (four 145-day-old chickens, four 204-day-old chickens and four 300-day-old chickens) were purchased from the Poultry Breeding Center of Ningbo Zhenning Animal Husbandry Co., Ltd. (Zhejiang, China) in Zhejiang Province. All the chickens were fed according to the same housing and feeding conditions. Ovarian tissues were collected from four randomly selected chickens at the 145-day-old, 204-day-old and 300-day-old, stages, respectively, and were immediately frozen in liquid nitrogen.\nThe animal study protocol followed the Chinese Animal Welfare Guidelines, and was approved by the Animal Welfare Committee of Zhejiang University (approval number: ZJU20190149).\nThe collected hen ovarian tissues were delivered to the Novogene Bioinformatics Technology Co., Ltd. (Beijing, China), who conducted all the library preparation and sequencing. For the RNA sample preparations, a total of 3 g of RNA per sample was employed as the input material. Following the manufacturer’s instructions, sequencing libraries were created using the NEBNext ®  UltraTM RNA Library Prep Kit for Illumina ®  (NEB, California, USA), and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer(5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure was ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 250~300 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, California, USA). Then, 3 µL USER Enzyme (NEB, California, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. Then, PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. Finally, the PCR products were purified (AMPure XP system, California, USA) and the library quality was assessed on the Agilent Bioanalyzer 2100 system.\nIn order to control the quality, in-house Perl scripts were used to process the raw data (raw readings) in the FASTQ format [ 14 ]. Clean data (clean reads) were acquired at this stage by eliminating adapter-containing reads, ploy-N containing reads, and low-quality reads from the raw data. The clean values for Q20, Q30, and GC content were computed at the same time [ 15 ]. All of the subsequent analyses were based on clean, high-quality data. The annotation data for the reference genome and gene models were acquired directly from the genome website ( http://ftp.ensembl.org/pub/release-105/fasta/gallus_gallus/ ,  http://ftp.ensembl.org/pub/release-105/gtf/gallus_gallus/ , accessed on 10 June 2020). Hisat2 v2.0.5 was used to create a reference genome index, and Hisat2 v2.0.5 was used to align paired-end clean reads to the reference genome [ 16 ]. FeatureCounts v1.5.0-p3 was used to count the reads numbers mapped to each gene [ 17 ]. Then, the FPKM of each gene was calculated based on the length of the gene and the reads count mapped to this gene [ 18 ].\nThe DESeq2 R package (1.16.1) was used to perform differential expression analyses on four biological replicates per condition [ 19 ]. Using a model based on the negative binomial distribution, DESeq2 provides statistical methods for identifying differential expression in digital gene expression data. The false discovery rate was controlled by adjusting the  p -values using the Benjamini–Hochberg method [ 20 ]. Genes discovered by DESeq2 with an adjusted  p -value of less than 0.05 were labeled as differentially expressed. The read counts for each sequenced library were modified using the edgeR software package through one scaling normalization factor prior to differential gene expression analysis. The edgeR R package (3.18.1) was used to perform differential expression analysis of two situations. The Benjamini–Hochberg method was used to alter the  p -values. The threshold for substantially differential expression was established at a corrected P-value of 0.05 and an absolute fold-change of two.\nThe clusterProfiler R package was used to perform gene ontology (GO) enrichment analysis of differentially expressed genes, with gene length bias adjusted [ 21 ]. Differentially expressed genes substantially enriched GO keywords with a  p -value of less than 0.05.\nKEGG ( http://www.genome.jp/kegg/ , accessed on 15 August 2020) is a database resource for deducing high-level functions and utilities of biological systems, such as the cell, the organism, and the ecosystem, from molecular-level data, particularly large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies [ 22 ]. To examine the statistical enrichment of differentially expressed genes in KEGG pathways, we utilized the clusterProfiler R program.\nUsing total RNA, 1 µg of total RNA was reverse transcribed to cDNA using the Prime Script RT kit (Takara, Hangzhou, China). qPCR was performed using SYBR Green PCR Master Mix (TaKaRa, HangZhou, China) on a StepOnePlus Real-Time PCR System. The 2 −ΔΔCt  method was used to calculate the relative expression levels of genes [ 23 ], using β-actin as a control. Three biological replicates were used to analyze all mRNA expression. The primers used in the qPCR were designed using the NCBI website ( https://www.ncbi.nlm.nih.gov/tools/primer-blast/ , accessed on 11 March 2022) ( Table S1 ).\n\nA total of 701,634,762 clean reads were obtained from the 12 libraries with an average of 584,695,63.5 clean reads in each sample (the numbers of reads ranging from 51,709,354 to 66,121,522). The Q30 and GC content were 92.32% to 94.86% and 50.4% to 51.78%, respectively. More than 83.85% of the reads could be mapped to the chicken genome (NH145d, with 83.85% to 92.24% alignment; NH204d, with 88.69% to 91.12% alignment; and NH300d, with 88.64% to 91.52% alignment) and the unique map reads ranged from 82.17% to 90.94% ( Table 1 ).\nA total of 8433 DEGs were acquired from the three groups, with 1176 up-regulated genes and 621 down-regulated in NH145d vs. NH204d, 1653 up-regulated genes, and 1955 down-regulated in NH145d vs. NH300d, and 1868 up-regulated genes, and 1160 down-regulated in NH300d vs. NH204d. In addition, the two comparison groups included 73, 1004 and 1030 DEGs, respectively; and 65 genes were co-expressed among the three comparison groups ( Figure 1  and  Figure 2 ,  Table S2 ).\nDEGs from the three comparison groups were enriched for GO and KEGG analysis, in order to further elucidate the contribution of the specific signaling pathways in ovarian development. In GO analysis, the DEGs were annotated into three ontologies of the GO database: biological process (BP), cell component (CC), and molecular function (MF). In NH145d vs. NH204d, 3624 DEGs were annotated to 669 GO terms, covering 343 biological processes (BP), 78 cellular components (CC), and 248 molecular functions (MF). The items of evident enrichment included immune response, chemokine activity, chemokine receptor binding, cytokine activity, plasma membrane, and G-protein coupled receptor binding. In NH145d vs. NH300d, 7060 DEGs were assigned to 857 GO terms, including 448 BP, 118 CC, and 291 MF. The items of significant enrichment included the activity of passive transmembrane transporter activity, chemokine activity, chemokine receptor binding, and growth regulation. In NH300d vs. NH204d, 7244 DEGs were classified into 858 different GO terms, including 455 BP, 113 CC, and 290 MF. The GO terms included multi-organism cellular process, multi-organism process, and extracellular matrix structural constituent ( Figure 3 ,  Table S3 ).\nFor KEGG analysis, in NH45d vs. NH204d, 632 genes were enriched in 118 pathways, of which 16 pathways were significantly enriched, mainly in the cell adhesion molecules (CAMs), cytokine–cytokine receptor interaction, herpes simplex virus 1 infection, intestinal immune network for IgA, and the calcium signaling pathway. In NH145d vs. NH300d, 987 genes were enriched in 143 pathways, of which nine pathways were significantly enriched, mainly in the neuroactive ligand–receptor interaction, cytokine–cytokine receptor interaction, vascular smooth muscle contraction, and ABC transporters. In group NH300d vs. NH204d, 1140 genes were annotated in 148 pathways, of which four pathways were significantly enriched, namely progesterone-mediated oocyte maturation, RNA degradation, oocyte meiosis, and homologous recombination, respectively. Among the significantly enriched pathways, four pathways overlap in the two comparison groups NH145d vs. NH204d and NH145d vs. NH300d, namely cytokine–cytokine receptor interaction, the calcium signaling pathway, regulation of the actin cytoskeleton, vascular smooth muscle contraction, and neuroactive ligand–receptor interactions ( Table 2  and  Figure 4 ,  Table S4 ).\nWe performed RT-qPCR analysis of the expression of the nine genes at different times. As shown in  Figure 5 A, the qPCR expression results were consistent with the trend of RNA-seq results.\n\nImproving egg production is an important goal in poultry breeding. The ovary is a critical organ associated with poultry reproduction and is an essential tissue for finding candidate genes relating to egg production. By investigating ovarian gene expression patterns during development, we can contribute to improving egg production performance and to a better understanding of the reproductive physiology of poultry. In a comparative transcriptome analysis of ovarian tissues from high- and low-laying black and white Muscovy ducks, nine genes including  TGFβ2, NGFR, CEBPD, CPEB2, POSTN, SMOC1, FGF18, EFNA5  and  SDC4  were considered to be closely associated with egg production [ 24 ]. Transcriptome analysis of duck ovaries at different egg production rates revealed that genes such as  MC5R, APOD, ORAI1  and  DYRK4  were more active in the ovaries of high-laying ducks, and that the steroid biosynthetic pathway, calcium reabsorption pathways regulated by endocrine and other factors, circadian rhythms, neuroactive ligand–receptor interaction pathway, fatty acid biosynthesis, and calcium signaling pathways had more important roles in high-laying ducks [ 7 ]. Zhang et al.’s [ 25 ] comparative transcriptomic analysis of ovaries from high- and low-egg-laying Lingyun black-bone chickens found four genes,  FOXA2, MED37D, HSP70 , and  RXFP2 , and three signaling pathways that may be related to egg production, namely the longevity-regulating pathway, multiple species pathway, embryonic signaling pathway and PPAR signaling pathway.\nThe low-density lipoprotein (LDL) receptor-related protein 8 ( LRP8 ) is a member of the LDL receptor family that has a role in endocytosis and signal transduction;  LRP8  paracrine interaction regulates follicular growth [ 26 ]. This gene is a component of the selenium delivery pathway to spermatogenic cells of mice [ 27 ],  LRP8  deficiency causes infertility in male mice, and it is necessary for sperm maturation [ 28 ]. This gene is also related to the development of the ovaries and follicles; it can affect the reproduction of female animals by regulating ovarian or follicular development [ 29 , 30 ]. In addition, this gene has also been reported to be a marker of fertility and follicle maturation [ 31 ]. In this study,  LRP8  was found to be significantly expressed at the peak of egg production, suggesting that LRP8 may have a positive role in regulating follicle growth in the ovary. The bone morphogenetic protein 6 ( BMP6 ) is expressed in the ovary and is recognized as an autocrine/paracrine regulator of follicular and luteal cell proliferation and steroidogenesis [ 32 ].  BMP6  is abundantly expressed in primordial, primary and secondary follicles [ 33 ]. During the transition from primordial to primary/secondary follicles,  BMP6  increases and enhances the growth of cultured secondary follicles [ 34 ]. In infertile women with endometriosis, the expression of the  BMP6  gene is relatively reduced [ 35 ]. A study found that  BMP6  regulates follicle dynamics and granulosa cell differentiation in the turkey ovary [ 36 ]. Zona pellucida glycoprotein 4 ( ZP4 ) is a member of the zona pellucida (ZP) family, which is an extracellular glycoprotein matrix that surrounds oocytes [ 37 ]. Mutations in  ZP4  are associated with abnormal zona pellucida and with female infertility [ 38 ]. It also plays a vital role in fertilization by functioning as a “docking point” for spermatozoa binding, which is followed by the stimulation of the acrosome response in the zona-attached sperm [ 39 ]. In summary,  LRP8, BMP6 , and  ZP4  were significantly upregulated in the NH204d group compared to the NH145d group, suggesting that  LRP8, BMP6  and  ZP4  may have important regulatory roles in egg production mechanisms during peak egg production.\nThe versican ( VCAN ) protein is an extracellular matrix proteoglycan that stabilizes hyaluronan in the expanded extracellular matrix [ 40 ]. It has been found that the in vitro matured mouse mound-oocyte complex is deficient in  VCAN , which may contribute to the poor health of in vitro matured oocytes and embryos [ 41 ].  VCAN  expression has been reported as a marker of pregnancy [ 42 ]. Furthermore, Shen et al. [ 43 ] found the level of  VCAN  expression in cumulus cells to be positively correlated with the early embryo morphology score. The α1(IV) chain encoded by the  COL4A1  (type IV collagen α1 chain) gene is ubiquitously expressed and essential for the basement membrane’s stability [ 44 ]. It is strictly related to ovarian function and follicular development [ 45 ].  COL4A1  has an essential role in regulating ECM remodeling during initial ovarian development in the bovine ovary and in regulating bovine fetal ovaries throughout pregnancy [ 46 ]. The inhibin β A subunit ( INHBA ) is a member of the transforming growth factor-β (TGFβ) superfamily, which encodes the βA-subunit of several activin and inhibin complexes [ 47 ].  INHBA  regulates proliferation, apoptosis, and hormone synthesis in granulosa cells. Knocking out  INHBA  and inhibin subunit βB ( INHBB ) in granulosa cells impacts activin and inhibin production, resulting in sterile mice with a complex ovarian phenotype [ 48 ]. In the present study,  VCAN, COL4A1  and  INHBA  were upregulated in the NH145d group compared to the NH300d group, and the results suggest that the three genes may contribute to egg production performance in the early stages of laying, and may play a regulatory role in follicle development as well as ovarian matrix development.\nLysyl oxidase ( LOX ) activity catalyzes the final enzymatic reaction required in the biosynthesis of cross-linked mature collagens and elastin and, as such, is critical in the formation and deposition of the extracellular matrix (ECM) [ 49 ]. It has been reported that  LOX  is identified as part of the coordinated regulatory loop of the rat ovarian endocrine, paracrine and autocrine levels during follicular development [ 50 ]. Liu et al. highlighted the role of local cortisol in the amnion in downregulating  LOX  gene expression through a negative response element in the  LOX  promoter, with impaired  LOX  functions contributing to fetal membrane rupture [ 51 ]. Harlow et al. [ 52 ] found an essential role in regulating follicle development by  LOX  in the ovary. Pentraxin 3( PTX3 ) is a prototypic long pentraxin, a secreted protein member of the pentraxin family [ 53 ]. It plays a crucial role in organizing the cumulus oophorus extracellular matrix and in in vivo fertilization. Inactivating  PTX3  in mice can disrupt cumulus oophorus formation around the oocyte and reduce the fertilization rate. Zhang et al. [ 54 ] suggest that the gene expression in cumulus cells indicates pentraxin 3 as a possible marker for oocyte quality; Li et al. [ 55 ] found that microRNA-224 delays oocyte maturation through targeting  PTX3  in cumulus cells. Indian hedgehog ( IHH ), a member of the hedgehog (HH) family, plays a significant role in regulating numerous developmental processes, and is hormonally controlled and associated with co-maturation of the theca interna in the mammalian ovary [ 56 , 57 ]. Ovaries lacking the  IHH  show loss of the theca layer, blunted steroid production, impaired folliculogenesis, and failure to form corpora lutea [ 58 ].  IHH  in epithelial cells generally acts as a paracrine growth factor for stromal cells in early gestation, and plays a role in peri-pregnancy preparation for uterine implantation [ 59 ].  IHH  is also an essential mediator of progesterone signaling in the uterus, and expression of this factor is critical in mediating the communication between the uterine epithelium and stroma required for embryo implantation [ 60 ]. In this study,  LOX, PTX3 , and  IHH  were significantly increased in the NH204d group compared with the NH300d group, which suggests that the three genes are involved in the regulation of laying in the peak period, and may play an important role in ovarian hormonal regulation, promoting oocyte maturation.\nGO (gene ontology) annotation and KEGG (Kyoto Encyclopedia of Genes and Genomes) analyses were used to further elucidate the biological roles of the DEGs. The GO results in the three comparison groups show that the biological processes and molecular functions in the ovary were allocated more DEGs than the cellular components. Chemokine activity and chemokine receptor binding were significantly enriched in NH145d vs. NH204d and NH145d vs. NH300d groups. Chemokines exert their effect by binding to their appropriate receptors, the expression levels of which may modulate their action [ 61 ]. The expression of chemokine and chemokine receptors in the endometrium suggests that autocrine and paracrine interactions involving these chemokines participate in endometrial physiology [ 62 ]. Chemokines are expressed in ovarian follicles and oocytes, potentially regulating follicular function and affecting reproductive lifespan [ 63 ]. The extracellular matrix is a fundamental structure present in all tissues of animal organisms. It exhibits various functions and is composed of many substances of different origins [ 64 ]. In addition, it has been shown that extracellular matrix components directly affect biological processes within the ovary, such as folliculogenesis, ovulation, and steroidogenesis [ 65 ]. Findings reported in the past few years indicate that the cumulus matrix plays a key role in the early events of in vivo fertilization [ 66 ]. Furthermore, the extracellular region and plasma membrane also play essential roles in the ovary; enrichment of the plasma membrane in the ovary is consistent with a recent report on the white Muscovy duck [ 67 ]. In a study of human sperm proteome profiles, GO analysis indicated that most of the differently identified sperm proteins were enriched in extracellular membrane-bounded organelles [ 68 ]. Moreover, “immune response” and “immune system process” were the biological processes with the highest DEGs, and they may play an essential role in the ovary of Ninghai chicken.\nFor KEGG, the signaling pathways, including cell adhesion molecules (CAMs), neuroactive ligand–receptor interaction, the calcium signaling pathway, and cytokine–cytokine receptor interaction, were the four most essential pathways associated with the ovary in Ninghai chickens. Cell adhesion is critical for ovary function and follicle development via interactions with the extracellular matrix and direct cell–cell contacts. The CAMs pathway has previously been shown to mediate a wide range of physiological functions, including cell-to-cell recognition, cell-to-matrix adhesion, and the development of early vertebrate embryos [ 69 ]. Furthermore, the calcium signaling pathway and neuroactive ligand–receptor interaction have been reported in more poultry ovaries; calcium (Ca 2+ ) is a signaling molecule that regulates a wide range of biological functions [ 70 ], and the calcium signaling pathway plays an essential role in the egg production of ducks, while the pathway was also enriched in egg production in goose ovaries [ 71 ]. Ye et al. [ 72 ] found that the neuroactive ligand–receptor interaction pathway and the calcium signaling pathway were the key pathways controlling duck brooding. Another study also shows that this pathway affects goose egg production through ovarian metabolic function. In addition, the enrichment of the cytokine–cytokine receptor interaction pathway in the ovary in this study is consistent with the findings of Jinghai Yellow chicken, which also included the calcium signaling pathway and neuroactive ligand–receptor interaction [ 38 ]. In summary, several reports have shown the importance of cell adhesion molecules (CAMs), neuroactive ligand–receptor interaction, calcium signaling pathways, and cytokine–cytokine receptor interaction in chicken egg production [ 24 , 73 ]. In addition to these four pathways, we identified common KEGG pathways involved in reproduction during three periods ( Figure 5 B), including ECM-receptor interactions, steroid biosynthesis and progesterone-mediated oocyte maturation, three of which were found to be closely related to egg production. However, the other six pathways were rarely reported, and although not significantly enriched during all three periods, they may play an important role in poultry during egg production. In general, we have provided some new insights into the pathways associated with poultry egg production.\n\nIn this study, we identified 8433 differentially expressed genes through transcriptome sequencing analysis of the ovaries of Ninghai chickens at different egg-laying periods. We predicted nine genes involved in ovarian development— LRP8, BMP6, ZP4, VCAN, COL4A1, INHBA, LOX, PTX3 , and  IHH —through the calcium pathway and neurotransmitter receptor interaction pathways. The results of this study will help in the investigation of the regulatory network and molecular mechanisms of ovarian development in Ninghai chickens, and will provide valuable information for understanding the mechanisms of ovarian action and the molecular approaches towards improving egg production performance in chickens.","source_license":"CC-BY-4.0","license_restricted":false}