RNA seq analyses of chicken reveals biological pathways involved in acclimation into different geographical locations

Transcriptome expression reflects genetic response in diverse conditions. In this study, RNA sequencing was utilized to profile multiple tissues such as liver, breast, caecum, and gizzard of Korean commercial chicken raised in Korea and Kyrgyzstan. We analyzed ten samples per tissue from each location to identify candidate genes which are involved in the adaptation of Korean commercial chicken to Kyrgyzstan. At false discovery rate (FDR) < 0.05 and fold change (FC) > 2, we found 315, 196, 167 and 198 genes in liver, breast, cecum, and gizzard respectively as differentially expressed between the two locations. GO enrichment analysis showed that these genes were highly enriched for cellular and metabolic processes, catalytic activity, and biological regulations. Similarly, KEGG pathways analysis indicated metabolic, PPAR signaling, FoxO, glycolysis/gluconeogenesis, biosynthesis, MAPK signaling, CAMs, citrate cycles pathways were differentially enriched. Enriched genes like TSKU, VTG1, SGK, CDK2 etc. in these pathways might be involved in acclimation of organisms into diverse climatic conditions. The qRT-PCR result also corroborated the RNA-Seq findings with R2 of 0.76, 0.80, 0.81, and 0.93 for liver, breast, caecum, and gizzard respectively. Our findings can improve the understanding of environmental acclimation process in chicken.


Scientific Reports
| (2020) 10:19288 | https://doi.org/10.1038/s41598-020-76234-8 www.nature.com/scientificreports/ environmental changes 11 . Transcriptome analysis has been widely applied to explore and identify differentially expressed genes (DEGs) involved in adaptation of chicken in various climatic stresses 12 . Measuring the gene expression also facilitates the pathways involved during environmental stresses 13 . The availability of the chicken reference genome sequence provides an opportunity for detailed analysis of stress resistant genes through their transcriptome 14 . Multiple reports confirm that even mild change (32-35 °C) in environmental temperature may negatively affect fertility in chicken [15][16][17] . We considered four tissues including breast, liver, cecum, and gizzard of Korean commercial chicken (Hanhyup-3) raised in two diverse climatic conditions in Korea and Kyrgyzstan respectively, for transcriptome analysis through RNA-Seq. We investigated GO and KEGG pathways to identify significantly enriched GO terms and  www.nature.com/scientificreports/ pathways involved in acclimatization into diverse environmental conditions. Our integrated study also provides tissue level insights into the molecular mechanisms involved in response to such different conditions.

Results
Transcriptome of four tissues; liver, breast, cecum, and gizzard of ten chickens (4 tissues /bird) from each location i.e. Korea and Kyrgyzstan have been compared. High-throughput RNA-Seq using Illumina HiSeq 2500 was performed and different quality measures of raw data have been measured through Illumina package bcl2fastq such as total bases, read count, GC (%), AT (%), Q20 (%), Q30 (%). Around 67 to 108 million raw reads per sample were generated, comprising GC percentage between 45 and 48% of generated transcripts. Complete raw data statistics of both Korean and Kyrgyzstan paired end data is provided in the supplementary file 1 and 2 respectively. After trimming adapter sequences and low quality reads through reads trimming, the quality of the reads were accessed and trimmed sequences were considered for further downstream analysis.
Mapping, differential expression, and clustering analysis of transcriptomes. Genes Fig. 2A,B respectively. Fig. 3A 18 . A circos image has been generated to show the overlap between DE genes among four tissues (Fig. 3B), blue curved link genes shows that belong to the same enriched ontology term, and the inner circle represents gene lists, where hits are arranged along the arc. Genes that hit multiple lists are colored in dark orange, and genes unique to a list are shown in light orange (Fig. 3A). The functional annotations and annotated genes were categorized into three groups such as molecular function, cellular component, and biological process. In case of liver, all DEGs were assigned to 27 GO terms among three major categories. Most significant (corrected p-value < 0.05) GO terms such as cellular process (GO: 0009987), cell (GO: 0005623), and binding (GO: 0005488) were enriched in biological process, cellular component, and molecular function respectively. Likewise in case of breast, cecum and gizzard, significant GO terms are metabolic and cellular process (GO: 0008152 and GO: 0009987), cell and organelle (GO: 0005623, and GO: 0043226), catalytic activity and binding (GO: 0003824 and GO: 0005488) were enriched in biological process, cellular component, and molecular function respectively (Fig. 3C).

