Identification of protein-damaging mutations in 10 swine taste receptors and 191 appetite-reward genes

Sequencing statistics

We initially selected 459 kb of target genomic DNA (gDNA) covering the exons from the 12 TASRs and 201 AR genes. After genome capturing, sequencing and read mapping, we successfully covered 372 kb of target sequence with a read depth at each nucleotide position (DP) above 1,000 in the 16 libraries as a whole. This corresponds to 81 % of the initial target size and 201 genes fully or partially sequenced to a DP > 1,000. The poorly sequenced genes include TAS2R3, ENSSSCG000000029894 (a swine ortholog of human TAS2R16), and 10 AR genes. The list of successfully sequenced genes is shown in Additional file 1. After genome capture, sequencing and read mapping, 162,848,637 reads mapped to the target gDNA regions.

Variant identification in TASR and AR genes

Likewise, we sequenced 357,844 (80.1 % of the initial selection) bp of exonic sequence from AR genes at a DP > 1,000, and identified 2,570 variant positions in AR exons, most of which were bi-allelic SNVs. 1,092 (42.4 %) variants were not annotated in dbSNP (Additional file 2). Four of these positions were multi-allelic with each allele predicted to have a different effect on protein sequence. Of the 2,574 assigned variant effects, 25 had the alternative allele fixed in the 16 pools. From the remaining 2,549 variants, 350 were classified as deleterious (H and Mdel) and affected 123 genes (Table 1 and Additional file 2). As for TASRs, H and Mdel as a whole tended to have lower pMAF than Mtol and L (Table 2). We identified 74 H variants in AR genes, 29 of which had pMAFs ≤ 0.01 and are thus considered rare. These H rare variants affected 21 genes (Table 3).

We also plotted the distribution of the AR gene variants predicted to have a strong impact on protein (H and Mdel) and those with a mild impact (Mtol and L) along the protein sequence divided in 10 consecutive position bins of equal amino acid length. We noted that the strong impact variants tended to be more abundant at the end of the protein (bin9 + bin10), and that these also tended to have, on average, larger pMAF. This trend was not observed in the set of mild impact variants (Fig. 1).

Per breed variant distribution

We also compared the percentage of variants that were H or Mdel per breed. No obvious differences were observed (p-val = 0.07) and the Duroc and Landrace were at the top (12.7 %) and bottom (9.4 %) ends, respectively (Table 4).

Altogether, 26 TASR variants were present in a single breed at pAAF ≥ 0.1 and might thus be breed-specific. Not surprisingly, the Asian pool, which involved 15 Chinese Meishan and 7 Vietnamese pigs, showed the highest proportion of pool-specific variants, with 12 being present only in this group (Additional file 3). The wild boar also displayed several allelic particularities, with nine unique variants. Noteworthy, seven of these wild boar variants, all with similar pAAF, mapped to TAS2R1 (Additional file 3). Overall, five (one H and four Mdel) breed-specific variants were predicted to have a strong impact on protein sequence (Additional file 3). The H variant is a stop gain in TAS1R1 that is present in 17 % of the Mangalitza genomes, respectively. We also identified two variants, a synonymous (L) and a non-synonymous tolerated (Mtol) SNVs that affect TAS2R1, that whilst being present at very high frequencies or even fixed in all breeds (pAAF ≥ 0.5), were absent in the Asian pool (Additional file 3).

We detected 306 AR coding variants that were uniquely present or absent in one breed (Additional file 3). Of these, 35 variants involving 32 genes were predicted to be of functional importance (Additional file 3). As in TASRs, the breed specific H and Mdel variants in the AR genes were more abundant in the Asian and the wild boar with 20 and five unique features, respectively. It was noteworthy that two of these H and Mdel variants, both in the Asian pools, were close to fixation (pAAF ≥ 0.9). These variants map to the serotonin receptor, HTR3C, and to the cytochrome P450 gene, CYP2A6 (Additional file 3).

