A meta analysis of genome-wide association studies for limb bone lengths in four pig populations

Abstract

Background

Limb bone length is an economically important trait in pigs, because it is negatively correlated with backfat thickness, and is also a determinant to the yield of hip and loin. Moreover, abnormal growth of the limb bone leads to leg structural weakness. Until now, the genetic architecture of the pig lime bone length remains poorly understood. The object of this study was to map genomic loci for limb bone length by genome-wide association study (GWAS) on 4 pig populations.

Results

We measured the lengths of five limb bones including scapula, humerus, ulna, femur and tibia that were dissected from the right-side carcass of 925, 331, 314 and 434 animals from White Duroc × Erhualian F₂ intercross, Erhualian, Laiwu and Sutai populations, respectively. We genotyped the 2004 pigs for 62,163 single nucleotide polymorphisms (SNPs) on the Porcine SNP60 BeadChip, and performed GWAS and a GWAS meta analysis in the 4 populations. In total, we identified 12 and 4 loci associated with the limb bone lengths at suggestive and genome-wide significant levels respectively, of which 4 loci were reported for the first time. The most prominent locus was identified in a 924-kb (kilo base pairs) linkage disequilibrium block on Sus Scrofa chromosome (SSC) 7, and High Mobility Group AT-hook 1 (HMGA1) appears to be a strong candidate gene in this region. Another promising locus is located in the middle of SSC4, and Pleiomorphic Adenoma Gene 1 (PLAG1) is a functionally plausible candidate gene underlying the locus. Because the lengths of the 5 limb bones are highly correlated to each other, most of significant loci were associated with all of the 5 traits; however, several loci showed specific effect on the length of one limb bone, such as the locus at the proximal end of SSC2 associated with only the scapula length.

Conclusion

To our knowledge, this study was the first GWAS meta analysis for limb bone lengths in pigs. As expected, the meta analysis is more powerful to identify genomic loci. A total of 16 loci were identified in this study, including four novel loci. HMGA1 and PLAG1 are two appearing candidate genes for pig limb bone lengths, which warrant further investigations.

Background

Limb bone length is a major factor influencing body height in pigs. It has been reported that body height is negatively correlated with backfat thickness, and is also a determinant to the yield of hip, loin, picnic shoulder and shoulder butt [1, 2]. Abnormal growth of the limb bone could lead to leg structural weakness, and unsound legs would cause huge economic loss, especially in breeding pigs [3].

In humans, body height is positively and highly correlated with limb bone length, with correlation coefficients from 0.66 to 0.88 [4, 5]. The heritability of human height is estimated to be above 0.8, but the top 697 associated SNPs collectively explain only 19.6 % of height heritability in human genome-wide association studies, known as the missing heritability [6, 7]. The domestic pig is an important large-animal model for studying complex traits in human as the pig shares many physiological similarities with human. Deciphering the genetic basis of limb bone growth in pigs will provide insight to understanding of the complex genetic architecture of human height.

Limb bone length is a multi-factorial trait in pigs. To date, a total of 39 quantitative trait loci (QTL) for limb bone length have been identified across the pig genome [8–10]. With the availability of Illumina PorcineSNP60 Beadchip [11], genome-wide association studies have been conducted on a variety of traits to improve the resolution of traditional QTL mapping [12]. GWAS meta analysis based on multiple populations can further increase the detection power and reduce false-positive findings by utilizing information from multiple independent studies [13]. To our knowledge, only one GWAS has been performed for pig limb bone length in an F₂ population produced by a Large White × Minzhu intercross[10].

The aim of this study was to identify genomic loci for the lengths of limb bones including scapula, humerus, ulna, femur and tibia by GWAS and a meta GWAS analysis in four pig experimental populations: a White Duroc × Erhualian F₂ intercross, Chinese Erhualian, Laiwu and Sutai pigs.

Results

Descriptive statistics of phenotypic traits

Descriptive statistics of phenotypic traits measured in the four pig populations are shown in Table 1. In general, the limb bones are longer (P < 0.05) in males than in females across those populations. An additional table shows the phenotypic correlation coefficients among the 5 measured traits related to limb bone length [see Additional file 1]. All of the phenotypic correlation coefficients were positive and highly significant (P < 0.0001).

Table 1 Descriptive statistics and the differences between sexes of the limb bone lengths

Distribution of the Z value in the meta analysis

