Original Article

Stage specific transcriptome analysis of liver tissue from a crossbred Korean Native Pig (KNP × Yorkshire)

Himansu Kumar1,http://orcid.org/0000-0002-4335-4517, Krishnamoorthy Srikanth1,†, Woncheol Park1, Kyung-Tai Lee2, Bong-Hwan Choi1, Jun-Mo Kim3, Dajeong Lim1,*, Jong-Eun Park1,*http://orcid.org/0000-0003-0718-3463
Author Information & Copyright
1Division of Animal Genomics and Bioinformatics, National Institute of Animal Science, RDA, Wanju 55365, Korea
2Division of Animal Genetics and Breeding, National Institute of Animal Science, RDA, Cheonan 31000, Korea
3Department of Animal Science and Technology, Chung-Ang University, Anseong 17546, Korea

These authors contributed equally.

*Corresponding author: Jong-Eun Park Division of Animal Genomics and Bioinformatics, National Institute of Animal Science, RDA, Wanju 55365, Korea Tel: +82-63-238-7309, E-mail: jepark0105@korea.kr Dajeong Lim Division of Animal Genomics and Bioinformatics, National Institute of Animal Science, RDA, Wanju 55365, Korea Tel: +82-63-238-7306, E-mail: lim.dj@korea.kr

© Research Institute of Veterinary Medicine, Chungbuk National University. All rights reserved.

Received: Nov 15, 2018; Revised: Dec 28, 2018; Accepted: Dec 28, 2018

Abstract

