The objective of this study was to investigate the accuracy of imputation to WGS in two pig lines when using a multi-line reference population and the numbers of sequenced animals that belonged to the target lines were small, and subsequently, to investigate the effect of using the resulting imputed WGS genotype or dosage scores in a GWAS. Imputation from 80 K or 660 K to WGS using a multi-line reference population resulted in only 40 to 50% of the SNPs having a Beagle R2 higher than 0.60 in the iWGS data. Nonetheless, when accounting for imputation inaccuracy by filtering iWGS genotypes or by using dosage scores, the number of QTL detected and the estimated proportion of phenotypic variance explained by these were larger compared to when conducting the GWAS using lower density SNP chip genotypes, especially when using iWGS dosage scores. In the following section, the factors that influenced imputation accuracy are discussed, followed by a discussion on the effect of using iWGS genotype or dosage scores on the results of GWAS.
Factors affecting accuracy of imputation
The average Beagle R2 found in this study was 0.39 for the DL-line and 0.49 for the LW-line, which were relatively low compared to other studies that have used multi-line reference populations but with a larger size (more than 242 individuals) [5,6,7], but similar to the accuracy of 0.46 obtained with 90 sequenced Holstein bulls by van Binsbergen et al. [4]. In our study, 168 individuals were included in the reference population and not all haplotypes present in our target populations were represented in the reference population, which increased imputation errors [4, 7]. Van Binsbergen et al. [4] investigated imputation accuracy in three scenarios that differed in the number of animals in a single breed reference population. They showed that for imputation from 50 K to iWGS, the accuracy was 0.37 with a reference population of 45 Holstein cows, and that it increased to 0.46 when the reference population increased to 90 Holstein cows. A similar increase of imputation accuracy was observed when imputing from 777 K to iWGS. These imputation accuracies were comparable or even higher than those found in our study, although we used a larger reference population. However, our reference population consisted of animals from 10 lines and the two target populations were only represented by 12 DL-line animals and 32 LW-line animals. The latter is probably the main reason for the disappointing imputation accuracy obtained in this study.
Our hypothesis was that adding animals from other lines and using a multi-line reference population would improve imputation accuracy, which was demonstrated previously especially for SNPs with a low MAF in the target population but that are segregating in other breeds [4, 6, 31,32,33]. For example, Bouwman and Veerkamp [6] showed that adding 60 individuals from the Jersey, Brown Swiss and Nordic Red Dairy cattle breeds to a reference population of 20 Holstein individuals increased the imputation accuracy of 777 K genotypes to iWGS for Holsteins from 0.71 to 0.83. In our study, the original small reference population was augmented by adding animals from other lines but imputation accuracy was still low for many chromosome regions. The genetic distances of the target populations with the other lines in the reference population were maybe too large because of, for example, different breeding goals. Given the small number of reference animals from the target lines (12 DL-line and 32 LW-line), accuracy of imputation would probably benefit more from additional sequenced DL-line and LW-line animals than from additional animals of other lines.
The imputation accuracy was lower for the DL-line than for the LW-line. Apart from the fact that there were only 12 DL-line individuals in the reference population compared to 32 LW-line individuals, the difference between starting density (i.e. 80 K genotypes for the DL-line and 660 K (imputed) genotypes for the LW-line) and target density (i.e. WGS) also contributed to the lower accuracy for the DL-line versus the LW-line. With a lower density SNP panel, the LD between the SNPs on the genotyping panel and WGS is lower and there is less information to identify shared haplotypes precisely, which increases the uncertainty of imputed genotypes [4, 33]. Stepwise imputation from 80 to 660 K to iWGS, instead of from 80 K to iWGS directly, as suggested by van Binsbergen et al. [4], and as performed for the LW-line, could improve imputation accuracy for the DL-line. However, this was not possible here because the number of available 660 K genotypes for the DL-line was not sufficient.
Other possible reasons for the difference in imputation accuracies between the two lines are differences in population structure or genetic architecture between the lines. For example, selection or different effective population sizes could have resulted in different LD decays and numbers of independent chromosome segments between the lines [34]. Populations with a smaller number of independent chromosome segments are expected to have a higher imputation accuracy because they are expected to have less LD decay across the genome and a greater number of shared haplotypes. In this study, the DL-line had a smaller number of independent chromosomal segments (280.6) than the LW-line (782.3), which was expected to increase imputation accuracy for the DL-Line. However, due to other factors such as the starting density and smaller representation in the reference population, we did not observe a higher imputation accuracy for the DL-Line.
In addition to the above-mentioned factors, imputation accuracy could also be affected by the unequal distribution of SNPs on the genotyping array along the genome and by mapping errors. The latter complicate imputation because incorrect positions of SNPs lead to errors in haplotypes and LD structure in the region they are incorrectly mapped to. Thus, an improved reference genome should increase imputation accuracy [1, 35].
GWAS with iWGS genotype and dosage scores
Iso-Touru et al. [36] and Daetwyler et al. [1] identified new QTL when using iWGS genotype and dosage scores compared to SNP panel genotypes. Similarly, we also found new QTL regions when iWGS genotype scores were used. For the DL-line and the LW-line, 48.9 and 64.4%, respectively, of the QTL detected with iWGS genotype scores were not identified with the 80 K SNP genotypes. The QTL that were identified with all genotype densities were also reported in other GWAS studies that used DL-line or LW-line SNP genotypes [37,38,39,40,41]. In all these studies, the most significant QTL region was located on chromosome 7 and can be linked to the Vert[n]in gene. This gene is important for vertebral development and is positively correlated with number of teats in pigs [37, 39, 41]. Several new QTL regions identified with iWGS genotype and dosage scores include candidate genes for number of teats. For example, at approximately 125 Mb on chromosome 2, a QTL region was identified for the LW-line based on iWGS that includes the casein kinase 1 gamma 3 gene, which plays a regulatory role in the Wnt signalling pathway, which is essential for early mammary gland formation [42,43,44]. For the DL-line, we found a QTL region on chromosome 2 at about 76.5 Mb that was close to another casein kinase 1 gamma gene (casein kinase 1 gamma 2).
Compared to the DL-line, more QTL regions were identified for the LW-line when 660 K genotypes and genotype scores were used but fewer QTL regions were identified when 80 K genotypes or dosage scores were used for GWAS. This was unexpected because the DL-line has a smaller population size and more DL-line iWGS genotype scores were removed because of low imputation accuracy. Therefore, the power to detect associations was expected to be lower in the DL-line than in the LW-line [45]. The larger number of QTL regions identified for the DL-line with 80 K genotypes and dosage scores might be because it has a smaller number of independent chromosomal segments (280.6) compared to the LW-line (782.3). A smaller number of independent chromosomal segments is expected to increase imputation accuracy because of less LD decay across the genome. As a result, power to detect associations might have been higher in the DL-line when using the 80 K genotypes and dosage scores. For 660 K and genotype scores, the power was lower in the DL-line compared to the LW-line because many SNPs with low imputation accuracy were removed. However, another possibility is that the smaller number of independent chromosomal segments for the DL-line, nay have resulted in the identification of a larger number of false positive QTL.
Of the 94 QTL regions found for the DL-line with iWGS genotype scores, 24% overlapped with QTL regions identified for the LW-line with iWGS genotype scores. Although power was limited in each line and we did not expect to detect each QTL in each line, differences might also be caused by differences in genetic architecture between the lines, such as MAF and LD patterns [29]. Comparing the MAF of the most significant SNPs that were identified in the DL-line with the same SNPs in the LW-line (see Additional file 2: Figure S3), clearly showed that the QTL regions had different MAF between the two lines. For example, a SNP on chromosome 8 located at 104.1 Mb had a MAF of 0.005 and 0.36 in the DL-line and LW-line, respectively. Assuming that the SNP is associated with the phenotype in both lines, its low MAF reduces the power to detect it in the DL-line.
Although a large number of new QTL were found using iWGS genotype, whether they are indeed new associations or artefacts of the definition of the QTL region can be questioned. Here, a QTL region was defined as the 0.5-Mb region to the left and right of the most significant SNP in a region because multiple studies have found that the average LD in commercial pig lines decreases below 0.3 when the SNPs are more than 0.5 Mb apart [28,29,30]. This definition may increase the number of QTL regions for iWGS because iWGS consists of many SNPs that are in very high LD. In addition, LD decay can vary across and within chromosomes [29] and, therefore, this definition can be too strict for some QTL regions and neighbouring regions could be merged into one region or the other way around. To test this, we increased the size of the QTL region to 1 Mb to the right and left of the most significant SNP in that region. This resulted in a reduction in the number of QTL detected with iWGS genotype scores from 104 to 64 for the DL-line and from 94 to 56 for the LW-line. Although this is a significant reduction, it is still larger than the number of QTL regions based on the original QTL region definition found with 80 K and 660 K.
In addition to the definition of the QTL region, the newly identified QTL regions could also be an artefact due to the significance threshold that we applied in this study. Here, the significance threshold was set by applying a Bonferroni correction using the number of independent chromosomal segments instead of the commonly used total number of SNPs, which does not take LD between SNPs into account. The number of independent chromosomal segments has often been used for Bonferroni correction for GWAS in human [26, 46, 47], plant [48,49,50], and animal [27, 51] genetics. However, so far there has been little consensus about the most appropriate approach for testing significance for GWAS. This should be a topic for future research.
Finally, even with low imputation accuracy, GWAS using iWGS can be beneficial in several cases. For example, the QTL regions identified in this study could help to pre-select SNPs to improve the accuracy of genomic predictions, as shown by [52,53,54]. In addition, QTL regions identified could be used as indicators for fine-mapping of possible interesting regions, even if they are based on poorly imputed WGS data. For example, SNPs that are in the detected QTL regions could be added to high-density SNP chips, or animals with extreme phenotypes could be sequenced for the detected QTL regions.
Dosage scores
In this study, imputation accuracy by SNP, measured with Beagle R2, was relatively low and, therefore, many SNPs were inaccurately imputed. When using iWGS genotype scores, quality control measures removed around 50% of the SNPs from the analysis and we assumed that the remaining SNPs were called without imputation error. Using iWGS dosage scores in the GWAS model is another way of analysing poorly imputed data. Dosage scores are a posterior probability of having one of the three genotypes and, thus it accounts for the uncertainty of imputation to iWGS. In this study, we identified 56.7 and 26.9% more QTL regions for the DL-line and the LW-line, respectively, when dosage instead of genotype scores were used. Moreover, the most significant SNPs in the QTL regions explained more of the phenotypic variance when using dosage scores. Genotype scores may not capture as much phenotypic variation because some information is lost due to inaccurate imputation. For example, for the DL-line, a clear peak was detected on chromosome 2 at approximately 56 Mb with dosage scores but not with genotype scores (Fig. 7). In this region, many SNPs did not pass the quality control for the iWGS genotype score scenario, because imputation accuracy for this region was low (Fig. 8). In addition, the Manhattan plot based on iWGS genotype scores (Fig. 7) showed a very odd pattern, with an 8-Mb region from 56 to 64 Mb including SNPs that had the same significance level. This region displays high LD (i.e., with an average (± SD) LD r2 of 0.46 (± 0.07)), and thus the SNPs in this regions had approximately the same imputation accuracy and dosage scores, and thus, the same significance level.
Genomic inflation factor
In this study, the genomic inflation factors were lower when using dosage scores compared to using 80 K and 660 K genotypes and iWGS genotype scores. The scenario using dosage scores included SNPs with a low imputation accuracy, and hence SNPs with a low MAF. In the other scenarios, there is ascertainment bias, which is caused by the preference of SNPs on a chip that are more likely to be common [55]. In addition, SNPs selected for iWGS genotype scores have ascertainment bias, because accurately imputed SNPs in general have a higher MAF. By (indirectly) selecting SNPs with a high MAF, it is easier to detect effects for these SNPs, and their surrounding SNPs, and therefore it is likely that more significant SNPs are found than expected based on the theoretical distribution of the test statistic distribution under the null hypothesis. Linkage disequilibrium between the SNPs will also result in more significant SNPs within the region that surrounds a causal variant. Both the higher frequency of SNPs with high MAF and long-range LD increase genomic inflation factors [56,57,58] and the rate of false positives. Use of dosage scores leads to less inflation of the test statistic because the ascertainment bias is partly removed, leading to less biased SNP effects compared to the use of lower density SNP chips and iWGS genotype scores.