Journal of Animal Science and Technology
BioMed Central
Research

Whole-transcriptome analyses of the Sapsaree, a Korean natural monument, before and after exercise-induced stress

Ji-Eun Kim1, Junkyung Choe1, Jeong Hee Lee1, Woong Bom Kim1, Whan Cho1, Ji Hong Ha2, Ki Jin Kwon2, Kook Il Han3, Sung-Hwan Jo1
1SEEDERS Inc, Daejeon, 34015 Republic of Korea
2The Korean Sapsaree Foundation, Gyeongsan, Gyeongbuk 38412 Republic of Korea
3Department of Animal training and event, Daekyeung University, Gyeongsan, Gyeongbuk 38547 Republic of Korea

© Kim et al. 2016. Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Received: Sep 4, 2015 ; Accepted: Mar 29, 2016

Published Online: Apr 15, 2016

Abstract

Background

The Sapsaree (Canis familiaris) is a Korean native dog that is very friendly, protective, and loyal to its owner, and is registered as a natural monument in Korea (number: 368). To investigate large-scale gene expression profiles and identify the genes related to exercise-induced stress in the Sapsaree, we performed whole-transcriptome RNA sequencing and analyzed gene expression patterns before and after exercise performance.

Results

We identified 525 differentially expressed genes in ten dogs before and after exercise. Gene Ontology classification and KEGG pathway analysis revealed that the genes were mainly involved in metabolic processes, such as programmed cell death, protein metabolic process, phosphatidylinositol signaling system, and cation binding in cytoplasm. The ten Sapsarees could be divided into two groups based on the gene expression patterns before and after exercise. The two groups were significantly different in terms of their basic body type (p ≤ 0.05). Seven representative genes with significantly different expression patterns before and after exercise between the two groups were chosen and characterized.

Conclusions

Body type had a significant effect on the patterns of differential gene expression induced by exercise. Whole-transcriptome sequencing is a useful method for investigating the biological characteristics of the Sapsaree and the large-scale genomic differences of canines in general.

Electronic supplementary material

The online version of this article (doi:10.1186/s40781-016-0097-1) contains supplementary material, which is available to authorized users.

Keywords: NGS; Bioinformatics; RNA-Seq; Transcriptome; Exercise; Physical stress; Sapsaree

Background

The Sapsaree (Canis familiaris) is a native Korean dog that is distributed throughout the Korean peninsula, and is very friendly, protective, and loyal to its owner. The Sapsaree population size decreased dramatically during the Korean War in the 1950s, and the breed was considered endangered. A program of systematic mating and reproduction to preserve the Sapsaree from extinction and maintain a pure pedigree generated the current population of about 4,000 individuals, including the 500 dogs now living at the Sapsaree Breeding Research Institute in Gyeongsan, Gyeongbuk province [1]. The Sapsaree was registered as a Korean Natural Monument (number: 368) in 1992. Although recent studies have shed light on the origin and various morphological and behavioral traits of the Sapsaree, the genetic backgrounds, abundant genetic polymorphisms, and novel genes of the Sapsaree are still not completely understood [2].

RNA sequencing (RNA-Seq) is one of the most useful next-generation sequencing tools for investigating the landscape of a whole transcriptome. Using RNA-Seq, researchers have identified differentially expressed and novel genes, unraveled the expression profiles underlying phenotypic changes, and discovered unannotated, transcriptionally active regions that cannot be detected by conventional gene prediction [3]. Furthermore, the recent use of RNA-Seq has captured the scale and complexity of organ-specific and tissue-specific transcriptomes, making RNA-Seq the technique of choice for investigating gene expression during complex phenomena such as stress [4].

Exercise is usually recognized as a stress factor, as is any environmental change that activates or pressures cells and tissues. And the level of response to physical exercise is different from individuals to individuals. So, it is needed to research the biological responses under physical stress. Several RNA-Seq studies of exercise-induced stress have been performed in equines [3, 4]. There has been little progress, however, in canine-based RNA-Seq studies, which have been limited to certain diseases. Therefore, we performed a large-scale analysis of whole-transcriptome data to investigate the gene expression levels before and after exercise in the Sapsaree. This study will provide the basic approach for identified the biologic characteristics and breeding the working dogs.

Methods

Morphological traits and blood sample collection

