Transcriptome analysis reveals key genes and pathways associated with egg production in Nandan-Yao domestic chicken.

BACKGROUND
Egg production is a very important economic trait in chicken breeding, but its molecular mechanism is unclear until now. Nandan-Yao chicken (Gallus gallus domesticus) is a native breed in Guangxi province, China, which is famous for good meet quality, but with low egg production.


METHODS
To explore the molecular regulation related to egg production, high egg production (HEP) and low egg production (LEP) were divided according to the total egg number at 55 weeks, and the concentration of serum sex hormones was tested to evaluate the physiological function of ovary and uterus. RNA sequencing (RNA-Seq) was used to explore the transcriptome from the ovary and uterus of Nandan-Yao chicken.


RESULTS
The levels of serum sex hormone showed that concentrations of estradiol (E2), follicle-stimulating hormone (FSH), and luteotropic hormone (LH) were significantly higher in HEP than those in LEP (P < 0.01), while the concentration of testosterone (T) was significantly lower in HEP (P < 0.01). RNA-Seq analysis identified 901 and 2763 differentially expressed genes (DEGs) in ovary and uterus, respectively. Enrichment analysis showed that DEGs were significantly involved in the regulation of tight junction in the ovary (P < 0.05), while in uterus, DEGs were mainly enriched in the phagosome, ECM-receptor interaction, cell adhesion molecules (CAMs), focal adhesion, cardiac muscle contraction, cytokine-cytokine receptor interaction, and the regulation of MAPK signaling pathway (P < 0.05). Protein network interaction and function analyses revealed that FN1, FGF7, SOX2 identified from the ovary, and UQCRH, COX5A, FN1 from the uterus might be key candidate genes for egg production in Nandan-Yao chicken.


CONCLUSIONS
Our study provided key candidate genes and pathways involved in the egg-laying process of Nandan-Yao chicken and could help to further understand the molecular mechanisms of chicken reproduction.