We performed hierarchical clustering using an Unweighted Pair Group Method with Arithmetic Mean (UPGMA) based on the 2,523 (217 TASR and 2,306 AR) variants that were present in at least one of the ten breeds and with pAAF information available in all these populations. One-hundred variants were excluded from the purebred comparison since they were unique to the F₂ animals. The resulting phylogenetic dendrogram is in line with what has been shown in other studies (Fig. 2). Briefly, the Western (European and USA) breeds cluster together and the Asian pool form a separate branch. Within the Western cluster, the Duroc is the only member of a distant branch whilst the European commercial breeds (Large White, Landrace and Pietrain) belong to another sub-group and the more ancient breeds Iberian, Majorcan Black, and Mangalitza cluster together with the wild boar.

pAAF and phenotype relationships in the F₂ groups

Selection and genotype-based validation of a set of protein-damaging variants

With the double aim to validate the subset of variants with more likely damaging impact and, to evaluate by genetic association, their impact in pig breeding, we designed a TaqMan genotyping assay targeting 128 putative polymorphisms. We originally aimed at including all the H variants found in the study and to fill the remaining of the assay with several M variants to cover all TASRs and several AR genes. However, the assay had design requirements that not all the variants fulfilled. Hence, we could only include nine H variants in TASR, 69 H variants in AR genes, 33 M variants in TASRs and 17 M variants in AR genes. To validate the selected variants, we genotyped 237 (Additional file 5) of the 304 sequenced pigs. We could not genotype the remaining 67 pigs as all their DNA was used for the sequencing step. 9 % of the assays failed to either amplify or clearly differentiate the genotype clusters, and 46 % did not show the alternative allele in any sample. 57 variants identified in our variant calling pipeline were validated by the TaqMan assay. The 57 polymorphisms (3 TASR H, 12 AR H, 31 TASR M and 11 AR M) affect eight TASRs and 23 AR genes (Additional file 6). All genotyped variants showed to be in Hardy-Weinberg Equilibrium across breeds (data not shown). Although none of the four genotyped rare polymorphisms (all H) showed a homozygote state for the minor allele, one H SNV (chr16_45548628_G_A), unique to the Vietnamese pigs and with a MAF slightly above the rare variant threshold (MAF = 0.013), was found in the homozygous state in two animals. This variant was predicted to cause a stop codon gain in the very last amino acid of the canonical HTR1A protein. Another introduction of a stop codon (chr8_38038074_G_A) was found in the middle of the Gamma Aminobutyric Acid receptor GABRG1, and was present only in Meishan pigs showing a MAF = 0.022. We next compared the predicted and the matched genotype-observed frequency of the minor allele and detected a strong correlation (r² = 0.93) between the two (Fig. 3).

We also carried a genetic association analysis with the genotypes and phenotypic values for retroperitoneal fat and daily gain in the F₂. 36 of the 57 variants segregated in the F₂ animals. The array contains two of the variants (chr2_3566521_A_C and chr18_7019623_G_A) that showed significant pAAF differences between the F₂ groups with the Fisher test. For daily gain, only the marker chr2_3566521_A_C was nominally significant (p-val = 0.025) (Additional file 7) but it did not reach significance after Bonferroni correction for multiple testing (p-val ≤ 0.05/36 = 0.0013). For retroperitoneal fat, two markers, chr9_46397500_A_G and chr18_7019387_A_T, were nominally significant (p-val = 0.028 and 0.044, respectively) but again, none reached statistical significance after Bonferroni correction (Additional file 7). Chr18_7019623_G_A, which was significant at the pAAF comparison between both F₂ groups, did not reach significance for any of the two traits (p-val = 0.09 and 0.1 for daily gain and retroperitoneal fat, respectively).

Technical considerations

Differences in DNA extraction methods and inaccuracies with respect to quantification might have a negative impact on the even distribution of sequencing reads among the genomes within a pool thereby reducing the accuracy of the pAAF/pMAFs as a proxy of the real allelic frequency. Nonetheless, the comparison of predicted versus observed allele frequencies in the 57 genotyped variants showed that pAAF/pMAFs were very good predictors and supply additional confidence of the quality of our results (Fig. 3).