The dogs were handled in accordance with Article 23 ‘Experiments with Animal’ of Korea’s Animal Protection Law, 2015, and the Korean Sapsaree Foundation cooperated and approved all animal care procedures prior to the initiation of the experiment. Ten Sapsarees housed in controlled environmental conditions were used in the experiment. The basic morphological traits of each Sapsaree including body height (BH), body length (BL), depth of chest (DC), body weight (W), and hair color were measured before the experiment. Then, the dogs were exercised one-on-one by trainers for 1 h in 5 min intervals, with the trainers leading the dogs in quick step 10 times around a course on a 20 m × 20 m square field located at the Korean Sapsaree Foundation facility in Gyeongsan. The exercise course included hurdling (40 cm, 70 cm, and 80 cm), high jumping (50 cm), A-shape hurdling, bridge-shape hurdling, and seesaw hurdling.

Blood samples were drawn under a veterinarian’s supervision from a cephalic vein before and immediately after exercise, resulting in a total of 20 samples. Immediately after collection, portions of the blood samples were dispensed into serum separator tubes (SST) (1.5 ml in each tube [Greiner Bio One, Kremsmuenster, Austria]) to measure hormone levels, and into PAXgene blood RNA tubes (2.5 ml in each tube, [PreAnalytiX, Hombrechtikon, Switzerland]) for RNA sequencing according to the respective manufacturers’ protocols [5, 6]. The blood for serum was allowed to clot at room temperature before centrifugation. The serum was then separated and stored frozen at -70 °C until completion of the case enrollment and sample collection.

Measurement the stress indicator in blood

Four substances related to physical stress were chosen based on previous studies: cortisol, aspartate aminotransferase (AST), creatine kinase (CK), and creatinine. The serum CK, AST, creatinine levels are commonly used to muscle damage indicators and renal function [710]. And a recent study indicated that physical stress resulted in immediate increase in the plasma concentrations of cortisol [11]. The level of target substances in serum from the upper side of the SST tubes were measured using a BS-400 chemistry analyzer (Mindray, Shenzhen, China) and an Immulite 1000 Immunoassay System (Siemens, New York, USA) [12, 13].

RNA sequencing and bioinformatic analysis

