Phenotypic Data Analysis
We performed a GWAS with different quantitative traits such as height, length, LHR, and body-weight in Jindo dog. We investigated 757 Jindo dogs with an average raw height of 48.51 cm, ranging from 33.5 to 59 cm and an average length of 51.42 cm ranging from 37 cm to 62 cm. The average LHR is 105.82, ranging from 86.4 to 130.7. The average raw weight in this individual dataset is 16.67 kg, ranging from 6.6 kg to 33.3 kg.
The size variance was noticed high because of different ages of dogs were included in our current study for example a medium dog (younger than 16 months), adult dog (up to 8 years), and senior dogs (8 years or older). The standard error of the mean value was 0.13, 0.13, 0.18, and 0.16 for height, length, LHR, and body weight, respectively. The distribution of dogs was done according to bodyweight quartiles such as Q1 (the lowest 25 % of registered weight) − 8.4, Q2 and Q3–25 % below (12.4) and above the median (20.8), Q4 – dogs with the highest 25% of registered weight values (24.6). A total of 9 (Q1), 46 (Q2), 603 (Q3), and 99 (Q4) dogs were in each body-weight quartile. The ratios of genotypic and phenotypic variance explained by all the SNPs were found to be 0.33, 0.08, 0.08, and 0.28 for height, length, LHR, and body weight, respectively. Figure 1 shows an example of phenotype measurement of male and female Jindo dogs. A descriptive statistical summary of the phenotypes included the minimum, maximum, and mean, as well as the standard deviations and standard error are given in Table 1.
Table 1
Descriptive summary of phenotype data.
Traits | | Phenotypic Data | | Variance |
| 492 Female, 265 Male | | | | |
| Mean | Median | Max | Min | SD | SE | h2 | Vp | V(G) |
Height (cm) | 48.61 | 48.5 | 59 | 33.5 | 3.69 | 0.13 | 0.33 | 13.26 | 4.31 |
Length (cm) | 51.42 | 51.3 | 62 | 37 | 3.62 | 0.13 | 0.08 | 106.68 | 8.81 |
LHR | 105.82 | 105.6 | 130.7 | 86.4 | 5.23 | 0.18 | 0.08 | 361.40 | 32.14 |
Weight (Kg) | 16.67 | 16.4 | 33.3 | 6.6 | 3.41 | 0.16 | 0.28 | 11.34 | 3.19 |
SD − standard deviation, SE − standard error, Vp − phenotypic variance, V(G) − genotypic variance, h2 − heritability, LHR − Length to height ratio. |
Genome-wide Association Study
Before carrying out the GWAS analysis, we analyzed the patterns of population structure through principal component analysis (PCA) and the population found to form a good clustering pattern (Fig. 2). We performed GWAS with a total of 757 dogs and 118,879 SNPs that passed after QC to discover significant loci associated with traits. The total genotyping rate was 0.988285. The Manhattan plot represents our mixed linear model-based GWAS result (Figs. 3, 4, 5 and 6). The association result revealed that a total of 40 SNPs in the genome were associated with all the four traits of height, length, LHR and body-weight at a genome-wide significant level (p < 5E-05). Among these SNPs, 10, 6, 13, and 11 SNPs were significant at the genome-wide level (p < 5E-05) for the trait of height, length, LHR and body-weight, respectively. A strong signal was identified only for height and body-weight trait, where one SNP (TIGRP2P201033_rs9187576) on Chr15 and 4 SNPs on Chr7 and Chr15 passed the Bonferroni corrected significance threshold (p < 4E-07). Most of these SNPs were located in the intronic regions and intergenic regions. No SNPs crossed Bonferroni corrected significance threshold (p < 4E-07) for the trait of length and LHR. The topmost significant (p < 5E-05) five associated SNPs and their annotated genes are presented in Table 2. For the height trait, top five SNPs (TIGRP2P201033_rs9187576, BICF2G630358507, BICF2G630358640, BICF2S23724992, and BICF2S2377250) were annotated to four genes (HHIP, LCORL, NCAPG, and P4HA1), respectively. For the length trait, the SNPs (BICF2P772349, BICF2P580549, BICF2S2366810, BICF2P517149, and BICF2G630778876) were annotated to four diferent genes such as FGFR3, PLCH1, IGF1, and SCAND3, respectively. For LHR, the SNPs (TIGRP2P119932_rs8658799, BICF2G630268912, BICF2G630490347, BICF2P491431, and BICF2G630507820) were successfully annotated to seven genes (BEGAIN-DLK1, TDRD1, EFEMP1, ZFHX3-PSMD7, and SLC23A2) genes, respectively. For body-weight, the SNPs (BICF2P279062, BICF2P979272, BICF2S23655947, BICF2S2336786, and BICF2P110929) were successfully annotated to seven genes (PTPN2, RASAL2, PARPBP-IGF1, IGF1, ENSCAFG00000001744-SNORD19) genes, respectively. The clear deviation between observed and expected p-values on the Q-Q plot for all the traits signifies good correlation results in the present study (Figs. 3, 4, 5 and 6).
Table 2
Top 5 SNPs associated with height, length, LHR, and body weight in Jindo dog.
Trait | SNP ID | Chr | Position | p-value | Freq | Gene | Type |
Height | TIGRP2P201033_rs9187576 | 15 | 43493060 | 1.49E-07 | 0.11 | HHIP | intron_variant |
| BICF2G630358507 | 3 | 91180492 | 6.91E-07 | 0.47 | LCORL | intron_variant |
| BICF2G630358640 | 3 | 91314456 | 2.1E-06 | 0.33 | NCAPG | intron_variant |
| BICF2S23724992 | 4 | 23625990 | 4E-06 | 0.10 | P4HA1 | intron_variant |
| BICF2S2377250 | 3 | 93890772 | 6.65E-06 | 0.11 | - | - |
Length | BICF2P772349 | 3 | 62323843 | 4.42E-06 | 0.29 | FGFR3 | Upstream_gene_variant |
| BICF2P580549 | 23 | 49471629 | 5.42E-06 | 0.27 | PLCH1 | Intron_variant |
| BICF2S2366810 | 31 | 40675954 | 3.71E-05 | 0.25 | - | - |
| BICF2P517149 | 15 | 41228320 | 3.98E-05 | 0.24 | IGF1 | Intron_variant |
| BICF2G630778876 | 35 | 25711203 | 4.18E-05 | 0.17 | SCAND3 | Intron_variant |
LHR | TIGRP2P119932_rs8658799 | 8 | 68874072 | 7.88E-07 | 0.36 | BEGAIN-DLK1 | Intron_variant |
| BICF2G630268912 | 28 | 25034623 | 8.78E-07 | 0.28 | TDRD1 | Intron_variant |
| BICF2G630490347 | 10 | 56634077 | 8.98E-07 | 0.39 | EFEMP1 | Intron_variant |
| BICF2P491431 | 5 | 79708771 | 1.20E-06 | 0.11 | ZFHX3-PSMD7 | Intergenic_region |
| BICF2G630507820 | 24 | 16702908 | 4.09E-06 | 0.07 | SLC23A2 | Intron_variant |
Weight | BICF2P279062 | 7 | 78202038 | 8.91E-08 | 0.03 | PTPN2 | intron_variant |
| BICF2P979272 | 7 | 21586876 | 1.54E-07 | 0.37 | RASAL2 | intron_variant |
| BICF2S23655947 | 15 | 41134849 | 2.12E-07 | 0.05 | PARPBP-IGF1 | intergenic_region |
| BICF2S2336786 | 15 | 41239301 | 2.86E-07 | 0.07 | IGF1 | intron_variant |
| BICF2P110929 | 14 | 10389714 | 1.06E-06 | 0.03 | ENSCAFG00000001744-SNORD19 | intergenic_region |
LHR − Length to height ration |
Genotype-phenotype Correlation
We acquired height, length, LHR, and body-weight measurements on genotyped Jindo dogs.
We then compared the traits (height, length, LHR, and body-weight) distributions between the different genotype classes at the topmost variants, TIGRP2P201033_rs9187576, BICF2P772349, TIGRP2P119932_rs8658799, and BICF2P279062, respectively. The box plots show the distribution of phenotypes among the different genotypes (Fig. 7). The solid lines in the box plot are the medians of phenotype’s per-genotype group (AA, AG, GG, CC, CG, and CG) of all four variants were described in Table 3. In addition, we calculated the mean height, length, LHR, and body weight per genotype (Table 3). Although the differences between all three genotype classes were not significant, it showed that the A-allele for SNP TIGRP2P201033_rs9187576 was correlated with a mean reduction of the height trait in Jindo. For SNP BICF2P772349 and TIGRP2P119932_rs8658799, the A-allele showed to be correlated with increased length and LHR in Jindo dog, respectively. The trait-increasing allele A for SNP TIGRP2P119932_rs8658799 and G for SNP TIGRP2P201033_rs9187576 was found to be a major allele for both the SNPs. This A allele was detected as a minor allele in the case of BICF2P772349 SNP of length trait. However, the AA genotype for SNP BICF2P772349 was found to be present only in a small number of total animals.
Table 3
Effect of top SNP genotypes on the height, length, LHR, and weight trait in Jindo dog.
Traits | SNP name | Genotype | Counts | Mean | Median | SD |
Height (cm) | TIGRP2P201033_rs9187576 | AA | 16 | 48.1 | 48.1 | 3.96 |
| | AG | 144 | 48.37 | 48.50 | 3.82 |
| | GG | 597 | 48.70 | 49 | 3.60 |
Length (cm) | BICF2P772349 | AA | 2 | 53.26 | 53.00 | - |
| | AG | 126 | 51.08 | 52.00 | 3.49 |
| | GG | 629 | 51.44 | 52.00 | 3.64 |
LHR | TIGRP2P119932_rs8658799 | AA | 306 | 105.93 | 104.90 | 4.98 |
| | AG | 349 | 105.69 | 104.86 | 4.32 |
| | GG | 102 | 105.67 | 104.85 | 5.86 |
Weight (kg) | BICF2P279062 | CC | 95 | 16.14 | 16.05 | 3.36 |
| | CG | 335 | 16.96 | 17 | 3.72 |
| | GG | 327 | 16.53 | 16 | 3.08 |
LHR − Length to height ration |
Gene-set Enrichment And Pathway Analysis
Functional enrichment analysis was carried out to identify classes of genes that are over-represented in a large group of genes and may have a connection with the studied phenotypes. Many genes are in GO term and KEGG pathway categories. From the GWAS result, 1222, 1132, 1102 and 1039 SNPs (nominal p < 0.01) were used to annotate 842, 817, 752, and 718 genes correlated with height, length, LHR and weight trait (Supplementary Table S1). Subsequently, the genes were uploaded to the bioinformatics tool, DAVID for finding the enriched GO terms and KEGG pathways. The analysis revealed that a total of 72 GO terms and 37 KEGG pathways for height, 69 GO terms and 19 KEGG pathways for length, 53 GO terms and 18 KEGG pathways for LHR, 50 GO terms and 25 KEGG pathways for weight traits were significantly (p < 0.05) enriched. Out of the total enriched GO terms and KEGG pathways, the top 5 significantly enriched GO terms and pathways are presented in Table 4. For height, the top 5 GO terms were GO:1902711 − GABA-A receptor complex, GO:0004890 − GABA-A receptor activity, GO:0006816 − Calcium ion transport, GO:0030336 − Negative regulation of cell migration, and GO:0045665 − Negative regulation of neuron differentiation, and KEGG pathways were cfa04024 − cAMP signaling pathway, cfa04723 − Retrograde endocannabinoid signaling, cfa04261 − Adrenergic signaling in cardiomyocytes, cfa05033 − Nicotine addiction, and cfa04540 − Gap junction. For length, the top 5 GO terms were GO:0030425 − Dendrite, GO:0005096 − Gtpase activator activity, GO:0043025 − Neuronal cell body, GO:0032332 − Positive regulation of chondrocyte differentiation, and GO:0042593 − Glucose homeostasis, and KEGG pathways were cfa04922 − Glucagon signaling pathway, cfa04915 − Estrogen signaling pathway, cfa04024 − cAMP signaling pathway, cfa04918 − Thyroid hormone synthesis, and cfa04080 − Neuroactive ligand-receptor interaction. The top 5 GO terms (GO:0048706 − Embryonic skeletal system development, GO:0010595 − Positive regulation of endothelial cell migration, GO:0050840 − Extracellular matrix binding, GO:0048484 − Enteric nervous system development, and GO:1901166 − Neural crest cell migration involved in autonomic nervous system development) and KEGG pathways (cfa04911 − Insulin secretion, cfa04080 − Neuroactive ligand-receptor interaction, cfa04931 − Insulin resistance, cfa04915 − Estrogen signaling pathway, and cfa05202 − Transcriptional misregulation in cancer) were significantly enriched for LHR. For weight, the top 5 GO terms were GO:0030425 − Dendrite, GO:0005096 − Gtpase activator activity, GO:0017137 − Rab gtpase binding, GO:0045666 − Positive regulation of neuron differentiation, and GO:0050852 − T cell receptor signaling pathway, and KEGG pathways were cfa04921 − Oxytocin signaling pathway, cfa05100 − Bacterial invasion of epithelial cells, cfa04713 − Circadian entrainment, cfa04918 − Thyroid hormone synthesis, and cfa04520 − Adherens junction.
Table 4
Top 5 Gene Ontology and KEGG pathways significantly enriched using genes associated with height, length, LHR and body weight.
Trait | Category | Term_ID | Term | Count | % | P-value | Genes |
Height | GOTERM_CC_DIRECT | GO:1902711 | GABA-A receptor complex | 6 | 0.004585 | 6.E-04 | GABRA2, GABRB2, GABRA6, GABRA5, GABRA4, GABRG1 |
| GOTERM_MF_DIRECT | GO:0004890 | GABA-A receptor activity | 6 | 0.004585 | 7.E-04 | GABRA2, GABRB2, GABRA6, GABRA5, GABRA4, GABRG1 |
| GOTERM_BP_DIRECT | GO:0006816 | Calcium ion transport | 8 | 0.006114 | 8.E-04 | PPP3CA, RYR2, CAMK2D, PLN, PRKCB, CCL5, CAV1, NMUR2 |
| GOTERM_BP_DIRECT | GO:0030336 | Negative regulation of cell migration | 11 | 0.008407 | 9.E-04 | RNF20, KANK1, DACH1, STK24, DAG1, ADARB1, LDLRAD4, TBX5, SULF1, RECK, ROBO1 |
| GOTERM_BP_DIRECT | GO:0045665 | Negative regulation of neuron differentiation | 9 | 0.006878 | 9.E-04 | DDX6, ZHX2, MEIS1, ID4, OLIG2, CNTN4, PBX1, GLI3, CDK5RAP2 |
| KEGG_PATHWAY | cfa04024 | Camp signaling pathway | 23 | 0.017577 | 6.E-06 | GRIA1, RYR2, GABBR2, CAMK2D, CREBBP, SUCNR1, PPP1R12A, BDNF, HHIP, ADCY2, ADRB1, ATP2B1, GLI3, GNAI1, MAPK10, TIAM1, GRIN3A, PLN, FXYD2, CREB3L2, DRD1, CREB5, GRIA4 |
| KEGG_PATHWAY | cfa04723 | Retrograde endocannabinoid signaling | 16 | 0.012228 | 1.E-05 | GABRA2, GRIA1, GABRB2, GABRA6, GABRA5, PRKCB, GABRA4, ADCY2, GNAI1, GRM1, GABRG1, MAPK10, RIMS1, GRM5, PLCB1, GRIA4 |
| KEGG_PATHWAY | cfa04261 | Adrenergic signaling in cardiomyocytes | 17 | 0.012992 | 8.E-05 | RYR2, CAMK2D, PPP2R5A, ADCY2, ADRB1, ATP2B1, ADRA1B, GNAI1, PLN, PPP2R2B, PPP2R5E, KCNQ1, FXYD2, CREB3L2, PLCB1, CACNG3, CREB5 |
| KEGG_PATHWAY | cfa05033 | Nicotine addiction | 9 | 0.006878 | 1.E-04 | GABRA2, GRIA1, GABRB2, GRIN3A, GABRA6, GABRA5, GABRA4, GRIA4, GABRG1 |
| KEGG_PATHWAY | cfa04540 | Gap junction | 13 | 0.009935 | 2.E-04 | GUCY1A3, GUCY1B3, GUCY1A2, PRKCB, LPAR1, ADCY2, ADRB1, GNAI1, GRM1, TUBB6, GRM5, DRD1, PLCB1 |
Table 4
Trait | Category | Term_ID | Term | Count | % | P-value | Genes |
Length | GOTERM_CC_DIRECT | GO:0030425 | Dendrite | 21 | 2.781457 | 8.E-07 | GABRA2, EPHA4, SYT4, GABRA5, KCNB1, DENND1A, MAGI2, KLHL1, PDYN, GRM1, BRINP2, HTR7, RELN, PLCB4, BRINP1, GRM7, DNER, NCS1, MAPT, BRINP3, TP63 |
| GOTERM_MF_DIRECT | GO:0005096 | Gtpase activator activity | 21 | 2.781457 | 8.E-05 | RGS18, RABGAP1, DAB2IP, ARHGAP18, RASAL2, ARFGAP3, ARHGAP15, ARHGAP24, TIAM2, TBC1D1, TBCK, TBC1D2, SIPA1L1, TBC1D21, DLC1, TBC1D13, RAPGEF2, PLCB1, DOCK1, GARNL3, RGS7 |
| GOTERM_CC_DIRECT | GO:0043025 | Neuronal cell body | 15 | 1.986755 | 7.E-04 | GABRA2, DENND1A, SEZ6, KLHL1, GRIK2, PDYN, BRINP2, UCHL1, GRIN3A, BRINP1, DNER, RAPGEF2, APOB, ASIC2, BRINP3 |
| GOTERM_BP_DIRECT | GO:0032332 | Positive regulation of chondrocyte differentiation | 5 | 0.662252 | 2.E-03 | PKDCC, SOX9, SOX6, GLI3, SOX5 |
| GOTERM_BP_DIRECT | GO:0042593 | Glucose homeostasis | 3 | 1.675676 | 2.E-03 | TFAP2B, RFX6, MCU |
| KEGG_PATHWAY | cfa04922 | Glucagon signaling pathway | 6 | 1.6275 | 3.E-02 | PLCH1, ATP6V1G3, ACACA, CREB3L2, SLC2A2, CREB5 |
| KEGG_PATHWAY | cfa04915 | Estrogen signaling pathway | 4 | 2.385 | 2.E-03 | CREB3L2, ADCY8, GRM1, CREB5 |
| KEGG_PATHWAY | cfa04024 | cAMP signaling pathway | 20 | 2.649007 | 1.E-04 | ABCC4, BDNF, CACNA1D, ADRB1, ADCY8, GRIN2B, GLI3, GNAI1, PPP1CC, TIAM1, GRIN3A, PLN, FXYD2, ORAI1, PLCE1, FFAR2, SOX9, CNGA4, GRIA4, RAPGEF4 |
| KEGG_PATHWAY | cfa04918 | Thyroid hormone synthesis | 3 | 2.6375 | 3.E-02 | CREB5, PLCH1, SCAND3 |
| KEGG_PATHWAY | cfa04080 | Neuroactive ligand-receptor interaction | 17 | 2.251656 | 4.E-02 | GABRA2, GABRB2, GRID1, GABRA5, ADRB1, TACR1, NMBR, GRIK2, PRL, GRIN2B, GRM1, GABRG1, GRIN3A, HTR7, GRM4, GRM7, GRIA4 |
Moreover, we have highlighted here the significant enrichment pathways and GO terms that were overlapping among the traits. Among the total enriched GO terms and KEGG pathways, significantly (p < 0.05) enriched 9 GO terms and 7 pathways were found to be shared among traits (Fig. 8). The enriched GO terms are GO:0005634 − nucleus, GO:0047497 − mitochondrion transport along the microtubule, GO:0051569 − Thyroid hormone synthesis GO:0045595 − regulation of cell differentiation, GO:0051561 − positive regulation of mitochondrial calcium ion concentration, GO:0042593 − glucose homeostasis, GO:0050840 − extracellular matrix binding, GO:0015629 − actin cytoskeleton, GO:0005829 − cytosol. The KEGG pathways are cfa04911 − Insulin secretion, cfa04915 − Estrogen signaling pathway, cfa04925 − Aldosterone synthesis and secretion, cfa04922 − Glucagon signaling pathway, cfa04931 − Insulin resistance, cfa04918 − Thyroid hormone synthesis, cfa04022 − cGMP-PKG signaling pathway.