To the best of our knowledge, this is the first published study from a targeted genome enrichment experiment in swine. As this approach reduces the sequencing throughput requirements, we were able to sequence the target sub-genome (the exome of 201 TASRs and AR genes) in the largest number of pigs (n = 304) reported to date, in a single experiment. By mapping nearly 163 million reads on the target exons, we reached an average depth of sequencing of 72× for each of the haploid genomes. This read depth allowed us to detect all variants present in the pools regardless of their frequency. Given that we sequenced 304 animals harbouring in total 608 chromosome sets (2 alleles each), we were able to detect rare variants with MAF ≥ 0.0016. No previous high throughput sequencing study in pigs reached this sensitivity to detect rare variants. We identified a large number of variant events, most mapping to exon flanking regions including promoters, introns, and upstream and downstream segments. Although some of these non-coding variants could be functional, their regulatory relevance is difficult to predict and was not the subject of our study, which focused mainly, on the detection of protein-damaging coding variation.

Variant identification

Overall, we have identified 2,793 coding variants but 27 of these are fixed in our sample set and are thus likely to be either errors in or private to the reference sequence. We annotated 217 variants segregating in the 10 canonical swine TASRs in our samples. Da Silva et al. [18] described 279 coding variants in their list of 21 taste and nutrient receptor genes after sequencing 79 pigs. However, Da Silva and co-authors included 8 genes that we did not study. These are 5 likely-paralogs of canonical TASRs and 3 genes that are not canonical TASRs but that are related to tasting fat (GPR41, and GPR84) and amino acids (GPRC6A). On the contrary, we included 2 canonical TASRs (TAS2R4 and TAS2R40) not studied by Da Silva et al. Moreover, 4 we classified as AR genes (GPR40, GPR120, CASR, GRM1 and GRM4) that are shared in both studies as they are not canonical TASRs or are receptors belonging to the glutamate pathway. 138 of the 279 variants, 71 in canonical TASR and 67 in the other shared genes, are common to both studies. The difference in the catalogue of identified variants is likely to be mostly due to the fact that Da Silva et al. studied a different set of animals, some belonging to breeds we did not target (Hampshire, creole, Brazilian and Tamworth). 395 of the 2,766 “segregating” variants are predicted to have a high impact (H and Mdel) on the protein sequence of 133 genes and are thus very likely to disrupt or strongly alter their function (Table 1 and Additional file 2). Some of these variants showed allelic frequency relationship with daily gain and retro-peritoneal fat deposition in the F₂ resource with phenotypic records (Table 5 and Additional file 4). In keeping with previous findings in human [13], we also detected that H and M variants tend to be more abundant and have higher pAAFs in TASRs than in non-sensory genes (Tables 1 and 2) which indicates that TASRs are subject to balancing selection. TASRs are comparable to the major histocompatibility complex genes (HLA in humans and SLA in pigs) as both are expressed at the surface of cells to detect particular molecules that could be hazardous for the body. Whilst HLA and SLA detect antigens to promote an immune reaction, TASRs sense chemical compounds to stimulate appropriate responses (reward, acrosomal reaction, smooth muscle contraction or immune function among others) depending on the cell type that is involved. Thus, a healthy animal population needs to be highly polymorphic in these genes in order to maximize its adaptability and survival to variable environments facing multiple threads. As expected, H and Mdel variants had on average, lower pMAFs than Mtol or L variants (Table 2), as the former are more likely to be deleterious and consequently, subject to purifying selection. This is well exemplified by the recent exome sequencing experiments that have successfully identified rare deleterious variants causing rare and severe mendelian disorders in humans [14]. Indeed, H and Mdel variants are more abundant and frequent at the 3′-end of the protein, where they are less likely to have an impact of the function of the protein (Fig. 1).

H rare variants

Remarkably, we found 32 rare H variants which mapped to TAS1R1, TAS1R3 and 20 AR genes (Table 3). These mutations are predicted to fully disrupt the function of the affected gene and might thus have an impact on taste and AR perceptions. Eight of the 20 AR genes carrying rare H variants, have been directly associated to taste, food intake, body size, diabetes or triglyceride levels. Two of them, GRM4 and GABRR1, are glutamate receptors that have been associated to body size [21] and feeding behaviour [22], respectively. Natural knock-out pigs, (i.e., pigs homozygous or compound heterozygous for rare H variants) in these genes may have severe consequences in taste perception and feeding attitude. Reverse genetics experiments to study the phenotypic changes that occur in such natural knock-outs would be very informative to understand the importance of these mutations on feeding behaviour and on a broad range of other traits of relevance in animal breeding and bio-medical sciences. However, the identification of such homozygotes may require the screening of thousands of pigs. A convenient alternative would be to identify heterozygous animals and cross them to generate these natural knock-outs.

