Evolutionary and Functional Analysis of Korean Native Pig Using Single Nucleotide Polymorphisms.

Time and cost-effective production of next-generation sequencing data has enabled the performance of population-scale comparative and evolutionary studies for various species, which are essential for obtaining the comprehensive insight into molecular mechanisms underlying species- or breed-specific traits. In this study, the evolutionary and functional analysis of Korean native pig (KNP) was performed using single nucleotide polymorphism (SNP) data by comparative and population genomic approaches with six different mammalian species and five pig breeds. We examined the evolutionary history of KNP SNPs, and the specific genes of KNP based on the uniqueness of non-synonymous SNPs among the used species and pig breeds. We discovered the evolutionary trajectory of KNP SNPs within the used mammalian species as well as pig breeds. We also found olfaction-associated functions that have been characterized and diversified during evolution, and quantitative trait loci associated with the unique traits of KNP. Our study provides new insight into the evolution of KNP and serves as a good example for a better understanding of domestic animals in terms of evolution and domestication using the combined approaches of comparative and population genomics.

INTRODUCTION

Advances in next-generation sequencing (NGS) technologies have significantly reduced the time and cost of producing population-scale genome sequences of various species. Many population-scale genome projects, such as the 1,000 genomes project (Auton et al., 2015), 1,000 bull genomes project (Daetwyler et al., 2014) and 100,000 genomes project (Turnbull et al., 2018), have been launched to understand genetic diversity in a population as well as their phenotypic consequences. Genome projects of domestic animals have been of particular interest due to their importance as a food resource for humans. For example, variants associated with milk production and curly coat were identified by the sequencing of 234 cattle genomes in the 1,000 bull genomes project (Daetwyler et al., 2014), and genes related to lipid metabolism were identified in the sheep genome project (Naval-Sanchez et al., 2018). Evolutionary studies of the various distinctive traits of domesticated mammals have also been performed to acquire comprehensive insights into molecular mechanisms underlying species-specific traits (Chen et al., 2019; Gordon et al., 2007; Lee et al., 2017; Vieira and Rozas, 2011).

Pigs are a member of the order Artiodactyla, which represents even-toed ungulates such as pigs, deer, giraffes, sheep, and goats. They diverged from rodents and primates approximately 94 million years ago (Groenen et al., 2012), and were first domesticated from wild boar populations in Europe and Asia 10,000 years ago (Giuffra et al., 2000). Since then, pigs have played an important role as a major source of meat for humans, and subsequent stages of selective breeding have significantly improved economically important traits, such as meat quality, backfat depth, and growth rate. A number of population-scale studies have been performed on various pig breeds to understand genomic features, such as genes, regulatory elements and variants, and the underlying molecular mechanisms affecting such beneficial traits (Davoli and Braglia, 2007; Sachs and Galli, 2009; Zhu et al., 2017).

The Korean native pig (KNP) is a domestic pig breed in Korea that originated from northern China approximately 2,000 years ago (Kim and Choi, 2002). KNP has higher fat content in its muscle, and more chewy meat texture, but also has a lower growth rate and poorer productivity than other pig breeds (Choi et al., 2015). To uncover the genomic determinants of KNP-specific traits, it is necessary to compare the genomes of many different pig breeds as well as other related species. To this end, there have been genome-wide association studies aimed at identifying genomic elements, such as genes and single nucleotide polymorphisms (SNPs), that have relevance to economically important traits of KNPs, such as meat quality and productivity (Choi et al., 2015; Edea and Kim, 2014; Lee et al., 2012a; 2012b). These genomic bases for specific traits have gained their functions as a result of natural or artificial selection (selective breeding) during evolution. Therefore, studies designed to explore evolutionary perspectives can contribute to extending our understanding of KNP-specific genomic features and their functional consequences. This type of study was recently performed for Korean cattle (Lee et al., 2016) and the Korean bob-tailed native dog (Lee et al., 2017) to discover evolutionary changes of SNPs, genes, and interactions, which were associated with their characteristic traits. These studies also led to the development of a web interface for easily identifying the rewiring of protein-protein interactions that have occurred during evolution (Kwon et al., 2018).