Gene ontology analysis. GO enriched terms in all four tissues shown in
KEGG pathway analysis. KEGG pathway analysis has been performed to explore the molecular interaction networks within cells and putative biological functions of DE genes. As shown in Fig. 4, metabolic pathway

Discussion
Genes expressed in cells of specific tissue's represents important information for understanding the function of tissue and their physiology 19,20 . Our study cataloged the collection of expressed genes in four chicken tissues such as liver, breast, caecum, and gizzard, and reported significant changes among the tissues from each other. We identified the most abundant transcripts and the collection of genes showed tissue specific expression. GO enrichment and KEGG pathways analyses provided confidence that the dataset is of both location Korean and Kyrgyzstan is of high quality and useful for downstream analysis. Number of studies have been focused on the transcriptional changes on chicken, but little is known about the impact of environmental changes on commercially important Korean chicken breed 11,21 . www.nature.com/scientificreports/ Impact of changes on liver tissue between Korean and Kyrgyzstan environment. Liver is an important metabolic organ that plays a critical role in lipid synthesis, degradation, and transport. However, the molecular regulatory mechanisms of lipid metabolism remain unclear in chicken 22 . In response to environmental changes, 174 and 141 genes were identified as significantly up and downregulated respectively in liver tissue. Differences in expressed genes were found between the Korean and Kyrgyzstan environment, including highly expressed novel genes and pathways 7 . Significantly important changes were observed in pyruvate metabolism related genes such as LDHA, ACACA, PDHA1, PDHB, are upregulated, whereas PCK1, PCK2, LDHB and PC were down-regulated (Fig. 6)

Impact of changes in breast muscle transcriptome between Korean and Kyrgyzstan environment.
In recent years, the importance of transcriptome changes in breast muscle function of Korean chickens and of the corresponding effects on meat quality has increased. It has been reported that Korean chickens are more sensitive to heat stress during transport and at high ambient temperatures than other chickens 26  www.nature.com/scientificreports/ sues by catalyzing the biosynthesis of guanidinoacetate which is a precursor of creatine (Fig. 7). It has also been reported that the PDXK gene is involved in adipogenesis 27 . GO terms such as glycolytic pathway (PKM, BPI, TPI1, and PFKL etc.), gluconeogenesis (GOT2, GPI, TPI1, PGAM1, ENO3, and PGK1 etc.), myosin complex (MYH1E-F, MYL1 etc.). Likewise KEGG pathways investigation showed that PPAR signaling pathway, ribosome, steroid biosynthesis, metabolic pathways, biosynthesis of antibiotics, pentose phosphate pathways, glycolysis/gluconeogenesis, carbon metabolism, pentose phosphate pathways were enriched in Kyrgyzstan. APOC3 gene is downregulated in PPAR pathway, indicates that the level of triglycerides decreasing by inhibiting the hydrolysis of lipoprotein in plasma 28 .

Impact of changes on caecum between Korean and Kyrgyzstan environment. The cecum is
an important part of the animal digestive system. Recent transcriptome studies have been focused on digestive system of chicken, but the molecular mechanisms underlying chicken cecum metabolism remain poorly understood. To identify genes related to cecal metabolism and to reveal molecular regulation mechanisms, we sequenced and compared the transcriptomes of cecum. In response to environmental challenges, a total of 167 genes were identified as significantly differentiated in caecum tissue.  (CDK1, CDK2, CCNB2, CCNB3, CYCS, GTSE1), progesterone mediated oocyte maturation (CCNB2-3, CDK1-2, BUB1, MAD2L1), cell cycle, glutathione metabolism (PTTG1, BUB1B, CCNA2, MCM2, CCNB2-3, PLK1), pyrimidine metabolism (UCK2, DCTD, TK1, DGUOK, RRM1). All mentioned genes were downregulated as shown in Fig. 8.