Breed particularities

We found several variants with breed particularities that could both explain in part the population history of these breeds and be the result of genetic adaptation to the particular environment or artificial selection. Overall, the pools of commercial breeds (Large White, Landrace, Pietrain, Duroc) contained more variants than the pools of the ancient/traditional breeds (Iberian, Majorcan Black, Mangalitza, Bazna). The sequencing of the commercial breeds involved a larger number of samples and consequently, more genetic variation could be captured. This difference is possibly due also to first, population bottlenecks and a low effective population size within ancient breeds and second, to the introgression of Asian germplasm in the European commercial breeds [23, 24] (Table 4). The Duroc pool is somehow in between European commercial and traditional breeds, with a relative low number of variants (Table 4). However, this was expected as the pig reference genome sequence was obtained from a Duroc animal and is therefore expected to share more similarities with our Duroc pool. On the other hand, the Asian breeds showed the largest number of variants (Table 4) which was also expected as our pool was made of two different swine populations (15 Meishan and 7 Vietnamese pigs) and also, these animals diverged from the European counterparts around one million years ago [25]. It is worth pointing out that as we merged the Meishan and the Vietnamese pigs, we cannot have specific information for these populations. Although the comparison of both Asian breeds would have been interesting, our aim was not to study the genetic variation within these breeds but to identify a large catalogue of variants by sequencing a set of divergent populations. For this reason, and to optimise the high throughput sequencing resources, we merged all Asian animals.

Although we tried to minimise co-ancestry between samples, we do not have pedigree data and hence, we cannot exclude the possibility that the trend on the reduction in the genetic diversity observed in some breeds is due to close familiar relationships. This is particularly true for our Mangalitza samples, which come from a highly inbred population from Romania and might thus not represent the genetic diversity existing within the Mangalitza breed as a whole. To better determine genetic diversity within each breed, we compared the genetic variability in the synonymous sites, which are considered to be neutral in evolutionary terms, in our TASR and AR genes using the method described by Watterson. This method corrects the number of variants by the number of individuals that were analysed [20]. We further corrected these values by the expected number of silent sites in the sequenced cds. Genetic variability among western breeds is highly similar (0.00090–0.00131 variants per mutant site) being lowest in Duroc and highest in Pietrain. The genetic variability seems to be higher in the Asian pool (0.00182).

We also compared the relative abundance of H + Mdel variants in each breed. In contrast to the recent results by Bianco and co-authors [26] who detected a higher ratio of deleterious variants in Western than in Asian breeds, we did not identify large differences across pools. This inconsistency could be explained by the fact that we interrogated a particular set of genes whilst Bianco et al. assessed all annotated genes in the pig genome. Also, both studies screened different animals from distinct breeds. Furthermore, the two studies used different sequencing strategies. Whilst we deep-sequenced 304 pigs in a pool-based strategy, Bianco et al. performed low depth whole-genome shotgun sequencing of 128 pigs. Our approach is better suited for the identification of rare variants and this difference could have altered the catalogue of variants identified in both studies. Therefore, the comparison of the two datasets needs to be taken with caution. Of note, we observed that 25 % of the 28 variants identified in TAS2R1 were wild boar specific (Additional file 3). This could indicate a particular haplotype that might have been lost in the European domestic breeds by artificial selection. Alternatively, as our wild boar gDNA pool was made with samples from three different European locations (Catalonia, Belgium, Romania), these variants, which have a relatively low pAAF [0.16–0.18], could well belong to one of these populations with no germplasm contributed to the analyzed domestic breeds. Noteworthy, some of these TAS2R1 variants are predicted to have a strong effect on the protein and could thus indicate adaptive selection to particular foods or environments (Additional file 3). Especially relevant are the breed-specific H variants with relatively high frequency in the affected breed (pAAF > 0.1) (Additional file 3). The stop codon that prematurely truncates TAS1R1 right in the middle of the protein in 17 % of the Mangalitza genomes might impair the ability to sense umami and might thus affect food preferences in a similar way as described in giant pandas, which lack a functional TAS1R1 potentially disrupting preferences for protein-rich sources [27]. Furthermore, two H variants in HTR3C, and CYP2A6 are unique and almost fixed in the Asian breeds (Additional file 3). As the pAAF between Asian and European pools are so dramatically different, we hypothesize that this might reflect an AR adaptation to very different environments.