Background
Increasing egg production is one of the main purposes of breeding in laying hens [1]. The selection of egg-laying birds based on partial egg production can increase total egg production [2]. The interval of generations per time unit can be reduced by half if the selection is done considering partial egg production according to the previous study [3]. The selection of partial, or annual, egg counts or egg production rates is a common method for improving the laying performance of hens, which can produce some positive genetic progress but it is slow [4]. However, most of the chicken laying traits, including egg production, egg weight, and age of sexual maturity, are polygenic traits with low to moderate heritability, making it di cult to estimate the level of genetic improvement in each generation [5].
The quest for improved egg production is an important focus of poultry breeding. A large number of people have done a lot of research on egg production through different technologies. For example, the genome-wide association study identi ed many mutants associated with egg production by chickens [6].
Wang CQ et al (2019) using RNA-seq for high yield and low yield Dagu China chicken (CDC) of the hypothalamus and pituitary expression pro les for sequencing and analysis of the transcriptome, functional annotation and pathway enrichment analysis showed that DEGs mainly rich in sugar glycosaminoglycans biosynthesis, extracellular matrix, protein extracellular matrix and extracellular space [1]. N.Wu et al performed the miRNA analysis of ovarian tissues in chickens with low and high rates of egg production using high-throughput sequencing. They found some signi cantly differentially expressed ovarian miRNAs, such as gga-miR-1744-3p, ggamiR-1655-5p, gga-miR-1734, and gga-miR-7465-3p, in the high egg-laying chickens [7]. Single nucleotide polymorphism (SNP) is a polymorphism of the DNA sequence produced by a single base mutation (base replacement, insertion or deletion, etc.), which is the most common form of variation in the genome [8]. Bhattacharya T K et al. (2019) found that the polymorphism of GnRHI and GnRHII genes was signi cantly correlated with egg production and body weight of white leghorn laying hens [9]. Li G et al. (2011) analyzed the effect of FSHR gene polymorphism on egg production in the Beijing oil chicken population. The correlation analysis showed that the polymorphic loci of g310 A > g, g181 A > T, and g159 C > T were all correlated with the number of eggs in Beijing-You chicken population of different weeks (P < 0.05) [10]. Han HX et al. (2015) studied the relationship between BMP15 gene polymorphism and reproductive traits of LaiWu Black chicken and identi ed three single nucleotide polymorphisms (SNPs) [11]. Sex hormone regulation also affects the development of sex organs and egg production. Estrogen promotes the maturation of animal reproductive organs. Estrogen cooperates with FSH to promote follicular development, inducing the emergence of LH secretion peak before ovulation and promoting ovulation. Egg-laying traits are complex traits controlled by multiple genes. However, at present, the number of genes associated with egg production remains very limited, and little is known about the links among these genes.
RNA-seq is a methodology for RNA pro ling based on next-generation sequencing that enables us to measure and compare gene expression patterns at an unprecedented resolution [12]. In recent years, to further explore the molecular mechanisms that regulate important economic traits of poultry, a large number of poultry transcriptome studies have been performed using RNA-Seq [7,[13][14][15][16]. The development of RNA sequencing has provided a powerful, highly reproducible, and cost-e cient tool for exploring the molecular mechanism of traits.
The Nandan-Yao chicken is a native chicken breed in China. It has the characteristics of coarse food resistance, strong foraging ability, delicate and delicious meat, but its performance of egg production is low. There are large differences among individuals in Nandan-Yao chicken population, which provides good material for the study of the mechanism of egg production. High egg production has an important economic value in chicken breeding, the current study aims to nd some key genes and signaling pathways in uencing egg production of Nandan-Yao chicken. In the present study, we sampled the ovarian and uterine tissues of three LEP and three HEP Nandan-Yao chickens for RNA-sEq. These data will contribute to the investigation of the regulatory mechanism of egg production.

Results
The concentration of serum sex hormone Concentrations of E2, T, FSH, LH in Nandan-Yao chicken serum are shown below (Fig. 1). The results showed that the concentration of E2 in the HEP was signi cantly higher than that in the LEP ( P < 0.01), the concentration of T in the HEP was signi cantly lower than that in the LEP (P < 0.01), and the E2/T in the HEP was signi cantly higher than that in the LEP ( P < 0.01) (Fig. 1A). The concentration of FSH and LH in the HEP was signi cantly higher than that in the LEP (P < 0.01) (Fig. 1B). These results indicated there was a better physiological function in HEP compared with LEP.

Transcriptome data
The results of transcriptome sequencing and mapping information are shown in Table 1. Total sequencing data for each sample were 77,056,834 to 105,340,082. The average GC content was 47.90%, the percentage of Q20 base was more than 95.94%, the average mapping was 91.24%.

Analysis of Differentially Expressed Genes
In the ovary, 901 genes were identi ed to be signi cantly DEGs, among of them, 472 and 429 genes were down-regulated and up-regulated, respectively ( Fig. 2A). In the uterus, we identi ed 2763 DEGs, including 1786 down-regulated genes and 977 up-regulated genes, respectively (Fig. 2B). The 299 DEGs displayed differential expression both in ovary and uterus (Fig. 2C). Hierarchical cluster analysis was used for viewing similarities between different samples. Figure 3 shows that individuals in HEP and LEP were all clustered together which illustrates the accuracy of the sample and the reliability of the differentially expressed genes obtained by our bioinformatics analysis.

GO and KEGG analysis for DEGs
Gene Ontology(GO)enrichment and Kyoto Encyclopedia of Genes and Genomes KEGG analysis was used to evaluate gene ontology and signaling pathway of DEGs (Fig. 4B). GO enrichment of DEGs in ovary showed no signi cant enrichment terms(P < 0.05). GO enrichment of DEGs in uterus revealed that there were 82 signi cantly enriched terms in the three categories(P < 0.05). The top 10 signi cantly enriched GO terms were shown in Fig. 4A, including extracellular structure organization, extracellular matrix organization, anatomical structure morphogenesis, ion homeostasis, cation homeostasis, inorganic ion homeostasis, blood vessel development, chemical homeostasis, regulation of the multicellular organismal process, cellular ion homeostasis. KEGG enrichment of DEGs in ovary showed that there was one signi cant enrichment pathway(tight junction). KEGG enrichment of DEGs in the uterus was signi cantly enriched in the 8 signaling pathways of oxidative phosphorylation, phagosome, ECM-receptor interaction, cell adhesion molecules (CAMs), focal adhesion, cardiac muscle contraction, cytokine-cytokine receptor interaction, MAPK signaling pathway (P < 0.05).

Integration of PPI network and module analysis
The DEG network interaction analysis of ovary and uterus is shown in
The results showed that the differentially expressed genes had the same expression trends in QRT-PCR and RNA-seq, which validated their accuracy.

Discussion
Although many studies have been conducted on egg production, the molecular mechanisms are still unclear. In this study, we employed RNA-seq to sequenced transcriptomes of two tissues (ovary and uterus) in three high-and three low-egg production individuals.
In this study, we identi ed 901 and 2763 DEGs in ovary and uterus, respectively. Of these, 118 and 140 DEGs were up-regulated and down-regulated in the ovary by log 2 |foldchange|> 4 while 156 and 298 DEGs were up-regulated and down-regulated in the uterus. These results suggest that these DEGs may play important roles in the regulation of egg production. For example, SPP1 was a candidate gene related to the formation of the egg and oviposition [17,18]. Studies identi ed a candidate gene BPIFB3 closely located on chromosome 20 is involved in the development of the reproductive system based on the genome-wide association study results [19]. The Glycoproteins EDIL3 regulates vesicle-mediated eggshell calci cation in avian biomineralization [20].
The DEGs detected in the ovary were mainly enriched in the tight junction (P < 0.05), while the DEGs detected in the uterus were mainly enriched in the phagosome, ECM-receptor interaction, cell adhesion molecules (CAMs), focal adhesion, cardiac muscle contraction, cytokine-cytokine receptor interaction, and the regulation of MAPK signaling pathway (P < 0.05). Tight junctions in the mucosal epithelium have essential roles as a mucosal barrier to prevent the transmembrane transporter activity, P-P-bond-hydrolysis-driven transmembrane transporter activity, ATPase activity coupled to transmembrane movement of substances, ATPase activity coupled to the movement of substances. It is possible that a large amount of energy generation and the transmembrane movement of substances is caused by egg production.

Conclusions
The current study identi ed a series of key genes and signaling pathways associated with egg production by RNA-seq and bioinformatics analysis. These key genes may regulate egg production by taking part in tight junction, phagosome, ECM-receptor interaction, cell adhesion molecules (CAMs), focal adhesion, cardiac muscle contraction, cytokine-cytokine receptor interaction, and the regulation of MAPK signaling pathway. The current ndings provide very important information for understanding the molecular mechanism of egg production and developing effective molecular markers for the molecular breed in Nandan-Yao domestic chicken.

Ethics statement
All experimental and sample collection procedures were approved by the Institutional Animal Care and Use Committee (IACUC) of the College of Animal Science and Technology of Guangxi University (Guangxi, China), with approval number GXU2018-058.

Animals and Sample Collection
A total of 300 Nandan-Yao chickens were bred in poultry experimental farm of Guangxi University, which were purchased from Guangxi Guigang Gangfeng Agriculture and Husbandry Co. Ltd. All layers were housed in individual pen of the same feeding and management condition. At the end of 50wks, the chickens were divided into high egg production (HEP) and low egg production (LEP) according to the total egg number. We selected three chickens from each HEP and LEP group as the fnal experimental birds (Table S1). The average egg number was 82 ± 1.4 and 185 ± 7.8 (Mean ± S.E.M.) for LEP and HEP chickens, respectively. Ovary and uterus were collected, then quickly frozen in liquid nitrogen and transfer into − 80℃ for permanent preservation until RNA extraction.Blood samples for testing the concentration of serum sex hormone were collected and the serum was isolated from blood by centrifuging.

Test of serum hormone concentrations
A serum sex hormone was tested by using Chicken E2, T, FSH, LH ELISA Kit (Shanghai Enzyme Biotechnology Co., Ltd.). The concentration of the sex hormone was calculated according to the standard curve. The t-test was used to test the signi cant differences between LEP and HEP.

Total RNA extraction
The total RNA was extracted from ovary and uterus using Trizol reagent (Life Technologies, USA) according to the manufacturer's instruction and modi ed a little. The quality and quantity of total RNA were tested according to the results of agarose gel electrophoresis and testing of UV-Vis Spectrophotometer Q5000 (Quawell, USA).

RNA sequencing and quality control
The cDNA libraries were constructed and sequenced on an Illumina HiSeq 2500 (Illumina, San Diego, CA, USA) in Novogene Bioinformatics Technology Co., Ltd., Beijing, China. RNA-seq was conducted following the manufacturer's standard procedures. Raw data were rst processed through trim galore [30], and clean data were obtained by removing reads containing adapter sequences, reads containing poly-N, and low-quality reads from the raw data. Furthermore, use FastQC for quality control of ltered data quality parameters of Q20, GC content, and sequence duplication levels were used for data ltering [31].

RNA-Seq analysis
Reference genome and gene model annotation les were downloaded from the genome website (ftp://ftp.ensembl.org/pub/release-95/fasta/gallus_gallus/dna/). Paired-end clean reads were aligned to the reference genome using Hisat2v2.1.0 [32,33]. Hisat2 was selected as the mapping tool because it can be faster and more precise than other splice mapping tools as well as obtain soft-splices information than other no splice mapping tools [34]. The stringtiev2.1.1 reference annotation-based transcript assembly method was used to construct and identify both known and novel transcripts from Hisat2 alignment results [35]. The differentially expressed genes analysis was performed using the DESeq2 R package (1.18.0) [36] and the criteria for screening differentially expressed genes were p-value < 0.05.

Bioinformatics analysis
All DEGs were subjected to GO and KEGG pathway enrichment analysis by the R package clusterPro ler 3.14.3 [37]. The parameters were set as minGSize = 3, p-value cutoff = 0.2, and q-value cutoff = 0.2,ont = all. The STRING [38] database is used to explore the interaction between DEGs. A con dence score > 0. 9  Validation of RNA-Seq 13 genes were selected from DEGs randomly for accuracy testing of RNA-SEq. Total RNA was isolated by Trizol as above from ovary and uterus and the rst chain of cDNA was synthesized by reverse transcription kit (Takara, Dalian, China). The Quantitative Real-time PCR(QRT-PCR) was performed according to the following program: the volume of the reaction mixture was 20 µl, with 2 µl of cDNA, 0.5 µl of primers, 10 µl of SYBR (Takara, Dalian, China), and 7 µl of RNA-free water. Then the QRT-PCR was run at 95 °C for 3 min; followed by 35 cycles of 95 °C for 10 s. β-Actin was taken as a house keep gene to correct the gene expression levels. The 2 −ΔΔCT method was used to calculate relative expression levels. The t-test was used to test the signi cant differences among genes.      Protein-protein interaction (PPI) network for DEGs in ovary and uterus. A total of 77 and 428 nodes, 175 and 2962 interaction associations in ovary and uterus, respectively. Red and blue nodes present up-regulated genes and down-regulated genes, respectively. Diamond is the nodes with the highest degree scores, while the size of the node is the expression levels of genes.

Figure 6
Three protein-protein interaction (PPI) hub network modules in the uterus. The three signi cant modules, (A) module 1 (MCODE score =18.579), (B) module 2 (MCODE score = 10), and (C) module 3 (MCODE score = 6), were constructed using PPI network of DEGs. Red and blue nodes present up-regulated genes and down-regulated genes, respectively. The size of the node is the expression levels of genes.