Transcriptome of Chicken Liver Tissues Reveals the Candidate Genes and Pathways Responsible for Adaptation into Two Different Climatic Conditions

Simple summary In this study, RNA-Seq data between liver tissues of Hanhyup chicken maintained in Korea and Kyrgyzstan was compared. We then investigated the candidate genes involved in response to environmental changes such as altitude, humidity, temperature etc. We found that 315 genes were differentially expressed genes (DEGs) in the Kyrgyz environment, out of which 174 and 141 were up- and down-regulated respectively. GO and KEGG pathway enriched between the two climatic conditions have been investigated. The identified candidate DEGs and enriched pathways may involve in the acclimatization of Hanhyup chicken to the Kyrgyz environment. Abstract RNA sequencing was used to profile the liver transcriptome of a Korean commercial chicken (Hanhyup) at two different environments (Korea and Kyrgyzstan) to investigate their role during acclimatization into different climatic conditions. Ten samples from each location were analyzed to identify candidate genes that respond to environmental changes such as altitude, humidity, temperature, etc. Sequencing reads were preprocessed, aligned with the reference genome, assembled and expressions were estimated through bioinformatics approaches. At a false discovery rate (FDR) <0.05 and fold change (FC) ≥2, we found 315 genes were DE. Out of 315 DE genes, 174 and 141 were up- and down-regulated respectively in the Kyrgyz environment. Gene ontology (GO) enrichment analysis showed that the differentially expressed genes (DEGs) were associated with energy metabolism such as pyruvate and lactate metabolic processes, and glycerol catabolic process. Similarly, KEGG pathway analysis indicated pyruvate metabolism, glycolysis/gluconeogenesis, biosynthesis, citrate cycles were differentially enriched in the Kyrgyz environment. DEGs like TSKU, VTG1, SGK, CDK2, etc. in such pathways are highly involved in the adaptation of organisms into diverse climatic conditions. Our investigation may serve as a resource for the chicken industry, especially in exporting Hanhyup chicken from Korea to other countries.


Ethics Statement and Sample Collection
Geographical location and a huge difference in climatic conditions throughout the year make Korea and Kyrgyzstan incomparable countries. There exists a vast difference in altitude, with Korea located at 250 m above mean sea level and Kyrgyzstan at 2500 m above mean sea level. The average humidity in Korea is 70% while in Kyrgyzstan it's about 40%. The Hanhyup chicken breed was reared according to commercial parameters, in NIAS campus rearing facility in South Korea. The summer season was considered for sampling of chickens. We selected two months old healthy female chickens by providing the same nutrition to all animals as guided by Hanhyup breeder a Korean Poultry breeding company. The same breed of chickens was reared under similar guidelines at the Institute of Biotechnology, National Academy of Science of Kyrgyz Republic. Liver tissues from all ten samples from each group were dissected and stored at −80 • C for RNA extraction. All animal experiments were performed by following the protocols approved by the approval committee of the National Institute of Animal Science (NIAS), Rural Development Administration (RDA), South Korea under the approval number 2018-262.

RNA Extraction, Library Construction, and Sequencing
RNA from each liver tissues was isolated using the Trizol/chloroform/isopropanol method following the manufacturer's protocol (Invitrogen TM , Carlsbad, CA, USA) was used to extract. The nano photometer (IMPLEN, CA, USA), was used to check the purity and integrity of RNA. The quality of the isolated RNA was measured by Bioanalyzer 2100 system using RNA Nano 6000 Assay kit (Agilent Technologies, CA, USA). Samples with a RIN (RNA integrity number) value >8 were considered for library preparation. The library was prepared by random fragmentation of sample, ligation of adapter, and tagmentation with the help of NEBNext Ultra™ RNA Library Prep Kit for Illumina (New England Biolabs (NEB), Ipswich, MA, USA) by following manufacturers protocol. Obtained adapter-ligated fragments were then amplified by PCR and gel purification was performed. TruSeq PE cluster kit v4-cBot-HS (Illumia) was used to generate the clusters as per the guidelines provided by the manufacturer. The library was loaded into flow cell of Illumina HiSeq 2500 platform for sequencing and generated the paired-end reads. Details of reads of each sample are shown in Supplementary Table S1. The real-time analysis (RTA) software was used for base calling of illlumina sequencer generated raw images. The Illumina package bcl2fastq was used to convert base call to FASTQ reads.

Quality Analysis and Mapping of Reads
The FASTQC program was used to check the quality of the obtained raw reads (Andrews S (2010) FastQC: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/. Accessed 2019 May 6). Good quality reads were obtained by removing the low-quality reads, poly-N reads, and phred-score ≤20. The Q20, Q30, and GC% were calculated. Trimmomatic (v0.33) was used to remove adapters from 5' and 3' end [20]. Subsequently, filtered sequences were used for reference-based genome mapping. To assemble the reads, clean reads were mapped against the reference genome of chicken (galGal6) with the help of Tophat2 (v2.1.1) mapping program [21]. Cufflink and cuffdiff (v2.2.1) suites were used to estimate the fragments per kilobase million (FPKM) value of each transcript based on expression value count [22].