In this study, we used a Z test to detect the association between SNPs and limb bone lengths in the GWAS meta analysis. The Kolmogorov-Smirnov test results showed that all of the statistics followed the standard normal distribution except the one of the ulna length (P = 0.0448) [see Additional file 2].

SNP quality control results

After quality control, 39461, 28094, 36585 and 41320 SNPs were retained for GWAS on the F₂ (n = 925), Erhualian (n = 331), Laiwu (n = 314) and Sutai (n = 434) populations, respectively [see Additional file 3]. The average physical distances between adjacent SNPs were 74.7, 105.1, 80.7 and 71.2 kb in the four populations, respectively. In the meta GWAS analysis, a common set of 15429 SNPs in the four populations was used, with an average physical distance between adjacent SNPs was 189.9 kb.

Principal component analysis results

The principal component analysis showed that no population stratification existed in the F₂ population, but the other populations had a certain level of population stratifications, especially in Erhualian and Sutai pigs [see Additional file 4]. To remove the effect of population stratifications, a random polygenic effect estimated with a genomic kinship was included in the GWAS model, and the residual inflation was corrected by genomic control. The effect of the population stratification was completely corrected in the four populations [see Additional file 5].

Genetic distances among individuals

The individuals were clustered into four groups according to their genetic distances [see Additional file 6], and the clustering results completely agreed with the populations without exception. The Erhualian and the Sutai were clustered into a group and the rest two populations into another group.

Single-population GWAS results

In single-population analyses, we identified 3 and 4 chromosomal regions associated with limb bone lengths at suggestive and genome-wide significant levels, respectively (Table 2, Fig. 1). In the F₂ we detected two loci associated with all of the 5 traits at 1 % genome-wide significant level. One was around 35 mega base pairs (Mb) on SSC7, and the other was in the middle of X chromosome. For the analysis in the Erhualian, we detected three QTL. One for all of the 5 traits was mapped around 35 Mb on SSC7. Another one for the three front-limb bone lengths was detected around 66 Mb on SSC4. The third one for the ulna length was identified around 103 Mb on the X chromosome. For the Laiwu, only one locus at the distal end of SSC1 was significantly associated with the lengths of scapula and tibia. Another region around 77 Mb on SSC4 had an association with the lengths of femur and ulna almost reaching the suggestive significant level (Fig. 1). For the Sutai, we identified three chromosomal regions affecting limb bone lengths. One for the humerus length was mapped around 9 Mb on SSC3. Another locus around 65 Mb on SSC14 was associated with the scapula length. The third one located around 14 Mb on SSC17 and had a significant effect on ulna length. To our knowledge, there is no QTL reported for limb bone length on this chromosome.

Table 2 The chromosomal regions significantly associated with the limb bones lengths

Fig. 1

Manhattan plots of GWAS for limb bone lengths. The solid, dashed and dotted horizontal lines indicate the 1 % and 5 % genome-wide and chromosome-wise (suggestive) significant threshold values, respectively

Comparing single-population GWAS results

We didn’t identify any significant locus common across all of the four tested populations, but two significant loci each on SSC7 and X were shared by the F₂ and Erhualian populations (Table 2, Fig. 1). This may be due to the fact that Erhualian is one of the two founder breeds of the F₂ population. Another locus around 77 Mb on SSC4 was shared by the Erhualian, F₂ and Laiwu populations (Table 2, Fig. 1 and Additional file 7). Although this locus was not significant in the Laiwu population, its signal was clear and its lead SNP almost reached the suggestive significant level. This locus was also uncovered in the F₂, because it was concealed by the noises from the effects of the major loci on SSC7 and X. After removed those noises, the locus was finally detected.

GWAS meta-analysis results

Under the assumption of the linkage phases being different across populations, we identified a total of 15 chromosomal regions significantly associated with the measured traits. Those loci replicated all of the loci identified by the single-population analysis except the locus associated with the ulna length at the proximal end of SSC3. Four loci were found to be associated with limb bone lengths for the first time, which were the locus at the distal end of SSC8 for the lengths of humerus and ulna, the locus at the distal end of SSC13 for scapula length, the locus at the proximal end of SSC17 and the locus in the middle of SSC18 for the lengths of femur and tibia. We also performed a meta analysis across the 3 pure breeds, and the results were similar to those across 4 populations except that the loci’s significant levels were slightly low [see Additional file 8].

Under the assumption of the linkage phases being same across populations, 8 loci were detected including the three loci on SSC8, 17 and 18 reported for the first time [see Additional file 9], and all of them had been identified under the other assumption.