Korean Native Pig (KNP) has a uniform black coat color, excellent meat quality, white colored fat, solid fat structure and good marbling. However, its growth performance is low, while the western origin Yorkshire pig has high growth performance. To take advantage of the unique performance of the two pig breeds, we raised crossbreeds (KNP × Yorkshire to make use of the heterotic effect. We then analyzed the liver transcriptome as it plays an important role in fat metabolism. We sampled at two stages: 10 weeks and at 26 weeks. The stages were chosen to correspond to the change in feeding system. A total of 16 pigs (8 from each stage) were sampled and RNA sequencing was performed. The reads were mapped to the reference genome and differential expression analysis was performed with edgeR package. A total of 324 genes were found to be significantly differentially expressed (|log2FC| > 1 & q < 0.01), out of which 180 genes were up-regulated and 144 genes were down-regulated. Principal Component Analysis (PCA) showed that the samples clustered according to stages. Functional annotation of significant DEGs (differentially expressed genes) showed that GO terms such as DNA replication, cell division, protein phosphorylation, regulation of signal transduction by p53 class mediator, ribosome, focal adhesion, DNA helicase activity, protein kinase activity etc. were enriched. KEGG pathway analysis showed that the DEGs functioned in cell cycle, Ras signaling pathway, p53 signaling pathway, MAPK signaling pathway etc. Twenty-nine transcripts were also part of the DEGs, these were predominantly Cys2His2-like fold group (C2H2) family of zinc fingers. A protein-protein interaction (PPI) network analysis showed that there were three highly interconnected clusters, suggesting an enrichment of genes with similar biological function. This study presents the first report of liver tissue specific gene regulation in a cross-bred Korean pig.

Keywords: Korean Native Pig; Yorkshire Pig; Principal Component Analysis; edgeR; MAPK signaling pathway

Introduction

Korean Native Pig (KNP) which has a uniform black coat color and excellent meat quality is well known for full filling consumer’s demand for flavored pork meat [1, 2]. KNP produces white colored fat, solid fat structure and good marbling, however, its growth performance is low, while the western origin Yorkshire pig has high growth performance [3]. So a crossbred pig which could be high performing as well as high quality meat producing, was generated (KNP × Yorkshire) [4]. Liver is known to play a significant role in fat and lipid metabolism, although, meat quality of the pig depends upon multiple factors like intramuscular fat content, muscular pH, cholesterol level, drip loss, water holding capacity, texture, breed, feeding style, slaughtering and cooking loss [3, 5], fat content is the major deciding factor of pork meat quality[6, 7].

RNA-seq analysis of liver tissue of swine can be useful to decipher the liver fat metabolism and their significance to improve the meat quality. The pig reference genome Sscrofa10.2 still has missing annotation [8], and recently whole genome sequence using third generation sequencing; long read approach (PacBio); has been initiated for annotation of pig genome [9]. Since, RNA-seq is less prone to mapping artifacts, and can be directly mapped to the available reference genome, it can uncover the association between meat quality and liver fat metabolism. To the best of our knowledge, RNA-seq analysis has not been done so far for liver tissue of cross breed pig. In this study, we used transcriptome sequences from liver tissue of 16 crossbreed pig sample (KNP × Yorkshire), using illumina sequencing technology. Sequencing has been done at two different stages ;10 week and 26 week for comparative analysis. Until the age of 10 weeks, piglets are considered as weaners, and from 10 weeks to slaughtering time (around 26 weeks) pigs are considered as rearing pigs. We considered 10 weeks and 26 weeks, to infer gene expression changes, keeping in mind that during this time period gene expression of liver and fat metabolism will change due to the change in feeding regimen and also due to effects of age. After checking for the quality of the raw sequences we identified differentially expressed genes (DEGs). Gene ontology (GO) and KEGG pathway analysis performed to established the relation between fat metabolism and liver through transcriptome analysis. qRT-PCR was performed to validate the identified DEGs.

Materials and Methods

Ethical statement and pig rearing

A total of 16 crossbred pigs between Korean native and Yorkshire breed were used in this study. The experiment was reviewed and approved by the Institutional Animal Care and Use Committee (IACUC no. NIAS2016-848). Feed and water were supplied during the experiment. These piglets were randomly allocated into two groups (n = 8). The two groups were raised for 10 weeks (stage 1) and 26 weeks (stage 2), respectively. The animals were sacrificed at the end of the growth stage and the liver tissues were collected for RNA-seq analysis. The tissues were dissected into pieces (1 cm3), and snap frozen in liquid nitrogen and stored at –80℃ until use.

Isolation of RNAs, and RNA sequencing

Total RNAs of all the samples were extracted using TRIzol reagent (Invitrogen, United States USA), and quality checked using NanoDrop (Thermo scientific, United States USA). After performing quality control (QC), the samples were used for library preparation. The sequencing library is prepared by random fragmentation of the cDNA sample, followed by 5’ and 3’ adapter ligation. Adapter-ligated fragments were then PCR amplified and gel purified. For cluster generation, the library was loaded into a flow cell where fragments were captured on a lawn of surface-bound oligos complementary to the library adapters. Each fragment was then amplified into distinct, clonal clusters through bridge amplification. After cluster generation, the samples were sequenced in an Illumina Hiseq 4000 sequencer [10]. The sequencing was performed by Macrogen Inc., Korea. The quality of paired end reads of all 16 samples was checked by FASTQC (v0.11.4) program [11]. Low quality bases and adapter sequences were removed with Trimmomatic software [12]. Only good quality trimmed reads were considered for downstream analysis.

Quality analysis, mapping, and assembly of sequenced RNA reads

The reference genome of pig (sus scrofa) and annotation files were downloaded from NCBI (https://www. ncbi.nlm.nih.gov/genome?term=sus%20scrofa) as GCF_ 000003025.6_Ssrofa11.1_genomic.fa and indexed with Bowtie2 program (Bowtie2 v2.3.4.1) [13]. Paired end clean reads were mapped against the reference genome with HISAT2 (HIASAT2.1.0) alignment program. HISAT is a sensitive and fast alignment program for NGS reads for DNA as well as RNA [14]. Mapped sequences were taken into consideration for assembly through assembler cufflink (cufflink 2.2.1). Cufflink program was used to assemble the transcripts and find their abundance among the samples, with parameter ‘min-frags-per-transfrag = 0’ and –library-type’ and kept all parameters as default [15]. Cuffmerge, a tool of Cufflink suit was used to merged samples and get the merged gtf file for further downstream analysis.

Identification of differentially expressed genes

Differential gene expression analysis was performed with EdgeR, which uses a negative binomial model for dispersion estimate for calculating biological variation [16]. Genes with an adjusted p-value (FDR) lower than 0.001 and a two-fold change was considered as significantly differentially expressed in the pairwise comparison of the samples stage 1 and stage 2. The R package Heatmaps was used to generate the DEG heatmap.

Gene ontology and KEGG pathway

The GO database (http://www.geneontology.org/) was used for functional annotation of DEGs. The q value was calculated and GO terms with q < 0.05 were considered significantly enriched. The biological pathways and functions were analyzed on the basis of q value through DAVID, and ClusterProfiler. Similarly KEGG pathways were examined by keeping q value < 0.05 (http://www. genome.jp/kegg/).

Validation of differentially expressed genes by qRT-PCR

The quantitative real time PCR was used to validate the top five DE genes by using a reverse transcriptase kit (Takara, Japan). 1 μg of total RNA was used to synthesize as the template for quantitative PCR by the SYBR Premix Ex Taq kit (Takara, Japan). Quantitative PCR analyses were conducted using the CFX96 Touch real-time PCR system (Bio-Rad, CA, USA). The PCR conditions were as follows; reverse transcriptase at 60℃ for 90 minutes, followed by denaturing at 95℃ for 8 min and then 40 cycles of replication 94℃ for 15 s and 60℃ for 1 min; then melting curve was run from 65–95℃ at each amplification cycle for the three replicates.

Results and Discussions

Analysis of sequencing reads

Raw reads obtained from transcriptome sequencing, were processed using FastQC and trimmed by removing adaptor sequences with Trimmomatic software [11, 12]. Clean transcriptome sequences were retrieved by removing adapter’s sequences, and removing the low quality reads (Phred quality score ≤ 10). The percentage of GC, AT, Q20, Q30 of the cleaned reads were calculated and shown in Table 1 for all the sixteen samples.

jbtr-19-4-116-t1
Table 1. Sample wise quality analysis summary
Download Original Figure

Clean reads of all samples were considered for alignment with the reference genome of pig (GCF_000003025.6_ Ssrofa11.1_genomic.fna). HISAT2 was used for alignment, alignment percentage of all samples are between 89 % and 95 % as shown in Table 2. The aligned samples were counted using featurecounts (ver 1.6.3) in the Subread package [17].

jbtr-19-4-116-t2
Table 2. Read depth and alignment rate for each sample library
Download Original Figure
Differentially expressed gene analysis

Principal Component Analysis (PCA) showed that the samples clustered according to stages (Fig. 1A) and Volcano plot were generated to show the DEGs between stage 1 and stage 2 as shown in Fig. 1B. A total of 20,396 clean mapped transcripts were considered for further analysis, and a total of 324 genes were found to be significantly differentially expressed (|log2FC| > 1 & q < 0.01), out of which 180 genes were up-regulated and 144 genes were down-regulated. Hierarchical clustering of the DEGs showed that the samples clustered based on condition (stage 1 vs stage 2) in Fig. 2. In Table 3, we report DEG that plays an important role in fat and lipid metabolism for example PATZ1 gene is involved in synthesis of long-chain polyunsaturated fatty acids (LU-PUFAs) by participating with the FADS1 (fatty acid desaturase 1) and FADS2 genes [18]. Fatty acids metabolism depends upon competitive binding between the PATZ1 and SP1 (specificity protein 1) and SREBP1c (sterol regulatory element-binding protein 1c), which determines the enhancer activity and regulates the expression of FADS1 [18–20]. Similarly, ZNF202 was found to be upregulated, it acts as a transcriptional receptor that can bind to the regulatory part of multiple genes involved in lipid metabolism. ZNF202 a zinc finger protein was reported to play a role in lipid metabolism in humans [21]. Another upregulated gene HAND2 (heart and neural crest derivatives expressed 2) is significantly associated with adipose differentiation through NOTCH signaling pathway, which is involved in lipid metabolism [22].

jbtr-19-4-116-g1
Fig. 1. Comparative analysis of transcriptome profiles of porcine liver from 10 and 26 weeks. (A) Principal component analysis (PCA) of transcripts at stage1 (10 weeks) and stage 2 (26 weeks) (B) Volcano plot showing the threshold applied for identifying significant differentially expressed genes (DEGs). PCA, principal component analysis; DEGs, differentially expressed genes.
Download Original Figure
jbtr-19-4-116-g2
Fig. 2. Heatmap of differentially expressed genes based on normalized gene count data showing the change in expression between Stage 1 and Stage 2.
Download Original Figure
jbtr-19-4-116-t3
Table 3. List of differentially expressed genes
Download Original Figure

TCF4 is also responsible for lipid metabolism to reduce the fat deposition in the human body, and can interlinked with pig lipid metabolism [23]. It has been reported that the NHLH1 gene is also associated with lipid metabolism as well as neurogenesis [24]. The downregulated ZIC3 is involved in activation of zinc finger transcription factor during beginning of gastrulation [25]. Some other downregulated gene like PLAG1, that can affect the male fertility in cattle and in several species [26].

Pathway analysis

GO consists of three main components like, biological process, molecular function, and cellular component. In this study, our findings suggesting that the DE genes of porcine’s liver tissue of two stages are involved in fat and lipid metabolism. Gene enrichment analysis (at q < 0.05) of DE genes showed enrichment of multiple GO terms. We have used q value, since we have taken our gene list and compared with every single pathway that is annotated in DAVID. Genes involved in biological process (Fig. 3A) such as cell division, DNA replication, regulation of signal transduction by p53 class member, translational initiation, G1/S transition of mitotic cell cycle, nucleotide-excision repair, DNA gap filling, positive regulation of telomerase activity, protein phosphorylation, activation of GTPase activity, cellular components like kinetochore, ribosome, focal adhesion etc. and performing molecular function such as structural constituent of ribosome, Poly(A) RNA binding, DNA helicase activity, microtubule motor activity etc. were found to be enriched.

jbtr-19-4-116-g3
Fig. 3. Functional annotation analysis of DEGs, (A) Gene ontology enrichment analysis of DEGs classified under biological process, cellular component and molecular functions. Only GO’s found to be significant a q < 0.05 was plotted. (B) KEGG pathways enriched amongst the identified DEGs (q < 0.05). DEGs, differentially expressed genes.
Download Original Figure

We then investigated the significantly enriched KEGG pathways (q < 0.05) and found pathways like cell cycle, P13K-Akt signaling pathway, Gap junction, Ras signaling pathway, Phagosome, p53 signaling pathway, MAPK signaling pathway, FoxO signaling pathway etc. to be significantly enriched (Fig. 3B).

Protein-protein interaction analysis

We then performed, protein-protein interaction (PPI) network analysis with the significant DEGs identified in this study using STRING database and considered the nodes as genes and the edges are the interaction between the nodes. Three clusters involving genes that functions in cell cycle, DNA replication, and DNA repair, Fanconi Anemia pathway were found to have significant interactions, as shown in Fig. 4.

jbtr-19-4-116-g4
Fig. 4. Protein-protein interaction network analysis of DEGs. Three highly connected clusters like cell cycle, DNA replication, and DNA repair. DEGs, differentially expressed genes.
Download Original Figure
qRT-PCR analysis to validate the top five DEG genes

We have selected the top 5 DEGs like PKP1, SLC17A3, OR10P1, LRFN4, IL2 which are involved in fat metabolism of liver tissues of crossbred pig. The expression levels of five genes were analyzed in stage 1 vs. stage 2. As shown in Fig. 5A, qRT-PCR is also corroborating the findings of RNA-seq findings. The Fig. 5B is showing the correlation (r2) at value 0.9428 and ΔCT and LOG2FC of RNA-seq analysis. Many similar interpretations have been reported earlier for validation of DE genes through qRT-PCR [27–30].

jbtr-19-4-116-g5
Fig. 5. (A) Expression levels of selected genes from RNAseq analysis (black bar) and their validation by qRT-PCR (grey bar). (B) R2 is the correlation of expression levels between the two methods.
Download Original Figure

Conclusion

This study presents the first report of liver tissues gene regulation in a cross-bred Korean pig. A detailed and comprehensive investigation of mRNAs in stage specific liver tissue may serve as a good resource for further study. By analyzing the stage specific RNA-Seq data of pig liver tissues, we showed the expression profile of mRNAs and their association with fat metabolism. Through GO, pathway and PPI analysis we found genes and TFs that are involved in fat metabolism, cell cycle, DNA repair and replication and signal transduction to be significantly enriched in stage 2 pigs relative to stage 1. However, functional genetics experiments are still needed to validate the functional association of mRNAs presented in this study.

Acknowledgements

The authors are thankful to Rural Development Administration (RDA), Korea, for supporting this research. This work was carried out with the support of "Cooperative Research Program for Agriculture Science & Technology Development (Project title: Porcine epigenomic map construction and investigation of the imprinted genes, Project No. PJ01187601)" RDA, Korea and also supported RDA Fellowship Program 2018 of National Institute of Animal Science, RDA, Korea.

References

1..

Cho SH, Seong PN, Kim JH, Park BY, Kwon OS, Hah KH, Kim DH, Ahn CN. Comparison of meat quality, nutritional, and sensory properties of Korean native pigs by gender. Korean J Food Sci Anim Resour 2007;27:475-481.

2..

Oh HS, Kim HY, Yang HS, Lee JI, Joo YK, Kim CU. Comparison of meat quality characteristics between crossbreeds. Korean J Food Sci Anim Resour 2008;28:171-180.

3..

Kim GW, Kim HY. Physicochemical properties of M. longissimus dorsi of Korean native pigs. J Anim Sci Technol 2018;60:6.

4..

Lundstrom K, Enfalt AC, Tornberg E, Agerhem H. Sensory and technological meat quality in carriers and non-carriers of the RN− allele in Hampshire crosses and in purebred Yorkshire pigs. Meat Sci 1998;48:115-124.

5..

Warriss PD, Brown SN, Edwards JE, Knowles TG. Effect of lairage time on levels of stress and meat quality in pigs. Anim Sci 1998;66:255-261.

6..

O'Hea EK, Leveille GA. Significance of adipose tissue and liver as sites of fatty acid synthesis in the pig and the efficiency of utilization of various substrates for lipogenesis. J Nutr 1969;99:338-344.

7..

Enser M, Richardson R, Wood JD, Gill BP, Sheard PR. Feeding linseed to increase the n-3 PUFA of pork: fatty acid composition of muscle, adipose tissue, liver and sausages. Meat Sci 2000;55:201-212.

8..

Warr A. High quality re-assembly of the pig genome using PacBio sequencing. In: The 3rd Livestock Genomics Meeting; 2016 Sep 14-16. vol. 15.

9..

Robert C, Kapetanovic R, Beraldi D, Watson M, Archibald AL, Hume DA. Identification and annotation of conserved promoters and macrophage-expressed genes in the pig genome. BMC Genomics 2015;16:970.

10..

Caporaso JG, Lauber CL, Walters WA, Berg-Lyons D, Huntley J, Fierer N, Owens SM, Betley J, Fraser L, Bauer M, Gormley N, Gilbert JA, Smith G, Knight R. Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms. ISME J 2012;6:1621-1624.

11..

Andrews S. FastQC: a quality control tool for high throughput sequence data [Internet]. 2010 [cited 2018 Dec 24]. Available from: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

12..

Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 2014;30:2114-2210.

13..

Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods 2012;9:357.

14..

Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nature Methods 2015;12:357-360.

15..

Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc 2012;7:562-578.

16..

Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010;26:139-140.

17..

Liao Y, Smyth GK, Shi W. FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2013;30:923-930.

18..

Pan G, Ameur A, Enroth S, Bysani M, Nord H, Cavalli M, Essand M, Gyllensten U, Wadelius C. PATZ1 downregulates FADS1 by binding to rs174557 and is opposed by SP1/SREBP1c. Nucleic Acids Res 2017;45:2408-2422.

19..

Keskin N, Deniz E, Eryilmaz J, Un M, Batur T, Ersahin T, Cetin Atalay R, Sakaguchi S, Ellmeier W, Erman B. PATZ1 is a DNA damage responsive transcription factor that inhibits p53 function. Mol Cell Bio 2015;MCB:01475-01414.

20..

Reynolds LM, Howard TD, Ruczinski I, Kanchan K, Seeds MC, Mathias RA, Chilton FH. Tissue-specific impact of FADS cluster variants on FADS1 and FADS2 gene expression. PLOS ONE 2018;13:e0194610.

21..

Wagner S, Hess MA, Ormonde-Hanson P, Malandro J, Hu H, Chen M, Kehrer R, Frodsham M, Schumacher C, Beluch M , Honer C, Skolnick M, Ballinger D, Bowen BR. A broad role for the zinc finger protein ZNF202 in human lipid metabolism. J Biol Chem 2000;275:15685-15690.

22..

Keller M, Hopp L, Liu X, Wohland T, Rohde K, Cancello R, Klos M, Bacos K, Kern M, Eichelmann F, Dietrich A, Schon MR, Gartner D, Lohmann T, Dreßler M, Stumvoll M, Kovacs P, DiBlasio AM, Ling C, Binder H, Bluher M, Bottcher Y. Genome-wide DNA promoter methylation and transcriptome analysis in human adipose tissue unravels novel candidate genes for obesity. Mol Metab 2017;6:86-100.

23..

Tai CC, Ding ST. N-3 polyunsaturated fatty acids regulate lipid metabolism through several inflammation mediators: mechanisms and implications for obesity prevention. J Nutr Biochem 2010;21:357-363.

24..

Mullick A, Groulx N, Trasler D, Gros P. Nhlh1, a basic helix-loop-helix transcription factor, is very tightly linked to the mouse looptail (Lp) mutation. Mamm Genome 1995;6:700-704.

25..

Weber JR, Sokol SY. Identification of a phylogenetically conserved activin-responsive enhancer in the Zic3 gene. Mech Dev 2003;120:955-964.

26..

Utsunomiya YT, Carmo AS, Neves HH, Carvalheiro R, Matos MC, Zavarez LB, Ito PK, O'Brien AMP, Solkner J, Porto-Neto LR, Schenkel FS, McEwan J, Cole JB, da Silva MV, Van Tassell CP, Sonstegard TS, Garcia JF. Genome-wide mapping of loci explaining variance in scrotal circumference in Nellore cattle. PLOS ONE 2014;9:e88561.

27..

Nagalakshmi U, Wang Z, Waern K, Shou C, Raha D, Gerstein M, Snyder M. The transcriptional landscape of the yeast genome defined by RNA sequencing. Science 2008;320:1344-1349.

28..

Core LJ, Waterfall JJ, Lis JT. Nascent RNA sequencing reveals widespread pausing and divergent initiation at human promoters. Science 2008;322:1845-1848.

29..

Camarena L, Bruno V, Euskirchen G, Poggio S, Snyder M. Molecular mechanisms of ethanol-induced pathogenesis revealed by RNA-sequencing. PLOS Pathog 2010;6:e1000834.

30..

Sun L, Lamont SJ, Cooksey AM, McCarthy F, Tudor CO, Vijay-Shanker K, DeRita RM, Rothschild M, Ashwell C, Persia ME. Transcriptome response to heat stress in a chicken hepatocellular carcinoma cell line. Cell Stress Chaperones 2015;20:939-950.