Variant and phenotype relationships

In our study, we aimed not only at identifying damaging variants but also at checking whether we were able to detect potential relationships with production traits. This comparison was done with only 38 pigs, which would typically be a very small number of animals and as such, the power to detect genetic associations is low. Thus, it does not aim at finding statistically significant genetic associations but at detecting a trend that could indicate this association. Indeed, we identified several allele - phenotype relationships that could indicate real genetic associations (Table 5). We believe that comparing the pAAFs of the two F₂ pools was a good strategy for the identification of allele - phenotype relationships, as these animals share a common genetic background and the only criteria used to make the pools, was their extreme and opposed phenotypic characteristics. Four of the 11 (36 %) TAS2R4 variants segregating in our F₂ resource showed significant differences at the nominal level in pAAF between the obese and fast growing F2_F pool and the lean slow growing F2_L pool (Table 5). As this is a large percentage, we believe that this could be a real association and that a polymorphism, perhaps a regulatory variant not assessed here, is in part responsible of the phenotypic differences between the divergent groups of pigs. Five of the six genes with significant allelic differences at the nominal level (all with p-val ≤ 0.01) between the F₂ pools have been directly linked to both taste and growth (umami taste and GRM1 [28], fat taste and CD36 [29], taste in general and P2RX2/P2RX7 [30] and weight gain and GABRA6 [31]) (Table 5) and could thus indicate a difference in the eating behaviour. This would explain the difference in growth and fat deposition between the F₂ groups. Remarkably, three variants in CD36, a gene associated to fat taste, had significant pAAF differences between the high and the low fat deposition pigs after correction for multiple testing (Table 5).

Two of the markers, chr2_3566521_A_C and chr18_7019623_G_A, in the ALDH3B1 ortholog ENSSSCG00000026349 and in TAS2R41, respectively, were also genotyped in the F₂ animals with TaqMan probes using the OpenArray technology. Only chr2_3566521_A_C was nominally significant for daily gain (p-val = 0.025), but did not reach the significance threshold after multiple testing (Additional file 7). chr18_7019623_G_A was not significant, but a marker 236 bp upstream from this SNP, chr18_7019387_A_T, was nominally significant (p-val = 0.044) for retroperitoneal fat content (Additional file 7). We seek to screen larger populations with phenotypic data to determine whether these associations could become significant.

Selection and genotype-based validation of SNPs

We also developed a genotyping array to validate the variants with highest damaging potential regardless of their frequency. We believe that a proportion of the 59 non-polymorphic positions might be real low-frequency polymorphisms that were present only in animals that we did not genotype due to the lack of available gDNA. However, we cannot rule out the possibility that some of these variants could simply be false positives caused by the erroneous mapping of some reads to highly similar regions. For example, none of the 17 frame-shift variants we included could be confirmed as a polymorphism but the high false positive rate among indels in high throughput sequencing experiments is well documented. Most of the H variants failed to proof polymorphic but they had very low pMAFs whilst most of the M variants, which on average showed higher pMAFs, were confirmed (Additional file 6). None of the rare H polymorphisms displayed the homozygous state for the minor allele in the genotyped samples. However, we identified two H variants that approached the rare MAF threshold, which showed the three genotypic classes in our animals (Additional file 6). The Vietnamese rare variant in HTR1A is a stop codon that was homozygous in two Vietnamese pigs. Nonetheless, we cannot consider these two animals as natural knockouts as this mutation is located at the very last amino acid of the canonical HTR1A protein, and it is not expected to have an impact on HTR1A function. Therefore, the low frequency of the minor allele might not reflect the existence of purifying selection forces. The Meishan H variant in GABRG1 is also an introduction of a stop codon, but this is located in the middle of the gene and is thus likely to disrupt its function. This gene increases neuronal activity and is associated to eating disorders and anxiety [32] and even with alcohol dependence in humans. Thus, natural GABRG1 knock-out pigs carrying this stop codon may show both signs of anxiety and a particular eating behaviour.