Chromosomal regions with P < 0.001

Since Bonferroni-corrected thresholds are very conservative, using the Bonferroni correction will decrease the detection power. To overcome it, some loci with significant levels slightly lower (P < 0.001, corrected by genomic control) than the suggestive level were listed in an additional table [see Additional file 10].

GWAS results after fixed the lead SNPs on SSC7 and X in the F₂ population

When we included the lead SNPs on SSC7 and X as fixed effects in the GWAS model, we identified an additional chromosome-wise significant locus at 84.2 Mb on SSC4 and some weak signals on other chromosomes, such as SSC1 and 8 [see Additional file 7].

Linkage disequilibrium block on SSC7

The top SNP located in a 1941-kb linkage disequilibrium (LD) block in the F₂ and in a 924-kb LD block in the Erhualian [see Additional file 11]. The LD block of the F₂ completely encompassed the LD block of the Erhualian, so the locus was fine mapping in a 924 kb interval. Although the four populations shared a 93-kb LD block from 34556148 to 35814316 base pairs, we didn’t guarantee the interval containing the causative gene because there was no signal detected in the other two populations.

Two-point linkage analysis results

There were 20 unmapped SNPs associated with limb bone lengths at the suggestive significant level [see Additional file 12]. Those SNPs consolidated the associating signals on SSC7 in the F₂, on SSC4 and 7 in the Erhualian, and on SSC1, 4 and 7 in the meta analysis. The unmapped SNP ss107830715, which tightly links to the SNP ss107865741 at 15.10 Mb on SSC17, added an associating signal to the tibia length on SSC17 in the meta analysis.

Discussion

In this study, we conducted GWAS and GWAS meta analyses in four pig experimental populations. In total, we identified 16 loci significantly associated with the limb bone lengths, of which 12 confirmed the QTL previously reported on SSC1, 2, 3, 4, 7, 13, 14, 15, and X [8–10], and four loci were reported for the first time.

QTL results vs GWAS results in the F₂ population

Our previous QTL mapping study detected 35 QTL for the five traits on 11 chromosomes in the F₂ population [9]. In this study, we detected only 2 chromosomes associated with those traits in the same F₂ population. The lower power of GWAS could be explained by the following factors. First, the Bonferroni-corrected thresholds in GWAS are more conservative than those obtained from the permutation test used in QTL mapping studies. Second, the error variance is smaller in the QTL mapping study than in the GWAS. When perform genome scan for QTL, the detected QTL in the previous rounds are included in the model as fixed effects, so their variances are removed from the error variance [14]. However, these kinds of variances are not removed from the error variance in GWAS. Third, the statistic model for QTL mapping fixed the additive, dominant and imprinting effects of a locus, while the GWAS model included only the additive effect.

When we used a looser stringent threshold (P < 0.001, corrected by genomic control), we didn’t detected more locus in the F₂ [see Additional file 10]. But when we corrected the effects of the two major loci by including their lead SNPs as fixed effects in the GWAS model, we identified an additional chromosome-wise significant locus on SSC4 and some other weak signals on other chromosomes [see Additional file 7]. This indicates that the significant effects of the loci on SSC7 and SSCX ruined the signals on other chromosomes. Therefore, if a trait is affected by few major loci, corrections of effects of these loci are required to improve the detection power of GWAS.

Single-population GWAS vs GWAS meta analysis

To our knowledge, this is the first GWAS meta analysis of limb bone lengths in pigs. We used a common set of informative SNPs in the four populations to conduct the GWAS meta analysis. Although the SNP density in the GWAS meta analysis is lower than that used in a single population, the meta analysis is more powerful to identify chromosomal regions associated with limb bone lengths than the single-population GWAS. In the single-population GWAS, the detection power is limited by the small population size (In spite there are more than 900 F₂ individuals, the major loci on SSC7 and X ruin the signal on other chromosomes), and the locus with small effect can be unveiled by single-population analysis. The meta analysis utilizes information from multiple independent studies and can increase the detection power. It not only detected 8 more regions than the single population analyses, but also confirmed all loci identified in the four populations separately except the locus at the proximal end of SSC3.

Linkage phases being same vs different across populations

Considering most of SNPs on the porcine 60 K SNP Beadchip are not the causative mutation for the five traits, so it is more reasonable to assume the linkage phases between the SNP and the causative mutation being different than being same across the four populations. Under the “different” assumption, we identified 15 loci for the five traits. While under the “same” assumption, just 8 of the 15 loci were detected and no other locus was identified. Therefore, the former assumption can significantly improve the power.