Differential Gene Expression Analysis
ExDEGA (v.1.6.5) an excel based differentially expressed gene analysis tool, developed by ebiogen, South Korea, was used to estimate the DEGs. The p-values and FDR were calculated and assigned to each gene. Genes with fold change ≥2 and FDR < 0.05, were considered as DE genes.

GO, KEGG Pathway Analysis, and Network Representation
List of up-and down-regulated DEGs, were subjected to DAVID online tool for gene enrichment analysis with fold change 2, and FDR < 0.05) [23]. The Benjamini-Hochberg FDR method was used to determine the significance of enrichment. Background species: Gallus gallus, level of GO: GOTERM_BP, CC, MF_DIRECT (the categories denote DAVID defined defaults) were selected. Gene Ontology database was used to establish the association of gene products with DE genes in terms of molecular function, cellular component, and biological processes (GO: http://www.geneontology.org/) [24]. Further, to know the relation of differentially expressed genes with metabolic pathways, we used Kyoto Encyclopedia of Genes and Genomes Pathway (KEGG) database (http://www.genome.jp/kegg) pathway database [25,26]. The top ten enriched pathways related genes were used for the construction of a biological network with the help of Cytoscape (v3.2.1) and plugin ClueGO [27,28].

Raw Reads Quality, Sequencing and Mapping
With the help of the Illumina HiSeq 2500 platform, we found 7-10 million reads. GC percentage and Q20/Q30 ranged between 44.9-47.5% and 92-98% respectively, as shown in Supplementary Table S1. TopHat2 alignment tool was used to map 92-95% sequences with the reference genome of chicken (galGal6).

RNA-Seq of Liver Tissue
Cufflink suite was used to estimate the genetic expression of transcripts. The transcripts with ≥1 FPKM value were used to generate a scatter plot between Korean and Kyrgyzstan chicken liver tissues as shown in Figure 1a. The red and green dots are representations of up-and down-regulated genes respectively, whereas the black one indicate non-DEGs. Gene expression levels were estimated by using FPKM values (Roberts et al., 2011). The FPKM values were normalized based on the quantile normalization method using EdgeR (R Development Core Team, 2016) [29][30][31]. As shown in Figure 1a GO : GOTERM_BP, CC, MF_DIRECT (the categories denote DAVID defined defaults) were selected. Gene Ontology database was used to establish the association of gene products with DE genes in terms of molecular function, cellular component, and biological processes (GO: http://www.geneontology.org/) [24]. Further, to know the relation of differentially expressed genes with metabolic pathways, we used Kyoto Encyclopedia of Genes and Genomes Pathway (KEGG) database (http://www.genome.jp/kegg) pathway database [25], [26]. The top ten enriched pathways related genes were used for the construction of a biological network with the help of Cytoscape (v3.2.1) and plugin ClueGO [27,28].

Raw Reads Quality, Sequencing and Mapping
With the help of the Illumina HiSeq 2500 platform, we found 7-10 million reads. GC percentage and Q20/Q30 ranged between 44.9-47.5% and 92-98% respectively, as shown in Supplementary Table S1. TopHat2 alignment tool was used to map 92-95% sequences with the reference genome of chicken (galGal6).

RNA-Seq of Liver Tissue
Cufflink suite was used to estimate the genetic expression of transcripts. The transcripts with ≥1 FPKM value were used to generate a scatter plot between Korean and Kyrgyzstan chicken liver tissues as shown in Figure 1a Normalized read counts were used for differential expression analysis between Korean and Kyrgyz environement, by considering genes with fold change 2 and FDR < 0.05. We found 174 and Normalized read counts were used for differential expression analysis between Korean and Kyrgyz environement, by considering genes with fold change 2 and FDR < 0.05. We found 174 and 141 up-and down-regulated genes respectively. The hierarchical clustering method in MeV v4.9 tool was used for clustering genes based on Euclidean distance of z-score and average linkage [32]. Heat map along with the clustering results are shown in Figure 2. Sample wise transcriptome expression results with p-value are provided in the Supplementary File (Excel sheet 1 and 2).
To elucidate the biological functions of DEGs and their molecular interaction in the cells, we performed the biological pathway analysis. We identified 30 enriched KEGG pathways (p < 0.05). KEGG pathways clustered into seven groups such as cellular processes (blue), drug development (red), environmental information (purple), genetic information (yellow), human disease (grey), metabolism (green), and organismal systems (orange) as shown in Figure 4. Some of the highly enriched pathways include metabolic pathways, cell cycle, glycolysis/gluconeogenesis, pyruvate metabolism, PPAR signaling pathway, MAPK signaling pathway, FoxO signaling pathway as shown in Figure 4 and Supplementary File (Table S3).

Network Representation of Top Ten KEGG Pathways
We selected the top 10 KEGG pathways which are differentially enriched, for network representation by a hypergeometric test with Benjamini-Hochberg FDR correction (fold change 2, FDR < 0.05). The regulations of the genes involved in pathways are shown as a color gradient from red to blue indicates, up-or down-regulation, respectively. Pathways such as cell cycle, oocyte meiosis, progesterone mediated oocyte maturation, DNA replication, and glutathione metabolism, etc. having only upregulated genes. Whereas, pathways like Pyruvate metabolism, glycolysis/gluconeogenesis, biosynthesis of antibiotics and citrate cycle are having up-and down-regulated genes both ( Figure 5).

Discussion
In the present study, we performed a transcriptome analysis to decipher the changes in liver transcriptome response in the liver tissue of Hanhyup chicken introduced to Kyrgyzstan compared to those raised in Korea. The difference in average altitude, humidity, temperature etc. between the countries can impact the genetic expression of liver tissue of chicken during adaptation to the different climatic conditions. A total of 315 DE genes were identified in transcriptome analysis. On the basis of differential gene expression analysis of liver tissues, the top upregulated genes are TSKU, PDXK, GSTA3, SERPINB10, ANKRD22, S100A9, LYG2 and top down-regulated genes are VTG1, APOV1, PCK1, HMGCL.
Gene ontology and KEGG pathway analysis of these genes showed that 27 GO terms and 30 KEGG pathways were enriched. It also showed that the enriched genes were mainly involved in pyruvate metabolism, glycolysis/gluconeogenesis, biosynthesis of antibiotics, PPAR, glutathione metabolism, and citrate cycle pathways (p < 0.05). Pyruvate metabolism pathway plays an important role in the synthesis of glucose from lactate and dihydroxyacetone in chicken liver [33]. Many studies have reported that under any kind of stress, the chicken liver plays a significant role in maintaining homeostasis through glycolysis/Gluconeogenesis pathways [17,34,35]. Glucose is an essential nutrient for any tissue, and the liver provides glucose through gluconeogenesis [36]. In our study, enrichment of such pathway suggesting that the glucose level might have been modulated through gluconeogenesis during Hanhyup chicken acclimation to the Kyrgyz environment [37].
PPAR, a critical pathway in lipid metabolism was significantly enriched. Three genes were up-regulated (PPAR γ, Apc-CIII, LPL,) and six were down-regulated (FABP, FABP1, Thiolase B, PEPCK, CPT-1, Perilipin) in the PPAR pathway [36]. PPARγ gene acts as a regulator of adipocyte differentiation and is involved in multiple diseases such as obesity, diabetes, atherosclerosis, and cancer. The up-regulation of PPARγ gene in adipose tissue and the brain has been reported to be associated with stress-related behavior [38]. Overrepresentation of LPL, results in the deposition of lipid on the arterial wall [39]. It has been reported that the LPL gene is weakly expressed in the liver, but up-regulated when carbohydrate availability is less. It facilitates the generation of ketone into the liver and subsequently provides energy to the brain and muscle [40]. The down-regulated gene family FABP is involved in the regulation of inflammation and cellular metabolism via the PPAR pathway. FABP1 particularly activates the PPARγ, and also acts as downstream transcriptional target of PPARγ. Enrichment of PPAR pathway and involvement of DE genes in this pathway indicates that during adverse conditions such as less availability of carbohydrates, diseased conditions, and stress conditions etc., PPAR helps in adaptation.
Further, biosynthesis of antibiotics pathways, responsible for the natural synthesis of antibiotics, was also enriched. The enrichment of such pathways may indicate the defense against the stress during acclimatization in a new environment [41]. DEGs in these pathways such as UGP2, ALDOC, ASL2 and ACAA1, CAT, HADHA are up-and down-regulated respectively in Kyrgyz condition. In liver tissue, UGP2 encodes the enzyme which can activate the carbohydrate interconversions [42]. Similarly, the gene ACAA1 involved in lipid and fatty acids metabolism, may indicate that during a change of environmental conditions, these genes are actively involved in energy metabolic pathways to full fill the requirements of the liver tissue of chicken.
Likewise, liver is also a source of glutathione circulation, enrichment of glutathione metabolism in our study indicates the role of Hanhyup chicken liver during stress caused by environmental factors such as temperature and humidity [37]. The citrate cycle pathway is known for the oxidation of carbohydrates and fatty acids. It plays an important role in adapting to low levels of oxygen [43]. Désert, Colette, et al. in 2008 reported that the liver plays a central role during starvation in chicken [5]. Overall, we tried to pinpoint the various pathways and genes which are significant for adaptation of Hanhyup chicken to the Kyrgyz environment.

Conclusions
This study provides an insight into the genetic mechanisms involved in the adaptation of Hanhyup chicken, a Korean commercial chicken breed, to the Kyrgyz environment. The results of this study revealed that several metabolic pathways such as PPAR signaling, FoxO signaling, citrate cycle, gluconeogenesis, etc. are modulated during acclimation of Hanhyup chicken to Kyrgyz environment. Highly enriched genes such as PPAR γ, Apc-CIII, LPL, FABP, PCK1, UGP2, and GSTA3 are potential candidate genes for local adaptation and acclimation. However, part of it might be caused by genetic differences and breeding conditions also. A detailed functional experiment could be carried out to validate these candidate genes. This study will be beneficial for poultry industries, especially for understanding the adaptation response of chickens from one environment to another.