In this study, we identified SNPs and genes relevant to KNP-specific traits, and discovered the mode and extent of their evolutionary changes using NGS data of six mammalian species (cattle, sheep, horse, dog, cat, and human) and six pig breeds (Korean wild boars [KWB], Duroc [DUR], Landrace [LAN], Yorkshire [YOR], Minipig [MINI], and KNP). We predicted the evolutionary trajectory of KNP SNPs within the used species as well as pig breeds. In addition, inter-species KNP-specific genes were identified based on the unique presence of SNPs in KNP, and their association with the olfactory function was discovered by enrichment tests using the Gene Ontology (GO) and pathway databases. Finally, the potential effects of intra-species KNP-specific nsSNPs within the different pig breeds on KNP-specific traits were examined by using pig quantitative trait loci regions. This study provides new insight into the evolution of the KNP and serves as a good example of the use of both population and comparative genomic approaches for a better understanding of domestic animals in terms of evolution and domestication.

MATERIALS AND METHODS

Species selection

The six mammalian species, Bos Taurus (cattle), Canis lupus familiaris (dog), Equus caballus (horse), Felis catus (cat), Ovis aries (sheep), and Homo sapiens (human), were used for inter-species evolutionary analyses. These species were selected based on several criteria, including relatedness to the pig, the existence of chromosome-level assembly, and the availability of the dbSNP data. Additionally, human was selected to be used as an outgroup species for evolutionary analyses. Six pig breeds, KWB, KNP, DUR, YOR, LAN, and MINI, were used for intra-species population analyses.

Sequence and annotation data preparation