Comparing the results among traits

As the five traits related to limb bone lengths are highly correlated [see Additional file 1], it is expected to identify chromosomal regions associated with all these traits, such as the region around 35 Mb on SSC7 in the F₂ and Erhualian populations (Table 2). However, some chromosomal regions are associated with one of the five traits. For instance, the region at the proximal end of SSC2 is only associated with the scapula length. This indicates that both common and specific genes are involved in the growths of different limb bones.

Candidate genes underlying promising loci

The locus around 35 Mb on SSC7 showed the most significant effect on the length of all limb bones. This locus has a pleiotropic effect on multiple traits, including the gait scores of front and rear legs [15], fat acid compositions [16], skin thickness [17], feed efficiency [18], growth and fatness [19]. The top SNPs locate in a 924-kb LD block in the F₂ and Erhualian populations [see Additional file 11], and 14 genes locate in this region. Recently, the confidence interval of this multi-faced locus has been refined into a 750-kb region that harbors the high mobility group AT-hook 1 gene (HMGA1) [19], which influences human’s height through chromatin structure modification [20]. Another member of high mobility group AT-hook family HMGA2 was also reported associated with human’s height [21]. Therefore, HMGA1 is a strong candidate gene for limb bone lengths on SSC7.

On chromosome X, the top associated SNPs of the five traits covered a 64.1-Mb region (from 43.4 to 107.5 Mb). It was probably caused by the extremely low recombination rate and low SNP density in this region. This region encompasses a big recombinant cold spot [22], in which SNPs are in highly linkage disequilibrium with each other. Furthermore, the lack of recombination usually results in many discarded SNPs that violate the Hardy-Weinberg equilibrium. The exceptionally large interval makes it impossible to highlight functionally plausible genes for the SSCX locus.

A region around 82 Mb on SSC4 was detected to be associated with the limb bone lengths in the GWAS meta analysis (Table 2) as well as in the F₂ population [see Additional file 7]. Pleiomorphic adenoma gene 1 (PLAG1) locates about 600 kb away from the top SNP and has been reported to be associated with pig body length [23], bovine stature [24] and human height [25], so it is a strong candidate gene for this locus.

Conclusions

In this study, we performed the first GWAS meta analysis for limb bone lengths in pigs. As expected, the meta analysis is more powerful to identify genomic loci than the single population GWAS. A total of 16 loci were identified in this study, including four novel loci. HMGA1 and PLAG1 are highlighted as novel candidate genes for pig limb bone lengths, and they are worthy of further investigations.

Methods

Ethics Statement

All procedures involving in the experimental animal are in compliance with guidelines for the care and use of experimental animals established by the Ministry of Agriculture of China. The ethics committee of Jiangxi Agricultural University specifically approved this study.

Animals and Phenotypic Measurements

A total of 2004 pigs were used in this study including 925, 434, 331 and 314 animals from F₂, Sutai, Erhualian and Laiwu populations respectively. The F₂ and Sutai populations had been described in our previous publications [15, 26]. In brief, the F₂ population was derived from a cross between two White Duroc boars and 17 Erhualian sows [26]. The Sutai pig is a Chinese synthetic breed that is originated from a cross between Duroc boars and Taihu sows, and has been selected for prolificacy and growth more than 18 generations [15]. Erhualian pigs were purchased from Jiaoxi Erhualian Specialized Cooperative Society in Jiangsu province in 2012. To cover most lineages of the Erhualian breed, 331 Erhualian individuals were selected from 11 sire and 53 dam families. Laiwu pigs (n = 314) were bought from Laiwu Stock Seed Farm in Shandong province from 2012 to 2013. Those pigs were derived from 11 sire and 45 dam families and covered almost all lineages of the Laiwu breed. Both Erhualian and Laiwu are Chinese indigenous breeds, the former is famous for its high prolificacy, and the latter is well known for its high intramuscular fat content. Erhualian and Laiwu boars were castrated before weaning, and females were intact. The Erhualian and Laiwu pigs at the age of 60 ± 3 days were transported from their original farms to Nanchang Guohong Ecological Farm in Jiangxi province. During the fattening period, both Erhualian and Laiwu pigs accessed ad libitum to the fresh water and a consistent feed, which met the pig’s needs for energies, amino acids, vitamins and minerals. Both Erhualian and Laiwu pigs were harvested at the age of 300 ± 3 days.