Conclusion
Primary objective of this study was to explore the molecular mechanism involved in tissues such as liver, breast, cecum, and gizzard when chicken encounters the environmental changes. In conclusion, our results give new insight to understand the changes in transcriptome that may occur during acclimatization into different locations. GO terms analyses showed that the metabolic and cellular process, cell and organelle, catalytic activity and binding were significantly enriched in almost all tissues during acclimatization into Kyrgyzstan. Similarly, KEGG pathway investigation showed that pathways such as PPAR, MAPK, FoxO, and CAM are significantly enriched in Kyrgyzstan environment. Our findings can be a useful resource for the transcriptomic investigations of adaptation efficiency-related genes in chicken and may provide significant clue for understanding the molecular genetic mechanisms in other chicken breeds.

Materials and methods
Preparation of tissues. We considered 2 months old female Hanhyup-3 chickens obtained from Hanhyup breeder, a poultry breeding company in South Korea, and grown in rearing facility of NIAS South Korea at 28 °C, humidity 80%, and pressure 1008 mbar. Four tissues such as liver, breast, cecum, and gizzard were considered for sampling. Similarly, same breeds of chicken were grown in rearing facility of Institute of Biotechnology, National Academy of Science of Kyrgyz Republic at temperature 26 °C, humidity 40%, and pressure 1003 mbar. 10 samples for each tissue from each location were considered for the experiment. Dissected tissues samples were stored at − 80 °C for further use.
RNA isolation and quality assessment. The total RNA was extracted by using Trizol/chloroform/isopropanol method (InvitrogenTM, Carlsbad, CA, USA) from all four samples liver, breast, cecum, and gizzard. The purity and integrity of RNA was estimated by the nano photometer (IMPLEN, CA, USA). RNA quality was checked by Bioanalyzer 2100 system using RNA Nano 6000 Assay kit (Agilent Technologies, CA, USA).  Reads mapping, quantification of gene expression level, and differential expression analysis. Qualified reads were mapped through Tophat2 (v2.1.1) software with default parameters to the chicken reference genome (GalGal5) [32][33][34] . The mapped BAM file has been used to estimate the gene expression through fragments per kilo base of transcript per million fragments mapped (FPKM) through cufflinks (version 2.2.1) 35 . Differential expression (DE) analysis of Korean and Kyrgyzstan groups (10 biological replicates from each location) was performed using the DESeq2 R package (version: 1.28.1). The P values adjustment had been done using the Benjamini and Hochberg's approach for controlling the FDR. Genes with a P value < 0.05 and fold change ≥ 2 were used as threshold for DE by DESeq2 36 .
GO and KEGG enrichment analysis of differentially expressed genes. GO  Quantitative reverse-transcription-PCR (qRT-PCR). We carried qRT-PCR for top ten DE genes obtained from RNA-seq for the purpose of validation. Primers were designed via Primer3 (https ://bioin fo.ut. ee/prime r3-0.4.0/prime r3/input .htm) 39 . The mRNA levels of the DEGs were normalized by the housekeeping genes β-Actin and GAPDH (Glyceraldehyde 3-phosphate dehydrogenase). The relative gene expression values were calculated using the 2 − ΔΔCt method. Finally, the correlations between RNA-Seq of top 5 up and down regulated genes and the mRNA expression level from qRT-PCR were estimated.
Ethics statement. All experiments were performed in accordance with relevant guidelines and regulations.
The experiments were also conducted according to norms approved by the National Institute of Animal Science www.nature.com/scientificreports/