We have developed the first assay for massive genotyping of variants with high functional potential in swine TASR and AR genes, and we now seek to perform association analysis on pig populations with phenotypic records for eating attitude, feed intake, growth, obesity, but also on semen quality and fertility, infection and immunity and behaviour abnormalities such as stress-related stereotypies and tail biting. This approach will help us understanding the impact of genetic variation in TASR and AR genes on traits of interest for the pig breeding industry. This genotyping assay will also allow us to identify natural knock-outs that could then be used in reverse genetics studies. Nonetheless, we acknowledge that this array tags a very small proportion (eight TASRs and 23 reward) of the genes involved in taste, appetite and reward. Ideally, the list of variants and genes should be extended to achieve a comprehensive analysis of the impact of these gene pathways in pig breeding.

Selection of target regions

The genomic regions of TASR exons were selected using the genome annotations accessible via Ensembl’s Biomart (www.ensembl.org/biomart; version 72, June 2013). We selected two TAS1Rs and the 10 annotated TAS2R genes (Additional file 1). TAS1R2 could not be included due to its unknown genomic location at the time of selection.

In total, 166 genes from the AR circuitries were selected from a study aimed at understanding the genetic signatures in giant panda that confer the highly selective diet based on bamboo only [33]. In their study, these authors retrieved the 166 genes from 4 review articles on appetite and food intake behaviour as described in their material and methods. In addition, 35 AR genes were identified by searching NCBI’s PubMed database using the keywords: “appetite”, “food intake”, “dopamine”, “serotonin”, “glutamate receptor”, “epinephrine”, “norepinephrine”, “reward”. We used Ensembl’s Biomart to identify and select the genomic coordinates of the exons from the swine AR orthologous genes (Additional file 1).

Samples

We used gDNA from 266 pigs belonging to eleven breeds or populations including Large White (United Kingdom), Landrace (Denmark), Pietrain (Belgium), Duroc (United States), Iberian (mainland Spain and Portugal), Majorcan Black (Balearic islands), Bazna (Romania), Mangalitza (Hungary), Meishan (China), Vietnamese pot-bellied and to the European wild boar. The Large White, Landrace, Pietrain and Duroc have been subjected to strong selection pressures in the last decades and at present are the most commonly used in intensive production systems. On the contrary, Iberian (mainland Spain and Portugal), Majorcan Black (Balearic islands), Mangalitza (Hungary, although our animals come from a closed herd in Romania) and Bazna (a Romanian breed obtained in the 19th century by crossing Mangalitza and Berkshire) are local breeds with much lower selection intensity and with very localised geographic locations and relatively small productions. The Large White was strongly introgressed with pigs from Asian origin back in the 18th and 19th century. The Landrace was developed in the late 19th century by crossing native Danish breeds with Large White pigs. All samples were collected from farms in Catalonia, Canada and France. Samples were selected based on DNA availability in our DNA archive. Although we do not have pedigree information, we tried to minimise co-ancestry by selecting samples from different farms when possible with the exception of the Mangalitza, which samples come from a close and highly inbred population and are thus likely to have strong familiar relationships. No phenotypic data is available from these samples. Furthermore, DNA from 38 F₂ pigs from an experimental intercross created with the aim of studying obesity-related traits was used. This experimental intercross is described in detail by Kogelman et al. [34]. Briefly, the F₂ population was created by inseminating seven Large White and seven Duroc sows with Göttingen minipig semen from 14 males. Both Large White and Duroc breeds have been selected for leanness and growth traits during many years, while Göttingen minipigs are prone to obesity. 563 F₂ pigs were created. The animals were housed at a regular pig farm, and slaughtered at a commercial slaughterhouse under veterinary supervision. Tissue and blood samples were collected at slaughter. Extensive phenotypic collection was performed from birth to slaughter (242 ± 48 days) including obesity, obesity-related, and metabolic phenotypes; and measurements of fat compartments at slaughter [34]. The retroperitoneal fat compartment was removed from the carcass by blunt dissection and weighted on a bench scale. To calculate daily gain, body weight was measured individually at birth and at 7 month of age (220 ± 45 days). The 38 F₂ pigs used in this study were selected to represent two divergent groups with respect to average daily gain (g/day) and retroperitoneal fat content (kg). The two divergent groups, i.e. the fast growing, obese (F2_F) and the slow growing, lean (F2_L), showed highly significant differences (T-test) for both traits (Additional file 8). gDNAs from purebred animals were extracted using different tissues (blood and a variety of solid tissues) and protocols including standard phenol - chloroform - Isoamyl alcohool organic extraction and the Charge Switch gDNA Micro Tissue kit (Invitrogen). DNA from the F₂ pigs was extracted from EDTA stabilized blood using a salting out procedure. Samples were pooled in 16 tubes on a per-breed basis using semi-equal amounts of gDNA as measured by Nanodrop. The sample size of the pools ranged between 12 and 24 (Additional file 5). Likewise, the sample size of each breed ranged between 12 and 45 (Additional file 5).