Scapula, humerus, ulna, femur and tibia were dissected from the right side carcass after the animal sacrificed at 240 ± 3 days in F₂ and Sutai populations and at 300 ± 3 days in the other two populations. The lengths of these limb bones were measured according to the method described in Mao et al. [9]. Briefly, scapula length is the maximum straight line distance from the cavitas glenoidalis to the border of scapular cartilage. Humerus length is the length from the head to the trochlea. Ulna length is the distance from the olecranon process to the styloid process. Femur length is the length from the greater trochanter to the intercondyloid fossa. Tibia length is the distance from the intercondylar eminence to the medial malleolus.

Genotyping and Quality Control

Genomic deoxyribonucleic acid (DNA) was extracted from ear or tail tissues under a standard phenol/chloroform protocol. All DNA samples were qualified and diluted to a concentration of 50 ng/μl. All of the animals from the four populations were genotyped by the porcine 60 K SNP Beadchip on an iScan System (Illumina, USA) following the manufacturer’s protocol. The SNP positions were determined according to the pig genome assembly Sscrofa10.2. The unmapped SNPs were assigned to chromosome 0, and their positions were arbitrarily given by the rule described in the reference [15]. SNPs were excluded if their call rate < 95 %, minor allele frequency < 5 %, or severely deviating from Hardy-Weinberg equilibrium (P < 0.000001), and the animals with call rate <95 % were removed. A common set of SNPs that passed quality control across the four populations were used in the GWAS meta analysis.

Statistical methods

Descriptive statistics of the measured traits were calculated by the MEANS procedure. Phenotypic difference between sexes was tested by the TTEST procedure of SAS9.0 (SAS Institute Inc., USA). The CORR procedure was used to compute phenotypic correlation coefficients, and the MIXED procedure was engaged to determine the fixed effects and the covariates for the GWAS model. Sex and batch were included as fixed effects, and body weight at harvest was included as a covariate in the GWAS model for all traits. To remove the effect of the population stratification, a random polygenic effect estimated with a genomic kinship was included in the GWAS model [27, 28]. The genomic kinship was calculated based on the identity-by-state of the SNPs on autosomes [29, 30]. The genomic kinship was also used to estimate the genetic distances among individuals, and a phylogenetic tree was constructed by the heatmap function of the R program.

The GenABEL in the R package was used to conduct individual population GWAS using an additive model [31]. The family-based score test for association was used to detect the association between SNPs and traits, and the residual inflation was corrected by genomic control [32–34]. By applying Bonferroni correction, the P values for suggestive, 5 % and 1 % genome-wide significant levels were one, 0.05 and 0.01 divided by the number of informative SNPs, respectively.

The inverse variance pooling method was applied to perform a GWAS meta analysis of the four populations [35]. The weight (w _i ) for the ith population is calculated by the following formula:

$$ {w}_i=\frac{1}{s_i^2} $$

where s _i is the standard error of the allele effect in the ith population.

Then, the pooled estimate of the effect (β) and its standard error (s ²) are calculated by the following formula:

$$ \begin{array}{l}\beta =\frac{{\displaystyle {\sum}_{i=1}^4{w}_i{\beta}_i}}{{\displaystyle {\sum}_{i=1}^4{w}_i}}\\ {}{s}^2=\frac{1}{{\displaystyle {\sum}_{i=1}^4{w}_i}}\end{array} $$

where β _i is the allele effect in the ith population.

A statistic of Z-test was calculated by the below formula under the assumption of the linkage phases being same across populations, and a Kolmogorov-Smirnov test was used to determine whether the statistic follows the standard normal distribution.

$$ Z=\frac{\beta }{s}=\frac{{\displaystyle {\sum}_{i=1}^4{w}_i{\beta}_i}}{\sqrt{{\displaystyle {\sum}_{i=1}^4{w}_i}}} $$

We used β _i and its absolute value to calculate the pooling β and Z values under the assumptions of the linkage phases being same and different across the populations, respectively. A Bonferroni correction was used to obtain the P values for suggestive, 5 % and 1 % genome-wide significant levels in the meta GWAS analysis.

The linkage between the mapped and unmapped associative SNPs was determined by two-point linkage analysis using CRI-MAP in the F₂ population [36]. Two SNPs were defined as closely linked markers if the logarithm of the odd score was equal to or greater than 3. Haplotypes of target regions were inferred by Simwalk2.9 [37], and linkage disequilibrium blocks were defined using Haploview4.2 [38].