The genome sequences of the seven mammalian species were obtained from the National Center for Biotechnology Information (NCBI) RefSeq database (O'Leary et al., 2016). Genome sequence identifiers are GCF_000003055.6 (cattle, bosTau8), GCF_000002285.3 (dog, canFam3), GCF_000002285.3 (horse, equCab2), GCF_000181335.2 (cat, felCat8), GCF_000298735.2 (sheep, oviAri3), GCF_000003025.6 (pig, susScr11), and GCF_000001405.33 (human, hg38). Next generation sequencing (NGS) data of six pig breeds was obtained from previous studies (Choi et al., 2015; Kim et al., 2015). The NGS data was generated from a total of 59 samples, including 10 of KNPs, 10 of KWBs, 6 of DURs, 14 of LAN, 15 of YOR, and 4 of MINI. Also, the coding DNA sequences and amino acid sequences of pig (assembly version: Sscrofa11.1), which are required for gene annotation and functional enrichment analysis, were downloaded from the Ensembl database (Zerbino et al., 2018). Gene annotation information which is formatted as the GTF format was downloaded from the Ensembl database. And, for the functional enrichment analyses, GO and orthologous gene information were obtained by the Ensembl Genes 95 version of the Ensembl biomart (Kinsella et al., 2011).

Single nucleotide polymorphism identification

SNPs were detected for each pig breed separately by integrating NGS data from multiple samples. The SNP calling and filtration were performed by the guidelines of the GATK Best Practice pipelines (Van der Auwera et al., 2013). SNPs were called using BWA-mem v0.7.15-r1140 with an additional -M option (Li, 2013) for NGS read alignment on the reference pig genome (assembly version: Sscrofa11.1), Samtools v1.3.1 (Li, 2011) and Picard v2.17.11 (https://broadinstitute.github.io/picard/) for preprocessing of alignments, and GATK v3.8 (DePristo et al., 2011) for SNP calling and filtering. To obtain high-quality SNPs, we filtered the called SNPs by several criteria: QD (Variant Confidence/Quality by Depth) < 2.0, MQ (RMS Mapping Quality) < 40.0, FS (Phred-scaled P value using Fisher's exact test to detect strand bias) > 60.0, HaplotypeScore (Consistency of the site with at most two segregating haplotypes) > 13.0, MQRankSum (Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities) < –2.5, and ReadPosRankSum (Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias) < –8.0. Then, the transition to transversion (Ti/Tv) ratio was calculated to ensure the quality of the called SNPs by vcftools v0.1.17 with --TsTv-summary option (Danecek et al., 2011). For seven mammalian species, SNP datasets were downloaded from the dbSNP database (build 148 for cat; 150 for cattle and pig; 151 for dog, horse, sheep and human) (Sherry et al., 2001).

Identification of orthologous genomic positions for SNP evolutionary analysis

Pairwise whole-genome sequence alignments of five species (cattle, sheep, horse, dog, and cat) against pig were built by using the LASTZ program (Harris, 2007) with E = 30, H = 2000, K = 3000, L = 2200, M = 50, O = 400, T = 1, and Y = 9400 options. The pairwise whole-genome sequence alignment between pig and human was downloaded from the UCSC genome browser (Haeussler et al., 2019). By using the whole-genome sequence alignments and the liftOver tool (Kent et al., 2002), the orthologous genomic positions of KNP SNPs in other mammalian species were identified.

Phylogenetic information preparation

The phylogenetic information among the seven mammalian species (pig, cattle, sheep, horse, dog, cat, and human) was downloaded from the TimeTree website (Kumar et al., 2017). The phylogenetic relationships among the six pig breeds were estimated using the SNPhylo program (Lee et al., 2014) based on their SNP data. In the case of the six pig breeds, we used only the topology, not branch lengths, of the estimated phylogenetic tree.

Inference of SNP changes during evolution

The evolutionary changes, such as appearance and disappearance of KNP SNPs during evolution, were inferred by using the maximum parsimony algorithm which minimizes the total number of changes during the evolutionary process (Farris, 1970; Fitch, 1971). First, the existence or absence of SNPs at the orthologous genomic positions in other species were identified. Then, the existence or absence of SNPs at each SNP position in the common ancestors was predicted by considering the existence of SNPs at orthologous genomic positions in all species using the maximum parsimony algorithm. Human was used as an outgroup species to resolve ambiguous cases. Next, we inferred the appearance or disappearance of them at each branch by comparing the SNP existence between all ancestor-descendant pairs. The SNP appearance and disappearance were defined as SNP absence in ancestor and existence in the descendant, and SNP existence in ancestor and absence in the descendant, respectively. In the case of evolutionary analysis within six pig breeds, a similar analysis was performed by using the other six mammalian species (except pig) as outgroup species. The maximum parsimony algorithm and inference of evolutionary changes were implemented by in-house scripts using the Perl programming language.

Definition of the evolutionary speed in inter- and intra-species evolutionary analysis

In the inter-species evolutionary analysis, for comparing the evolutionary speed of each branch, the rate of SNP changes was defined as the number of appeared and disappeared SNPs per billion year. For intra-species evolutionary analysis, the ratio of SNPs changes was used, which is the fraction of SNP changes against the number of SNPs in the immediate common ancestor nodes.

Identification of positively selected genes

To identify positively selected genes (PSGs) in KNP, we generated multiple protein sequence alignments for the inter-species KNP-specific genes and all corresponding orthologous genes (Ensembl Genes 95) in other mammalian species using MUSCLE v3.8.425 (Edgar, 2004). Out of 1,747 inter-species KNP-specific genes, a total of 252 genes with one or more corresponding orthologous genes in each of the other mammalian species were used in this analysis. Gaps in the multiple protein alignment results were removed and the alignments were transformed to DNA sequence alignments by using trimAl v1.4. with -phylip_paml -nogaps -splitbystopcodon options (Capella-Gutierrez et al., 2009). Also, the DNA sequence alignments were filtered out by the length cutoff (<150 bp). The PSGs were identified by using the codeml program in the PAML package (Yang, 1997) with branch-site models which were specified by the following model parameters: model = 2, NSsites = 2, omega = 1, and fix_omega = 1 for a null model and fix-omega = 0 for an alternative model. The alternative and null models were used for modeling the dN/dS ratios on branches, and the dN/dS ratio was fixed to 1 for all branches in the null model. The Chi-square test was used to check the statistical significance of each gene, and q-values (False Discovery Rate-adjusted P value) were calculated for multiple test correction. For some KNP genes, there were multiple orthologous genes in other species. However, only one of the orthologous genes per species can be used in the test for PSGs. Therefore, the codeml program was run for each of different combinations of orthologous genes (one from each species) with the q-value cutoff 0.05, and the KNP genes with at least one of the combinations passing the q-value cutoff were considered as the PSG.

SNP annotation, and functional enrichment and pathway analysis

The KNP-specific SNP sets were annotated, and their functional effects were predicted by the SnpEff program v4.3t (Cingolani et al., 2012). The pig reference genome, CDS sequences, and gene annotation information of pig (assembly version: Sscrofa11.1) were used for SNP annotation. Non-synonymous SNPs (nsSNPs) were identified by selecting the SNPs which were classified into the four sequence ontology terms: missense variant, start loss, stop gained and stop loss. Functional enrichment analyses of the annotated genes were performed by enrichment test with GO information. The enrichment tests were conducted by the hypergeometric test, and q-values were calculated for multiple test correction with 0.05 as a cutoff. The statistical tests and corrections were performed with R (https://www.R-project.org). Additionally, we performed the pathway analysis using the additional pathway databases, KEGG (Kanehisa et al., 2019) and REACTOME (Fabregat et al., 2018), by g:GOSt, which is the functional profiling module of g:Profiler (Reimand et al., 2016), with default parameters.

Analysis with quantitative trait loci

Pig quantitative trait loci (QTL) data was obtained from the Animal QTL database (Release 39; pig assembly version SS_11.1) (Hu et al., 2019). To estimate the degree of association between the intra-species KNP-specific nsSNPs and QTLs, the number of intra-species KNP-specific nsSNPs which were placed in the QTL regions was normalized by million base pair. The example of intra-species KNP-specific nsSNPs was visualized by using the UCSC genome browser track hub (Raney et al., 2014) with pig RefSeq genes (O'Leary et al., 2016).

Data access

All SNPs of the six pig breeds detected in this study have been submitted to the European Variation Archive database (https://www.ebi.ac.uk/eva/; project ID: PRJEB34412). The intra-species KNP-specific nsSNP data within the pig QTL regions is available on My Hubs page on the UCSC track hub website (https://genome.ucsc.edu/cgi-bin/hgHubConnect) using the hub URL (http://bioinfo.konkuk.ac.kr/KNP_nsSNPs/hub.txt).

RESULTS

Evolutionary analysis of KNP SNPs among related species

A total of 36,283,251 SNPs was called from six pig breeds: KNP, KWB, DUR, LAN, YOR, and MINI (Table 1; Materials and Methods section). Among them, more than 80% of SNPs in each pig breed were found in the dbSNP database (build 150). The transition to transversion (Ti/Tv) ratio was ranged from 2.31 to 2.36, which is similar to the numbers found in other studies of pigs (Choi et al., 2015; Li et al., 2017). For evolutionary analysis, additional SNPs of other six mammalian species (cattle, sheep, horse, dog, cat, and human) were downloaded from the dbSNP database (Supplementary Table S1; Materials and Methods section).

Table 1

Statistics of SNPs of six pig breeds

Pig breeds	No. of identified SNPs	No. of SNPs in dbSNPa	Ti/Tv ratio
KNP	16,545,631	13,930,964 (0.84)	2.32
KWB	23,528,922	19,587,524 (0.83)	2.36
DUR	11,518,093	9,677,013 (0.84)	2.32
LAN	19,717,820	16,346,309 (0.83)	2.29
YOR	18,422,733	15,161,229 (0.82)	2.27
MINI	10,293,949	8,600,902 (0.84)	2.31
Total	36,283,251	29,077,408 (0.80)

The fractions of identified SNPs found in dbSNP (build 150; total 58,431,079 pig SNPs) are shown in parentheses. For calling SNPs, the pig reference assembly version susScr11 was used.

To understand the evolutionary patterns of KNP SNPs, we traced their evolutionary history based on SNP information in other six mammalian species shown above, orthologous genomic positions of KNP SNPs in the used species, and the predicted SNPs in the common ancestors of the used species by the maximum parsimony algorithm (Materials and Methods section). In this analysis, only 3,779,064 KNP SNPs, which have their orthologous genomic positions in all other six mammalian species, were used. As shown in Fig. 1, less than 4% out total of 3,779,064 KNP SNPs were observed at orthologous genomic positions in other species, and the percentage was decreased as a function of divergence from KNP. A total of 7,215 KNP SNPs in A1, 249,822 KNP SNPs in A2, 249,175 KNP SNPs in A3, 3,487 KNP SNPs in A4, and 198 KNP SNPs in A5 were predicted as being present in each ancestor, which also showed a similar decreasing pattern of the number of KNP SNPs as divergence time increases. Those ancestral SNPs were further compared to discover the evolutionary changes of KNP SNPs by calculating the rate of SNP changes, which was defined as the number of appeared and disappeared SNPs per billion year. For branches from A1 to KNP, a total of 15.366 KNP SNPs per billion year appeared from A1 to A2, and this rate increased to 56.955 from A2 to KNP. Among KNP, cattle, and sheep, the number of ancestral SNPs was almost preserved from the ancestor A2 to A3, and more than 40% of SNPs disappeared during evolution from A3 to cattle and sheep. As expected, a very small number of SNPs was predicted to be present in A4 compared with A2. A4 is the common ancestor of horse, dog, and cat, which are more divergent species than cattle and sheep from KNP. Interestingly, the number of SNPs was increased from A4 to horse, and from A5 to dog and cat, which was opposite to the pattern observed from A2 to KNP, cattle, and sheep. This was because compared with the number of orthologous genomic positions with SNPs in each species (horse, dog, and cat), only a small number of genomic positions have common SNPs in all three species, which led to the prediction of a small number of ancestral SNPs in A4 and A5.

Fig. 1

Evolution of KNP SNPs.

The number next to KNP represents the number of KNP SNPs used in this analysis, and the numbers next to the other descendant species show the number of SNPs found at orthologous genomic positions of KNP SNPs in those species. The numbers next to the names of ancestors (A1 to A5) are the number of KNP SNPs predicted to be present in those ancestors. The numbers on branches represent the rate of SNP changes which is the number of appearing and disappearing events per billion year on each branch. The rate of appearance and disappearance events are separately shown in parentheses with + and – sign respectively. The zero rate in parentheses with an asterisk (*) symbol means a very small rate that cannot be shown with a current precision of numbers.

Functional enrichment and pathway analysis of KNP genes with SNPs only in KNP

The SNPs in six mammalian species (KNP, cattle, sheep, horse, dog, and cat) were used to identify KNP genes with SNPs not found at orthologous genome positions in the other five species (hereafter “inter-species KNP-specific genes”; 1,747 genes found). The GO enrichment test was then performed with q-value cutoff 0.05 (Materials and Methods section). The inter-species KNP-specific genes were enriched in three major categories of function; olfaction, immune system, and translation (Supplementary Table S2).

The aforementioned 1,747 inter-species KNP-specific genes were further reduced to genes with only nsSNPs (total 1,016 genes found). The GO enrichment test was also performed as above, and a total of 34 functions in the three GO categories were significantly enriched for those genes (Table 2,Supplementary Table S3). In terms of the biological process category, olfaction-associated GO terms, such as G protein-coupled receptor signaling pathway (GO:0007186), signal transduction (GO:0007165), detection of chemical stimulus involved in sensory perception of smell (GO:0050911), response to stimulus (GO:0050896), and sensory perception of smell (GO:0007608), were significantly enriched. Pathway analysis using g:Profiler (Reimand et al., 2016) also confirmed that those genes are significantly related to olfactory function (Materials and Methods section; Supplementary Table S4). For example, olfaction-related pathways, such as olfactory transduction (KEGG:04740) in the KEGG database (Kanehisa et al., 2019), and olfactory signaling pathway (REAC:R-SSC-381753), G alpha (s) signaling events (REAC:R-SSC-418555), GPCR downstream signaling (REAC:R-SSC-388396), signaling by GPCR (REAC:R-SSC-372790), signal transduction (REAC:R-SSC-162582), regulation of IFNA signaling (REAC:R-SSC-912694), and interferon alpha/beta signaling (REAC:R-SSC-909733) in the Reactome database (Fabregat et al., 2018), were found to be significantly enriched. The 1,016 genes were also significantly enriched for the GO terms related to the immune system, such as natural killer cell activation involved in immune response (GO:0002323), T cell activation involved in immune response (GO:0002286), positive regulation of peptidyl-serine phosphorylation of STAT protein (GO:0033141), B cell proliferation (GO:0042100), defense response (GO:0006952), response to exogenous dsRNA (GO:0043330), and humoral immune response (GO:0006959).

Table 2

Significantly enriched functions (q-value cutoff 0.05) of 1,016 inter-species KNP-unique genes containing non-synonymous SNPs only in KNP not in the five other mammalian species

GO ID	GO term	No. of associated genes	q-value
GO:0007186	G protein-coupled receptor signaling pathway	523	< 1.0E-05
GO:0007165	Signal transduction	515	< 1.0E-05
GO:0050911	Detection of chemical stimulus involved in sensory perception of smell	514	< 1.0E-05
GO:0050896	Response to stimulus	490	< 1.0E-05
GO:0007608	Sensory perception of smell	489	< 1.0E-05
GO:0006412	Translation	105	< 1.0E-05
GO:0002323	Natural killer cell activation involved in immune response	9	4.0E-04
GO:0002286	T cell activation involved in immune response	9	6.9E-04
GO:0033141	Positive regulation of peptidyl-serine phosphorylation of STAT protein	9	6.9E-04
GO:0022900	Electron transport chain	15	8.5E-04
GO:0006952	Defense response	16	2.7E-03
GO:0042100	B cell proliferation	9	7.9E-03
GO:0006826	Iron ion transport	9	1.5E-02
GO:0043330	Response to exogenous dsRNA	9	1.5E-02
GO:0006959	Humoral immune response	9	1.7E-02
GO:0098664	G protein-coupled serotonin receptor signaling pathway	7	2.7E-02

The functions in the biological process category are shown (others are in Supplementary Table S3).

To find KNP genes strongly affected by selective pressure during evolution, PSGs in KNP were identified from 252 out of 1,747 inter-species KNP-specific genes, which have orthologous genes in all of the other six mammalian species (cattle, sheep, horse, dog, cat, and human; Materials and Methods section). A total of 48 PSGs were discovered (Supplementary Table S5), and the GO enrichment test showed that most of the PSGs were enriched in olfaction-related GO terms, such as sensory perception of smell (GO:0007608), response to stimulus (GO:0050896), and olfactory receptor activity (GO:0004984) (Supplementary Table S6).

Evolutionary analysis of KNP SNPs within pig breeds

To study the evolution of KNP SNPs during the domestication process of six pig breeds (KNP, KWB, DUR, LAN, YOR, and MINI), the presence and absence of KNP SNPs in ancestors were predicted using the maximum parsimony algorithm based on the phylogenetic tree obtained from the SNP data (Materials and Methods section; Fig. 2). In the phylogenetic tree, a subtree was created with KNP and the other four pig breeds (DUR, MINI, YOR, and LAN), and KWB was placed outside of the subtree. From 44% (MINI) to 72% (LAN) of KNP SNPs were shared with the other pig breeds. The existence information of KNP SNPs in other pig breeds was used to predict the presence of KNP SNPs in ancestors, and 11,223,549 KNP SNPs in P1, 14,163,099 KNP SNPs in P2, 11,466,536 KNP SNPs in P3, 9,497,632 KNP SNPs in P4, and 11,387,129 KNP SNPs in P5 were predicted as present. Note that the number of KNP SNPs in this analysis was smaller than the number in Table 1 because some of KNP SNPs which failed to produce unambiguous ancestral SNPs were filtered out.

Fig. 2

Changes of SNPs during pig domestication.

The numbers next to the names of ancestors (P1 to P5) and descendant pig breeds represent the numbers of existing SNPs. The numbers on branches represent the ratio of SNP changes which is the number of appeared and disappeared SNPs against the number of present SNPs of an immediate ancestor (aSNPs). The ratio of SNP appearance and disappearance are separately shown in parentheses with + and – sign, respectively. The numbers which are higher in disappeared SNPs than appeared SNPs are represented in bold, and the red-dashed rectangle highlights the group that has undergone strong negative selection pressure. The tree topology was predicted by SNPhylo. The zero rate in parentheses with an asterisk (*) symbol means a very small rate that cannot be shown with a current precision of numbers.

For profiling and comparing the evolutionary changes of KNP SNPs in the used pig breeds, we calculated the ratio of SNP changes, which is defined as the number of SNP changes against the number of SNPs in an immediate ancestor on each branch. Note that the rate per billion year, which was used in the inter-species evolution study, was not used here because the exact time of domestication of each breed is hard to obtain. For branches from P1 to KNP, the gradual increase of SNPs was observed (ratio +0.262 from P1 to P2, +0.152 from P2 to KNP). However, the number of ancestral SNPs was decreased from P2 to four pig breeds (DUR, MINI, YOR, and LAN), and the ratio of SNP disappearance from P3 to P4 (ratio –0.172) was almost 25 folds greater than the ratio from P3 to P5 (ratio –0.007). The biggest and smallest SNP changes were observed on branches from P4 to MINI (ratio 0.295) and from P1 to KWB (ratio ~0), respectively.

Functional enrichment and QTL analysis of KNP genes with SNPs absent in other pig breeds

From the comparison of SNP positions in six pig breeds, a total of 510,474 SNPs was identified as being present only in KNP. Among them, 4,304 were nsSNPs (hereafter “intra-species KNP-specific nsSNPs”). To examine the effect of these intra-species KNP-specific nsSNPs, we performed a functional enrichment analysis of 2,608 genes containing at least one of the 4,304 intra-species KNP-specific nsSNPs. Among the 2,608 genes, a total of 503 genes were associated with the 11 significantly enriched GO terms, and most of them were enriched in functions associated with olfaction (Table 3).

Table 3

Significantly enriched functions of 2,608 KNP genes with unique non-synonymous SNPs

GO ID	GO term	GO domain	No. of associated genes	q-value
GO:0007165	Signal transduction	BP	295	< 1.0E-05
GO:0007186	G protein-coupled receptor signaling pathway	BP	307	< 1.0E-05
GO:0007608	Sensory perception of smell	BP	277	< 1.0E-05
GO:0050896	Response to stimulus	BP	278	< 1.0E-05
GO:0050911	Detection of chemical stimulus involved in sensory perception of smell	BP	294	< 1.0E-05
GO:0005886	Plasma membrane	CC	310	< 1.0E-05
GO:0016020	Membrane	CC	454	< 1.0E-05
GO:0016021	Integral component of membrane	CC	451	< 1.0E-05
GO:0004930	G protein-coupled receptor activity	MF	289	< 1.0E-05
GO:0004984	Olfactory receptor activity	MF	294	< 1.0E-05
GO:0003735	Structural constituent of ribosome	MF	38	4.4E-02

BP, biological process; CC, cellular component; MF, molecular function.

To further investigate the association with interesting traits, the intra-species KNP-specific nsSNPs were compared with pig QTL based on the Pig QTLdb (Hu et al., 2019). We obtained a total of 4,241 intra-species KNP-specific nsSNPs positions within 4,864 QTLs (Supplementary Table S7). Based on the normalized number of intra-species KNP-specific nsSNPs per million base pairs in each QTL, the top-ranked QTLs, including HDL/LDL ratio QTL, feed conversion ratio QTL, intramuscular fat content QTL, marbling QTL, and muscle moisture percentage QTL, were examined. We found that the intra-species KNP-specific nsSNPs are highly related to meat quality, especially, marbling and muscle moisture level. The intramuscular fat content QTL region, which contains 21 intra-species KNP-specific nsSNPs and the TAS1R2 gene associated with dietary adaptation (da Silva et al., 2014), is visualized in Fig. 3 by using the UCSC genome browser track hub (Raney et al., 2014). In this region, there are 18 pig genes and 21 intra-species KNP-specific nsSNPs (Fig. 3A). Also, the example of the TAS1R2 gene shows an amino acid change caused by the intra-species KNP-specific nsSNP (Fig. 3B).

Fig. 3

Examples of intra-species KNP-specific nsSNPs.

(A) Pig RefSeq genes were represented by blue blocks with their gene symbols and coordinates on the top panel. The intra-species KNP-specific nsSNP information is represented by a vertical line containing a pair of two alleles (top: reference allele, bottom: alternative allele) under the gene annotations, and their relative fractions are shown by vertical lines with different heights for each allele. This track was displayed by the UCSC genome browser track hub. (B) In the TAS1R2 gene example, the top panel shows gene structure with the direction of transcription (blue arrow), and bottom panel represents the comparison of amino acids in the third exon among six different pig breeds. Two different amino acids caused by the intra-species KNP-specific nsSNP are shown together with a slash delimiter.

Commonality of nsSNPs and genes found in between inter- and intra-species analysis

The KNP-specific nsSNPs and genes obtained from the inter- and intra-species analysis were further examined to find their commonality. Among 7,992 inter-species and 4,304 intra-species KNP-specific nsSNPs, 662 common nsSNPs were found (Supplementary Fig. S1A). These common nsSNPs contributed to 331 common genes between 1,747 inter-species and 2,608 intra-species KNP-specific genes (Supplementary Fig. S1B). The 331 common genes include the gene (ENSSSCG00000000615) annotated in the backfat-thickness related QTL region, which is one of the KNP-specific traits (Paszek et al., 2001), and other genes (ENSSSCG00000001827, ENSSSCG00000040964) related to the KNP-specific phenotypes such as carcass length (de Koning et al., 2003; Falker-Gieske et al., 2019).

DISCUSSION

In this study, the evolutionary and functional analysis of KNP with the six mammalian species (cattle, sheep, horse, cat, dog, and human) and five different pig breeds (KWB, DUR, LAN, YOR, and MINI) were performed by using the SNP data. Through the evolutionary analysis, we discovered how KNP SNPs have generated during evolution within the used species as well as pig breeds. The inter-species KNP-specific genes with the uniquely present SNPs only in KNP and their functions associated with interesting traits were discovered using the GO and pathway databases. In addition, the relationships between intra-species KNP-specific nsSNPs and pig quantitative trait loci regions were examined.

In the evolutionary analysis, the changes of KNP SNPs during evolution were investigated by using the information of shared KNP SNPs in other mammalian species (Fig. 1) as well as other pig breeds (Fig. 2), and they were quantified in terms of divergence time (in the case of inter-species analysis) or the ratio of SNP changes (in the case of inter-pig breed analysis). To reduce confounding effect by the different number of known SNPs in the used species (or pig breeds), we only focused on KNP SNPs and their existence at the orthologous genomic positions in the other species (or pig breeds). If reliable evolutionary trajectories of SNPs of other species (or pig breeds) are obtained, the evolutionary history of KNP SNPs can be more characterized, and distinguishing evolutionary patterns compared with other species (or pig breeds) can be discovered.

A total of 1,016 out of 1,747 inter-species KNP-specific genes contained nsSNPs not found at orthologous genome positions in the other mammalian species (cattle, sheep, horse, dog, and cat). These genes were found in various data resources, such as the GO, KEGG and REACTOME database, as being significantly related to functions associated with olfaction, immune system and translation. These functions have been known as the characteristics of pigs associated with adaptation and domestication (Groenen et al., 2012; Vincent et al., 2012; Zhu et al., 2017). Pigs have the largest olfactory receptor gene families, including more than 1,300 functional olfactory receptor genes located on 16 pig chromosomes (Nguyen et al., 2012). Among the 1,016 inter-species KNP-specific genes, 177 genes were olfactory receptor genes (Supplementary Table S8) related with receptor signaling pathway and the detection of chemical stimulus involved in sensory perception of smell. In addition, immunity-related genes in pigs have highly expanded in the evolutionary process. For example, pigs have more than 50 type I interferon genes, which are significantly more than in human and mouse (Sang et al., 2014). Immune response is significantly associated with the feed efficiency and product quality of pigs (Horodyska et al., 2018). In our analysis, a total of 12 immunity related genes were observed in the above inter-species KNP-specific genes (Supplementary Table S9).

Out of the 1,747 inter-species KNP-specific genes, a total of 48 genes were predicted as PSGs. Among them, six genes (ENSSSCG00000036715, ENSSSCG00000038050, ENSSSCG00000037486, ENSSSCG00000033023, ENSSSCG00000038781, ENSSSCG00000040604) were classified as the G protein-coupled receptor and annotated with odorant binding GO term (GO:0005549) which was enriched in 48 PSGs. It suggests that KNPs have been characterized and diversified the functions associated with olfaction in the evolution process, which is in line with the characteristics of pigs found in previous evolutionary and comparative studies (Groenen et al., 2012; Zhu et al., 2017). Additional GO terms, such as single-stranded RNA binding (GO:0003727) in the molecular function category and transposition, RNA-mediated (GO:0032197) in the biological process category, were also enriched in the 48 PSGs. Seven genes (ENSSSCG00000032719, ENSSSCG0000003778, ENSSSCG000000364866, ENSSSCG00000036856, ENSSSCG00000040198, ENSSSCG00000040767, and ENSSSCG00000033328) annotated with those GO terms belong to the LINE-1 retrotransposable element ORF1 protein subfamily that is highly related to LINE-1 which is known as main repetitive element groups in pig (Groenen et al., 2012).

We identified 4,241 intra-species KNP-specific nsSNPs positioned within 4,864 pig QTL regions. The associations between the intra-species KNP-specific nsSNPs and QTL regions were evaluated by calculating the normalized number of intra-species KNP-specific nsSNPs per million base pairs. Interestingly, the top-ranked QTL regions were the intramuscular fat content QTL and muscle moisture percentage QTL that are known to having positive correlation with the meat quality-related features such as high-hat content in muscle, meat tenderness, juiciness, and taste (Grindflek et al., 2001). The intramuscular fat content QTL is associated with the fat content and fatty acid composition in muscle (Hocquette et al., 2010; Sanchez et al., 2007), and the muscle moisture percentage QTL is also related to the meat quality in pigs (Kärst et al., 2011; Wang et al., 2016; Watanabe et al., 2018). Those features may be related to the formation of the well-known features of KNP meat, which are high-fat content in muscle and chewy texture (Cho et al., 2015). In addition, QTLs, such as the feed conversion ratio QTL and the total number born alive QTL, were also highly ranked and seemed to be associated with the poor productivity of KNP. Of interest was the example of the TAS1R2 gene, which contains an amino acid change in KNP. The TAS1R gene family is known to be associated with taste and nutrient sensing as well as dietary adaptations (da Silva et al., 2014) (Fig. 3B), and therefore might also be related to the meat quality of KNP.

Lastly, the commonality of nsSNPs and genes between inter- and intra-species analysis was examined, and 662 nsSNPs within 331 genes were found as common. These nsSNPs and genes were observed in the KNP-specific traits-related QTL regions, such as the backfat-thickness related and carcass length QTL.

In summary, SNPs of KNP, other pig breeds, and related species were used to discover the mode and extent of the evolutionary changes of the KNP SNPs, inter- and intra-species KNP-specific nsSNPs, inter-species KNP-specific genes, their functions, and QTLs harboring those intra-species KNP-specific nsSNPs. As sequencing technologies become faster and cheaper, population-scale genomic studies have been actively performed for many species and breeds. Furthermore, comprehensive comparative and evolutionary analysis can help enhance a deeper understanding of species diversity and their history. For domestic animals in particular, it is very important to understand their genomic history and the characteristics which lead to phenotypic consequences, in order to preserve and improve their valuable traits. This study suggests a good example of both population and comparative genomic approaches for a better understanding of domestic animals in terms of evolution and domestication.

Supplemental Materials

Note: Supplementary information is available on the Molecules and Cells website (