The F₂ animals belong to an experimental population in Denmark and they were subject to animal care, maintenance and experimental work according to the ‘Animal Maintenance Act’ (Act 432 dated 09/06/2004) and the approval from the Danish Animal Experimentation Board (J.nr. 2007/561-1434). Specialized professionals at each institution obtained all the other blood and tissue samples following standard routine monitoring procedures and relevant guidelines. No animal experiment has been performed in the scope of this research.

Capture of genomic regions, library prep and high throughput sequencing

The 16 gDNAs pools were subjected to genomic capture and library preparation following Agilent’s SureSelect protocol for Illumina paired-end sequencing. Briefly, three μg of porcine gDNA were sheared on a Covaris™ E220 instrument. The fragment size (150–300 bp) and the quantity were confirmed with the Agilent 2100 Bioanalyzer 1000 chip. Fragmented DNA was end-repaired, adenylated and ligated to Agilent specific paired-end adapters. The DNA with adapter-modified ends was PCR amplified (six cycles, Herculase II fusion DNA polymerase). PCR product size and quantity were determined on the Agilent 2100 Bioanalyzer DNA 1000 assay and hybridized to the genomic capture baits for 24 h at 65 °C (Applied Biosystems 2720 Thermal Cycler). The hybridization mix was washed in the presence of magnetic beads (Dynabeads MyOne Streptavidin T1, Life Technologies) and the eluate was PCR amplified (16 cycles) in order to add the indexed tags using 6 bp SureSelectXT indexes for Illumina. The final library size and concentration was determined on an Agilent 2100 Bioanalyzer 1000 assay.

Each library was sequenced on an Illumina HiSeq 2000 instrument in a fraction of a sequencing lane following the manufacturer’s protocol, with paired end run of 2x101bp. Images analysis, base calling and quality scoring of the run were processed using the manufacturer’s software Real Time Analysis (RTA 1.13.48) and followed by generation of FASTQ sequence files by Illumina’s proprietary CASAVA software.

Read mapping and variant calling