Total RNA was extracted from the PAXgene tubes according to the manufacturer’s instructions [6]. Starting with the total RNA, mRNA was purified using poly (A) selection or rRNA depletion, converted into double-stranded cDNA, and amplified by PCR. To check the RNA quality, all the RNA samples were examined for RNA Integrity Number, and 28S to 18S rRNA value using Bioanalyzer. Next, the construction of library were used with the Illumina TruSeq RNA Sample Preparation Kit v2 (catalog #RS-122-2001, Illumina, San Diego, CA) following the manufacturer’s instructions [14]. And the library was quantified using the KAPA library quantification kit (Kapa Biosystems KK4854) following the manufacturer’s instructions [15]. The final individual libraries were sequenced using the Illumina Hiseq2000 platform, which created 100 bp paired-end (PE) RNA-sequencing reads.

To collect high-quality transcriptome data, we filtered the sequencing data by phred score (Q ≥ 20) and minimum length (≥25 bp) using the SolexaQA software [16]. The filtered reads were mapped to 48,370 reference mRNAs from Canis lupus familiaris using the bowtie2 software (mismatches ≤ 2) [17, 18]. The number of mapped reads for each mRNA was counted and then normalized using the DESeq packages in R [19]. Differentially expressed genes (DEGs) were selected by over 100 mapped read counts, a ≥ twofold change in reads coverage and a binomial test with a false discovery rate (FDR) ≤ 0.01 at the first. And then, the final DEGs were identified by a ≥ twofold change in reads coverage, a binomial test with a false discovery rate (FDR) ≤ 0.01, and a read count ≥ 1,000 either before or after exercise. The FDR was applied to identify the threshold p-value for multiple tests and was calculated using DESeq. Correlation analysis and hierarchical clustering was performed to group the genes according to patterns of expression using the AMAP library in R [20].

And we tested the stability of the gene expression levels before and after the exercise of housekeeping genes, such as HNRNPH1, GAPDH, RPL8, TAF4B, and TAF1 with t-test method supporting the statistical significance [21].

Functional enrichment analysis

Functional enrichment analyses were carried out using the Gene Ontology (GO) database and including all three GO categories (biological processes, cellular components, and molecular functions), providing a structured and controlled vocabulary to describe the gene products [22]. We also used the KEGG database to identify the biological mechanisms and metabolic pathways associated with the differentially expressed genes corresponding their enzyme commission numbers [23]. DAVID is a web-accessible annotation system (https://david.ncifcrf.gov/home.jsp) that provides a comprehensive set of functional annotation tools for investigators for understanding the biological meanings behind large lists of genes [24]. We used DAVID to analyze the clusters of differentially expressed genes annotated by the Entrez gene IDs of the genes with counts ≥ 2 and FDR ≤ 0.1 of each GO and KEGG term.

Statistical analysis

To identify correlations between the gene expression pattern and body weight, we carried out t-tests comparing W/BH, W/BL, and W/DC values between two groups of dogs that showed different gene expression patterns. To control for group differences in baseline weight, the weight was divided by the BH, BL, and DC. A p-value ≤ 0.05 was used as the cutoff for significance in all analyses using the t-test function of R.

Results and discussion

Morphological traits and blood sample collection

The basic morphological traits of the 10 Sapsarees measured before the experiment are shown in Table 1. All the dogs were males, 13 to 60 months of age. The BH measured from the ground to the top of the withers ranged from 49 cm to 62 cm. The BL measured from the point of the shoulder to the rear point of the croup ranged from 58 cm to 70 cm. The DC measured from the elbow to the top of the withers ranged from 21 cm to 28 cm. The W measured by weighing balance ranged from 18.7 kg to 30 kg [25].

Table 1. Summary of the morphological traits of 10 Sapsarees
Group Name Birth Month of agea Sex Hair colorb BHc (cm) BLd (cm) DCe (cm) Wf (kg) W/BH W/BL W/CL
I Cheongbaek 2011.02.15 46 Male BT 58 65 27 22 0.38 0.34 0.81
I Rookie 2010.03.23 57 Male W 62 67 25 24.5 0.40 0.37 0.98
I Chaeum 2010.10.24 50 Male Y 50 60 21 15.4 0.31 0.26 0.73
I Hwangryong 2011.11.08 37 Male DY 60 69 25 20 0.33 0.29 0.80
I Pyeonggang 2013.11.19 13 Male BT 49 59 21 18.7 0.38 0.32 0.89
Average - - - - 55.8 64 23.8 20.12 0.36 0.316 0.842
II Tong 2010.08.15 52 Male BT 58 66 26 31 0.53 0.47 1.19
II Huimang 2009.12.11 60 Male W 58 68 28 24 0.41 0.35 0.86
II Pyeongtan 2013.11.17 13 Male Y 54 63 22 25.5 0.47 0.40 1.16
II Hwangdol 2011.09.29 39 Male SY 62 70 27 30 0.48 0.43 1.11
II Bongsik 2010.07.22 53 Male DY 53 58 23 19 0.36 0.33 0.83
- Average - - - - 57 65 25.2 25.9 0.45 0.396 1.03

aMonths of age measured from birth to the date of exercise

bHair color : BT (Black & Tan), W (White), Y (Yellow), DY (Dark Yellow), SY (Strong Yellow)

cBH Body Height, dBL Body Length, eDC Depth of Chest, fW Body Weight

Download Excel Table

Physical stress indicators in serum

Compared with those before exercise, the concentrations of AST, CK, and creatinine were slightly increased after exercise, but the increases were not significant (Additional file 1: Tables S1 and S2). Cortisol, a key hormone from the adrenal glands, was significantly elevated after exercise in all individuals except for Cheongbaek (Additional file 2: Figure S1, Additional file 1: Table S2). Hence, it seems that cortisol can be used as a marker of exercise-induced stress.

RNA sequencing and Bioinformatic analysis

After the RNA quality check, the RNA sequencing generated over 75 Gbp (about 3.8 Gbp per sample) of data consisting of 100 bp paired-end reads (Table 2). Trimming resulted in reads with a mean length of 86.58 bp across all samples and a total combined length of about 57 Gbp, which was 74.9 % of the raw sequence.

Table 2. Summary of raw sequencing reads
Name Exercise Num. of reads (ea) Avg. length (bp) Total length (bp)
Cheongbaek before 31533452 100 3153345200
after 38271980 100 3827198000
Rookie before 41724150 100 4172415000
after 38515768 100 3851576800
Chaeum before 34658490 100 3465849000
after 33784972 100 3378497200
Hwangryong before 36358784 100 3635878400
after 41631454 100 4163145400
Pyeonggang before 38708814 100 3870881400
after 43781622 100 4378162200
Tong before 40314536 100 4031453600
after 39786602 100 3978660200
Huimang before 35384558 100 3538455800
after 32013742 100 3201374200
Pyeongtan before 41878502 100 4187850200
after 40729858 100 4072985800
Hwangdol before 39336568 100 3933656800
after 38413046 100 3841304600
Bongsik before 33239988 100 3323998800
after 37272096 100 3727209600
Total - 75733982 100 75733898200 (100 %)
Download Excel Table

Using bowtie2, 80.76 % of the filtered reads were successfully mapped to the current dog reference genes (Canis lupus familiaris, 48,370 mRNAs) [18, 26]. A novel bioinformatics pipeline for processing large amounts of transcriptome sequences was built. We calculated the expression levels of all the genes with mapped reads from the 10 individuals before and after exercise. By comparing the coverage before and after exercise, we identified 2,549 DEGs. The numbers of up-regulated and down-regulated genes were different in each individual. Pyeonggang, Huimang, and Pyeongtan had fewer than 10 DEGs, but each of the other dogs had more than 100 genes (Additional file 3: Figure S2). After filtering out genes with low levels of expression, 525 genes were grouped into two clusters (C1 and C2) of 276 and 249 genes, respectively, depending on their expression pattern (Fig. 1 and Additional file 4).

jast-58-0-17-g1
Fig. 1. HeatMap showing hierarchical clustering of differentially expressed genes regulated by exercise. The log2Ratio for each significantly differentially expressed gene was used. Each column represents a Sapsaree individual, and each row represents a differentially expressed gene. Expression differences are shown in different colors; red indicates up-regulation after exercise, and green indicates down-regulation after exercise. The 525 genes were grouped into two clusters (C1 and C2). The dogs were divided into two groups. Group I included Cheongbaek, Rookie, Chaeum, Hwangryong, and Pyeonggang. Group II included Tong, Huimang, Pyeongtan, Hwangdol, and Bongsik
Download Original Figure

And the reliability of the gene expression analysis based on RNA-Seq was confirmed by t-test which supported the statistical significance of the stability of housekeeping genes expression before and after the exercise (Additional file 1: Table S3).

Functional enrichment analysis

A total of 26 unique genes among the 525 differentially expressed genes were assigned to 11 functional groups based on GO assignments (Additional file 1: Table S4). The genes were involved in biological processes such as programmed cell death, negative regulation of cell motion, regulation of cellular component biogenesis, protein metabolic process, and others. The cellular components linked to the genes were the intracellular parts, such as cytoplasm. The molecular function assignments were mainly to the catalytic and binding activities, such as cation binding.

A further functional classification of all the differentially expressed genes was performed using the KEGG database. A total of 84 unique genes among the 525 differentially expressed genes were assigned to 30 metabolic pathway terms, including phosphatidylinositol signaling system, ribosome, proteasome, oxidative phosphorylation, and others (Additional file 1: Table S5). The functional annotation analyses indicated that the genes regulated under exercise-induced stress were mainly related to the metabolite pathways.

Correlation between gene expression pattern and body type

Based on the gene expression patterns induced by exercise, the 10 Sapsaree individuals were divided into two groups (Groups I and II; Fig. 1). Group I included Cheongbaek, Rookie, Chaeum, Hwangryong, and Pyeonggang, while Group II included Tong, Huimang, Pyeongtan, Hwangdol, and Bongsik. The differentially expressed genes formed two clusters. Cluster 1 (C1) was down-regulated after exercise in Group I but up-regulated after exercise in Group II. Cluster 2 (C2) was up-regulated after exercise in Group I but down-regulated after exercise in Group II.

We examined the cortisol change after exercise, age, and body type for differences between the two groups. Only the body type showed a significant difference between Group I and Group II. Using W/BH, W/BL, and W/DC to characterize the overall body type, we found that Group I had a heavier body type than Group II (Table 1). In addition, W/BH, W/BL, and W/DC were each significantly different (p ≤ 0.05) between the two groups based on independent t-tests (Table 3). Those findings suggest differences in gene regulation between heavier dogs and lighter dogs under exercise stress.

Table 3. Differences in body type between groups I and II
Mean of group Ib Mean of group IIc td dfe p-valuef
W/BHa 0.357 0.452 -3.386 8.124 0.0093
W/BLa 0.312 0.397 -3.211 9.413 0.0100
W/DCa 0.823 1.018 -2.598 8.364 0.0306

aW Body Weight, BH Body Height, BL Body Length, DC Depth of Chest

bMean of Group I: average value of Group I

cMean of Group II: average value of Group II

dt : t-value, edf : the degree of freedom value, fp-value : significance

Download Excel Table

Identification of candidate marker genes for different response patterns to exercise-induced stress

We identified 525 genes with differential expression before and after exercise. We manually curated seven genes (PTPRC, LOC102157092, RPL18, S100A8, LOC612054, LOC102151356, EIF1) that had distinctly different expression patterns before and after exercise in Groups I and II (Fig. 2). For example, PTPRC (Entrez ID : 490255) is described as ‘protein tyrosine phosphatase, receptor type, C’. Some genes had uncharacterized functions. Figure 2 shows the log2 fold changes of the seven selected genes after exercise for each of the dogs. There were clear differences in the changes in expression levels of the seven selected genes after exercise between the two groups of dogs. The seven selected genes could be used as markers of exercise stress and for grouping dogs according to body type.

jast-58-0-17-g2
Fig. 2. Representative genes that were differentially expressed between the two groups of dogs. The x-axis represents the seven selected genes, and the y-axis represents the log2 fold change (2FC) value of each gene. The 2FC values for each individual are represented by a colored bar. a shows two genes with negative 2FC values in Group I and positive 2FC values in Group II. b shows five genes with positive 2FC values in Group I and negative 2FC values in Group II. The gene identifiers are shown at the top of each histogram: PTPRC (490255, protein tyrosine phosphatase, receptor type c), LOC102157092 (102157092, complement receptor type 1-like), RPL18 (476422, ribosomal protein L18), S100A8 (490461, S100 calcium binding protein A8), LOC612054 (612054, uncharacterized LOC612054), LOC102151356 (102151356, uncharacterized LOC102151356), and EIF1 (403674, eukaryotic translation initiation factor 1)
Download Original Figure

Conclusions

This study provides the Sapsaree whole-transcriptome sequence data and identifies genes that are differentially expressed before and after exercise in the Sapsaree, which could be used as markers of exercise stress. The pattern of changes in global gene expression induced by exercise was different depending on the body type. RNA sequencing and gene expression analysis can be useful for grouping the Sapsarees by body type and response to exercise. This study provides a basis for future research investigating the biologic characteristics of the Sapsaree and the large-scale genomic differences of canines in general.

Appendices

Additional files
Additional file 1: Table S1
Concentration of four substances in blood at before and after exercise. Table S2. The test for statistical significance of the four substances in serum level of before and after exercise. Table S3. The test for statistical significance in the stability of housekeeping genes expression before and after exercise. Table S4. GO functional classification of DEGs in cluster. Table S5. KEGG metabolic pathway annotation of DEGs in cluster. (DOCX 34 kb)
jast-58-0-17-suppl1.docx
Additional file 2: Figure S1
Cortisol concentration change in blood before and after exercise. (PPTX 50 kb)
jast-58-0-17-suppl2.pptx
Additional file 3: Figure S2
MAplot showing the expression pattern and the number of genes in each Sapsaree before and after exercise. Expression pattern is shown in different colors dots; red dots indicates up-regulation and green dots indicates down-regulation. And right-top side showing the number of up- regulated and down-regulated genes. Pyeonggang, Huimang, and Pyeongtan had only a few differential gene expressed before and after exercise while others had many more. (PPTX 491 kb)
jast-58-0-17-suppl3.pptx
Additional file 4:
Gene expression value of the 525 DEGs. (XLSX 163 kb)
jast-58-0-17-suppl4.xlsx

Abbreviations

BH

body height

BL

body length

DC

depth of chest

W

body weight

RNA-Seq

RNA sequencing

DAVID

database for annotation, visualization, and integrated discovery

KEGG

kyoto encyclopedia of genes and genomes

Notes

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

KJE, CJK, and LJH wrote the manuscript. The manuscript was revised by JSH. KJE and CJK performed the bioinformatic and statistical analyses. JSH supervised and coordinated the study design and the data analysis. KWB, CW, HKI, HJH, and KKJ supported the exercise experiment of dogs and collection of blood samples. All authors read and approved the final manuscript.

Acknowledgements

This work was supported by ‘biological study for protection and preservation of the Sapsal Dog, Natural monument no. 368 (NRICH-1507-B12F)’, funded by the National Research Institute of Cultural Heritage, Republic of Korea.

References

1.

Ha JH, Alam M, Lee DH, Kim JJ. Whole Genome Association Study to Detect Single Nucleotide Polymorphisms for Behavior in Sapsaree Dog (Canis familiaris). Asian-Australas J Anim Sci. 2015; 28(7):936-42.

2.

Kim KS, Jeong HW, Park CK, Ha JH. Suitability of AFLP markers for the study of genetic relationships among Korean native dogs. Genes Genet Syst. 2001; 76(4):243-50.

3.

Park KD, Park J, Ko J, Kim BC, Kim HS, Ahn K, et al. Whole transcriptome analyses of six thoroughbred horses before and after exercise using RNA-Seq. BMC Genomics. 2012; 13:473.

4.

Capomaccio S, Vitulo N, Verini-Supplizi A, Barcaccia G, Albiero A, D’Angelo M, et al. RNA sequencing of the exercise transcriptome in equine athletes. PLoS One. 2013; 8(12):e83504.

5.

Greiner bio-one (http://www.greinerbioone.com/). Accessed 4 Apr 2016.

7.

Clarkson PM, Kearns AK, Rouzier P, Rubin R, Thompson PD. Serum creatine kinase levels and renal function measures in exertional muscle damage. Med Sci Sports Exerc. 2006; 38(4):623-7.

8.

Baird MF, Graham SM, Baker JS, Bickerstaff GF. Creatine-kinase- and exercise-related muscle damage implications for muscle performance and recovery. J Nutr Metab. 2012; 2012:960363.

9.

Piccione G, Giannetto C, Fazio F, Casella S, Caola G. A comparison of daily rhythm of creatinine and creatine kinase in the sedentary and athlete horse. J Equine Vet Sci. 2009; 29(7):575-80.

10.

Nathwani RA, Pais S, Reynolds TB, Kaplowitz N. Serum alanine aminotransferase in skeletal muscle diseases. Hepatology. 2005; 41(2):380-2.

11.

Mastorakos G, Pavlatou M, Diamanti-Kandarakis E, Chrousos GP. Exercise and the stress system. Hormones (Athens). 2005; 4(2):73-89.

12.

Mindray (http://www.mindray.com/en/products/34.html). Accessed 4 Apr 2016.

13.

SIMENS (http://www.healthcare.siemens.com). Accessed 4 Apr 2016.

16.

Cox MP, Peterson DA, Biggs PJ. SolexaQA: At-a-glance quality assessment of Illumina second-generation sequencing data. BMC Bioinformatics. 2010; 11:485.

17.

Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009; 10(3):R25.

18.

Lindblad-Toh K, Wade CM, Mikkelsen TS, Karlsson EK, Jaffe DB, Kamal M, et al. Genome sequence, comparative analysis and haplotype structure of the domestic dog. Nature. 2005; 438:803-19.

19.

Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010; 11(10):R106.

20.

Antoine Lucas. AMAP: Another Multidimensional Analysis Package. R package version 0.8-12. 2014. (http://CRAN.R-project.org/package=amap). Accessed 4 Apr 2016.

21.

Brinkhof B, Spee B, Rothuizen J, Penning LC. Development and evaluation of canine reference genes for accurate quantification of gene expression. Anal Biochem. 2006; 356(1):36-43.

22.

Gene Ontology Consortium. The Gene Ontology (GO) database and informatics resource. Nucleic Acids Res. 2004; 32(Database issue):D258-61.

23.

Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000; 28(1):27-30.

24.

da Huang W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009; 4(1):44-57.

25.

Han KI, Alam M, Lee YM, Lee DH, Ha JH, Kim JJ. A study on morphology and behavior of the Sapsaree: A Korean native dog (Canis familiaris). J Anim Sci Technol. 2010; 52:481-90.

26.

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