Reads were hard trimmed from the end of the read up to the first base with a quality of at least 10. Reads with at least 40 nt of length were mapped to Sus scrofa reference version 3 (http://hgdownload.cse.ucsc.edu/goldenPath/susScr3/bigZips/susScr3.fa.gz). As the list of TASRs was clearly shorter than the catalogue of AR genes and to minimize type I error (false negative variants), we applied a slightly different read mapping strategy for the two gene sets: 1) For TASR variant detection, reads were mapped first with the GEM toolkit [35] allowing up to four mismatches, and unmapped reads were then aligned to the swine genome using the more permissive BFAST read aligner [36]. We further manually curated the TASR variant list by removing these variants that clustered in high variant density regions as these indicate the presence of wrongly mapped reads - most of them probably aligned by Bfast - and thus false variant calls; 2) For the detection of variants in the AR genes, reads were mapped using only GEM as the manual curation would have been too labour intensive and prone to errors. Alignment (.bam) files containing only properly paired, uniquely mapped reads without duplicates were submitted to variant calling. Each pool was processed separately. The ploidy of the pool was calculated as two times the number of individuals in the pool and input as an optional argument in GATK 3.1 UnifiedGenotyper [37]. For variant calling, read numbers were down-sampled to 1,000 reads per position.

Single pool variant calls were merged into a multi-sample.vcf file using GATK CombineVariants [37]. To confirm that variant positions not called in certain pools had a homozygous reference genotype, a second round of single pool variant calling was performed restricted to the list of variant positions in the merged.vcf file. Subsequently, results were merged again. Functional annotations were added using snpEff [38] with the Sscrofa10.2.69 database, and variants were classified according to their predicted impact as High (H), Moderate (M), Low (L) and Modifier (Additional file 9). Porcine dbSNP version 138 and porcine SIFT scores and deleteriousness prediction were annotated using snpSift [39] and the Ensembl Variant Effect Predictor (VEP) online tool (http://www.ensembl.org/Tools/VEP). M variants were further classified as deleterious (Mdel) or tolerated (Mtol) according to SIFT predictions. Genes of interest and the original target region of the capture experiment were annotated using vcftools [40]. Base counts at variant positions in the merged.vcf file were annotated using GATK Variant Annotator [37]. We used the proportion of reads carrying each allele as an estimator of the allelic frequency for both the alternative allele (pAAF) and for the minor allele (pMAF).

Identification of breed-specific variants

We searched for allelic variants that were uniquely present or uniquely absent in only one breed. We chose these variants that were either breed-specific but with a pAAF in the specific breed ≥ 0.1 and the variants that having a pAAF ≥ 0.5 in all breeds, were absent in only one breed. These breed specific features were assessed on the 2,793 variants identified in TASR and AR genes.

Phylogenetic tree

We calculated the pair-wise Pearson correlation of the pAAFs from the eight European breeds, the wild boars and the Asian animals (Meishan and Vietnamese together) with the 2,523 variants with the alternative allele present in at least one breed and with allelic information in the ten breeds and the wild boar. We then constructed a phylogenetic tree using an UPGMA method. These calculations were done using an in-house developed R script (Do.upgma.pops.bootstrap; https://bioinformatics.cragenomica.es/numgenomics//people/sebas/software/software.html).

Assessment of allele frequency – phenotype relationships

pAAFs at each variant position between the two F₂ pools F2_F and F2_L were compared and the significance of these differences were determined using the Fisher exact test in an R environment. We excluded multi-allelic variants and variants with no read count in at least one pool from the analysis. After this filtering, we were able to compare 1,377 variants.

OpenArray design, genotyping and variant validation

In order to validate the most likely deleterious variants and to develop a genotyping assay to be used in future association studies in pig populations with relevant phenotypes, we developed a genotyping array containing 128 of such potential polymorphisms. We chose the TaqMan® OpenArray® Real-Time technology for Genotyping (Life Technologies) and designed the baits using the online Custom Assay Design tool. This design requires the absence of polymorphisms in the 2 bp window centred at the target variant and has also some constrains on the probe’s melting temperature. As a consequence, we could not design assays for all the variants using this approach. We first selected all the H variants in TASRs and AR genes and supplemented the assay with M variants in both gene groups.

The genotyping was performed in a QuantStudioTM 12 K Flex Real-Time PCR System (Life Technologies). This platform is a high performance, high-throughput technology based on real-time PCR, which enables to run up to 12,000 data points, including SNV, small insertions and deletions (indels), simultaneously.

250 ng of gDNA and master mix were loaded to the OpenArray plates using the AccuFillTM robotic system (Life Technologies), filled with an immersion fluid and sealed. OpenArray plates were genotyped according to the manufacturer’s recommendations. Genotype analysis was performed using both Taqman Genotyper version 1.3 and Symphoni Suite software (Life Technologies).

Genetic association of OpenArray genotypes with daily gain and retroperitoneal fat in the F₂

We used a univariate mixed model to determine genetic associations of the 36 polymorphisms genotyped and segregating in the F₂ animals with the phenotypic records for daily gain weight and retroperitoneal fat content. The analysis was calculated with the software GEMMA [41].