New insight into the SSC8 genetic determination of fatty acid composition in pigs

Background Fat content and fatty acid composition in swine are becoming increasingly studied because of their effect on sensory and nutritional quality of meat. A QTL (quantitative trait locus) for fatty acid composition in backfat was previously detected on porcine chromosome 8 (SSC8) in an Iberian x Landrace F2 intercross. More recently, a genome-wide association study detected the same genomic region for muscle fatty acid composition in an Iberian x Landrace backcross population. ELOVL6, a strong positional candidate gene for this QTL, contains a polymorphism in its promoter region (ELOVL6:c.-533C < T), which is associated with percentage of palmitic and palmitoleic acids in muscle and adipose tissues. Here, a combination of single-marker association and the haplotype-based approach was used to analyze backfat fatty acid composition in 470 animals of an Iberian x Landrace F2 intercross genotyped with 144 SNPs (single nucleotide polymorphisms) distributed along SSC8. Results Two trait-associated SNP regions were identified at 93 Mb and 119 Mb on SSC8. The strongest statistical signals of both regions were observed for palmitoleic acid (C16:1(n-7)) content and C18:0/C16:0 and C18:1(n-7)/C16:1(n-7) elongation ratios. MAML3 and SETD7 are positional candidate genes in the 93 Mb region and two novel microsatellites in MAML3 and nine SNPs in SETD7 were identified. No significant association for the MAML3 microsatellite genotypes was detected. The SETD7:c.700G > T SNP, although statistically significant, was not the strongest signal in this region. In addition, the expression of MAML3 and SETD7 in liver and adipose tissue varied among animals, but no association was detected with the polymorphisms in these genes. In the 119 Mb region, the ELOVL6:c.-533C > T polymorphism showed a strong association with percentage of palmitic and palmitoleic fatty acids and elongation ratios in backfat. Conclusions Our results suggest that the polymorphisms studied in MAML3 and SETD7 are not the causal mutations for the QTL in the 93 Mb region. However, the results for ELOVL6 support the hypothesis that the ELOVL6:c.-533C > T polymorphism has a pleiotropic effect on backfat and intramuscular fatty acid composition and that it has a role in the determination of the QTL in the 119 Mb region.


Background
One of the main sources of human-consumed meat is pork, which represents more than 40% of the meat produced worldwide [1]. The success of pig production is strongly related to improvements in growth and carcass yield. Meat-quality traits are essential for the processing industry and end-consumer acceptance [2], and, as a result, these qualitative traits have been the subject of many studies in breeding programs. Fat content and fatty acid (FA) composition in swine are becoming increasingly studied because of their effect on sensory and nutritional quality of meat. They determine important sensory and technological aspects of pork and meat products because of their influence on the melting point and oxidative status of porcine tissues [3]. In addition, the amount and type of fat in the diet have a major impact on human health. The high consumption of saturated fatty acids (SFA) raises plasma LDL-cholesterol, which is a major risk factor for arteriosclerosis and coronary heart disease (CHD) [4][5][6]. However, recent studies suggest that individual SFA have different physiological effects. Indeed, lauric acid (12:0), myristic acid (14:0) and palmitic acid (16:0) raise LDL and HDL cholesterol plasma levels, whereas stearic acid (C18:0) is considered neutral [7,8], although some epidemiologic evidence suggests that stearic acid (C18:0) is associated with CHD [9]. In contrast, cis-monounsaturated fatty acids (MUFA) and polyunsaturated fatty acids (PUFA) are beneficial for human health. PUFA have been shown to protect against CHD [10], whereas MUFA are also considered to have a hypocholesterolemic effect [11] and, in addition, to have a beneficial effect on insulin sensitivity [12].
The main goals of this work were: (1) to study the QTL architecture for FA composition on SSC8 in the F 2 generation of the IBMAP cross using a panel of 144 informative SNPs, and (2) to analyze additional positional candidate genes.

Animal samples
Animals used in this study belong to the IBMAP experimental population [19]. Two Iberian (Guadyerbas line) boars were crossed with 30 Landrace sows to generate the F 1 generation. Six F 1 boars were coupled with 67 F 1 sows to obtain 470 F 2 animals. In addition, gene-expression analyses were carried out on 56 females from a backcross (BC1_LD) generated by crossing five F 1 (Iberian x Landrace) boars with 23 Landrace sows. All animals were maintained under intensive conditions and feeding was ad libitum with a cereal-based commercial diet. The experiments were performed in Europe following national and institutional guidelines for the ethical use and treatment of animals in experiments. In addition, the protocol was approved by the Ethical Committee of the Institution (IRTA Institut de Recerca i Tecnologia Agroalimentàries). F 2 animals were slaughtered at an average age of 175.5 ± 0.3 days. However, tissues for RNA extraction were not isolated from animals of the F 2 generation. Backcross animals were slaughtered at an average age of 179.8 ± 2.6 days, and samples of liver and adipose tissue were collected, snap-frozen in liquid nitrogen and stored at -80°C until analysis.
Genomic DNA was extracted from blood samples of all animals by the phenol-chloroform method, as described elsewhere.

Traits analyzed
The composition of 10 FA in IMF and BF (taken between the third and the fourth ribs) tissues was determined by gas chromatography as described in [16,17,19]. Subsequently, the percentage of each FA, relative to the total FA, was calculated as well as the global percentages of SFA, MUFA, PUFA and related indices, including desaturation and elongation indices.

Genotyping and quality control
A total of 470 animals were genotyped for 144 SNPs located on SSC8; these include a selection of 142 informative SNPs derived from the Porcine SNP60K BeadChip [20] and two SNPs that corresponded to the previously detected polymorphisms in the FABP2 [21] and MTTP [22] genes. These SNPs [See Additional file 1: Table S1] were included in a custom-generated panel, genotyped using a Veracode Golden Gate Genotyping Kit (Illumina Inc.) and analyzed with a Bead Xpress Reader (Illumina Inc.). SNP positions were based on the whole-genome sequence assembly 10.2 build of Sus scrofa (http://www. animalgenome.org/repository/pig/). All genotypes were assigned using the GenomeStudio software (Illumina Inc.). Markers that had a minor allele frequency (MAF) lower than 5% and missing genotypes that had a frequency greater than 5% were removed using PLINK [23] software. In total, 133 SNPs (92%) passed this quality-threshold filter and were used in the subsequent analysis. Genotypes of all the parents were obtained with the 60 K SNP chip (Illumina) [17] or by pyrosequencing [21,22].
Fifty-six animals of the BC1_LD were genotyped for SNPs SETD7:c.-1034T > G and SETD7:c.700G > T and the two MAML3 microsatellites for gene-expression studies. In addition, a subset of 168 F 2 animals were genotyped for SNPs ELOVL6:c.-533C > T, SETD7:c.700G > T and the two MAML3 microsatellites for association studies. All parents and grandparents of these animals were also genotyped in the same way.

Association analysis
Association analysis was performed for FA composition and indices of FA metabolism in 470 F 2 animals. A mixed model that accounts for additive effects was performed using Qxpak 5.0 [24]: where y ijlkm is the l th individual record, sex (two levels) and batch (five levels) are fixed effects, β is a covariate coefficient with c being carcass weight, λ l is a -1, 0, +1 indicator variable depending on the l th individual genotype for the k th SNP, a k represents the additive effect associated with SNP, u l represents the infinitesimal genetic effect treated as random and distributed as N(0, Aσ u ) where A is a numerator of the kinship matrix and e ijlkm is the residual. A similar model that fitted different QTL effects was used to test the hypothesis of the presence of two QTL located in the studied regions with effects a 1 and a 2 on the same FA: The R package q-value [25] was used to calculate the false-discovery rate (FDR), and the cut-off of the significant association at the whole-genome level was set at the q-value ≤ 0.05. Version 2.15.2 of R [26] was used to calculate the descriptive statistics for the 10 analyzed traits and their related indices.
For linkage and linkage disequilibrium (LDLA) analysis, haplotypes were reconstructed using DualPHASE software [27], which exploits population (linkage disequilibrium) and family information (Mendelian segregation and linkage) in a Hidden Markov Model setting. Then, QTL fine-mapping was performed for the most significant traits C16:1(n-7), C18:0/C16:0, C16:1(n-7)/ C18:1(n-7), and the FA average chain-length (ACL) by applying the mixed model: in which b is a vector of fixed effects (sex and batch), h is the vector of random QTL effects corresponding to the K cluster defined by the Hidden State, u is the vector of random individual polygenic effects and e is the vector of individual error.
Amplification and sequencing of the pig MAML3 and SETD7 genes Genomic DNA samples from 10 individuals of the BC1_LD and two Iberian boars were used to amplify and sequence the proximal promoter and exon 1 of the MAML3 and SETD7 genes.
A 931-bp region of the MAML3 gene was amplified and sequenced in two overlapping fragments of 517 bp and 663 bp. Primers [See Additional file 2: All primers were designed using the PRIMER3 software [28] and were validated using the PrimerExpress 2.0 software (Applied Biosystems).
PCR (polymerase chain reactions) were carried out in a total volume of 25 μL containing 0.6 units of Ampli-Taq Gold (Applied Biosystems), 1.5 to 2.5 mM MgCl 2 depending on the primers [See Additional file 2: Table  S2], 0.2 mM of each dNTP, 0.5 μM of each primer and 20 ng of genomic DNA. The temperature profile was 94°C for 10 min and 35 cycles at 94°C for 1 min, 58°C to 62°C depending on the primers [See Additional file 2: Table S2] for 1 min and 72°C for 1.5 min, including a final step of 7 min at 72°C. Gradient parameters were determined based on size and GC content of the amplicon. The samples were then analyzed on 1.5% agarose gels. Purification was performed using an Exonuclease I and FastAP™ Thermosensitive Alkaline Phosphatase [29]. For the sequencing reaction, we used the Big Dye Terminator v.3.1 Cycle Sequencing Kit and an ABI Prism 3730 DNA analyzer was employed (Applied Biosystems). Polymorphisms were checked through the Seq scape v2.1.1 program (Applied Biosystems).

Detection of microsatellite polymorphisms
Based on the sequencing results of the promoter region and exon 1 of the MAML3 gene, two new microsatellites were identified. Both microsatellites were independently amplified using fluorescent primers [See Additional file 2: Table S2]. PCR were performed in a 25-μL reaction mix containing 20 ng of genomic DNA, 0.2 mM of each dNTP, 2.5 mM MgCl 2 , 0.5 μM of each PCR primer and 0.6 units of AmpliTaq Gold (Applied Biosystems). PCR were run as follows: 94°C for 10 min, 35 cycles of 94°C for 1 min, 58°C for 1 min, 72°C for 1.5 min and a final extension step at 72°C for 7 min. The two amplicons were mixed at a ratio of 1:3 (HEX:FAM) and analyzed using capillary electrophoresis on an ABI Prism 3730 DNA analyzer (Applied Biosystems) and the ROX-500 GeneScan Size Standard. The peak height of each product was determined using Peak Scanner 2 software (Applied Biosystems).

RNA isolation and cDNA synthesis
Total RNA was extracted from liver and BF tissues using the RiboPure kit (Ambion), according to the manufacturer's recommendations. RNA was then quantified using a NanoDrop ND-1000 spectrophotometer (NanoDrop products) and RNA integrity was assessed with an Agilent Bioanalyzer-2100 (Agilent Technologies). One μg of total RNA of each sample was reverse-transcribed using the High-Capacity cDNA Reverse Transcription kit (Applied Biosystems) in a reaction volume of 20 μL.

Gene-expression quantification
Fifty-six females of the BC1_LD were used to quantify gene expression. The expression of MAML3 and SETD7 was analyzed using the 48.48 microfluidic dynamic array IFC chip (Fluidigm) according to the manufacturer's instructions. Briefly, 2 μL of 1:5 diluted cDNA was preamplified using 2X Taqman PreAmp Master Mix (Applied Biosystems) and 50 nM of each primer pair in 5-μL reaction volume. The cycling program consisted of an initial step of 10 min at 95°C followed by 16 cycles of 15 s at 95°C and 4 min at 60°C. At the end of this pre-amplification step, the reaction products were diluted 1:5 (diluted preamplification samples). RT-qPCR on the dynamic array chips was conducted on the BioMarkTM system (Fluidigm). A 5-μL pre-mix sample containing 2.5 μL of Sso-Fast EvaGreen Supermix with Low ROX (Bio-Rad), 0.25 μL of DNA Binding Dye Sample Loading Reagent (Fluidigm) and 2.25 μL of diluted pre-amplification samples (1:16 or 1:64 from the diluted pre-amplification samples from liver and BF samples, respectively), as well as a 5-μL assay mix containing 2.5 μL of Assay Loading Reagent (Fluidigm), 2.25 μL of DNA Suspension Buffer (Teknova) and 0.25 μL of 100 μM primer pairs (500 nM in the final reaction) were mixed inside the chip using the IFC controller MX (Fluidigm). The cycling program consisted of an initial step of 60 s at 95°C followed by 30 cycles of 5 s at 96°C and 20 s at 60°C. A dissociation curve was also drawn for each primer pair.
Data were collected using the Fluidigm Real-Time PCR analysis software 3.0.2 (Fluidigm) and analyzed with the DAG expression software 1.0.4.11 [30] using standard curves for relative quantification. Relative standardcurves with a four-fold dilutions series (1/4, 1/16, 1/64, 1/256, 1/1024) of a pool of 10 cDNA samples were constructed for each gene to extrapolate the value of the quantities of each studied sample. Of the four endogenous genes tested (ACTB, B2M, HPRT1, TBP), ACTB and TBP had the most stable expression [31] in both tissues. The normalized quantity values of each sample and assay were used to compare our data.
Mean values between genotypes were compared using a linear model implemented in R, which performs a single stratum analysis of variance considering sex and batch as fixed effects. Differences were considered statistically significant at a p-value of 0.05.
Two regions that contain trait-associated SNPs (TAS) were clearly visualized in the association plots at around 93 Mb and 119 Mb for all of the above-mentioned FA and indices with the exception of the C20:2(n-6)/C18:2 (n-6) elongation ratio, which showed only one significant TAS region at 120.99 Mb (Table 1). For all significant traits, the 119 Mb TAS region showed a stronger signal than the 93 Mb region. The strongest effects of both TAS regions were found for palmitoleic acid (C16:1 (n-7)) content and C18:0/C16:0 and C18:1(n-7)/C16:1 (n-7) elongation ratios.
A combination of linkage disequilibrium and linkage analysis (LDLA) was then performed for the most significantly associated traits [See Additional file 2: Table S3]. With this haplotype-based approach, it is possible to simultaneously exploit linkage analysis and linkage disequilibrium. Several studies have shown the usefulness of this strategy for fine-mapping and QTL interval reduction [27,32]. The LDLA study identified the two TAS regions by association analysis, with the 119 Mb region showing the strongest statistical signal for all analyzed traits. Figure 1 shows the two genomic regions identified for the C18:0/C16:0 elongation ratio. Plots of the other three traits analyzed are shown in Additional file 3: Figure S1.
In order to determine whether one or two QTL were segregating on SSC8 for the BF FA and their indices, models fitting one QTL against a model considering two different QTL were tested. Results of the LR test indicated that the model with two QTL was the most likely for the 10 traits analyzed [See Additional file 2: Table S4].
Previously, a QTL scan for BF FA composition was performed with 369 animals from the same F 2 generation [16], but only six microsatellite markers were genotyped. A clear effect of SSC8 markers was observed only for percentages of palmitic (C16:0) and palmitoleic (C16:1(n-7)) FA and ACL. A suggestive effect on percentage of oleic acid (C18:1(n-9)) was also observed. However, the confidence interval for this QTL was greater than 30 cM. Two other studies of our group analyzed positional candidate genes for this QTL, i.e. FABP2 [21] and MTTP [22], but the localization of the QTL was not refined. In addition, QTL for IMF palmitic (C16:0) FA composition have been reported in a Duroc x Large White F 2 cross [13] and for stearic (C18:0) FA in a White Duroc x Erhualian F 2 cross [14].
Here, two QTL at approximately 93 Mb and 119 Mb were detected and affected the BF composition of the six FA and the four indices mentioned above in the 470 F 2 animals. The palmitoleic (C16:1(n-7)) FA QTL on SSC8 have been shown to be segregating in different crosses of the IBMAP population, and both QTL have a pleiotropic effect on BF and IMF FA deposits.

Gene annotation and identification of polymorphisms in positional candidate genes
Gene annotation of the two TAS genomic regions allowed us to identify genes related to FA metabolism. In the first region, the genes mastermind-  MAML3 is a member of the Mam gene family, which plays an essential role in the stabilization of Notch transcriptional activation complexes [33]. This Notch signaling pathway mediates short-range communication between cells, and it has recently been associated with the regulation of lipogenesis and gluconeogenesis in liver [34]. A 931-bp fragment of the pig MAML3 gene that covers part of the promoter region and part of exon 1, was amplified from genomic DNA and sequenced. Two novel microsatellites were found: MAML3_MS1, a (CA) n tandem repeat located in the promoter region and MAML3_MS2, a (CGG) n tandem repeat identified in exon 1. The variability of both microsatellites is described in Table 2.
The product of the SETD7 gene is a histone methyltransferase that specifically monomethylates Lys-4 of histone H3 [35] and, thus, it is involved in the epigenetic transcriptional regulation of genes, activating genes such as collagenase or insulin [36]. To identify polymorphims in the porcine SETD7 gene, a 839-bp fragment of the SETD7 promoter and exon 1 was amplified from genomic DNA and sequenced. In addition, the identification of polymorphisms in the entire coding region of the SETD7 gene was performed using RNA-Seq data [37] with the Integrative Genomics Viewer (IGV) software (http://www.broadinstitute.org/igv/). Alignment and analysis of these sequences led to the identification of nine polymorphisms ( Table 3). Two of these polymorphisms were used to genotype BC1_LD animals, one located in the promoter region (SETD7:c.-1034T > G) and one nonsynonymous polymorphism in exon 6 (SETD7:c.700G > T), which determines an amino-acid change of valine to leucine. Apart from the fact that these SNPs are located in the SETD7 gene, they were selected because they showed divergent allelic frequencies between the Iberian and Landrace IBMAP founders i.e. the SETD7:c.-1034 T and SETD7:c.700 T alleles were fixed in the Iberian boars. Complete linkage disequilibrium between the two SNPs was observed in the genotyped BC1_LD animals and, thus, only SETD7:c.700G > T was further genotyped in 168 animals belonging to the F 2 generation.  In the second region, the ELOVL6 gene was identified at position 120.12 Mb. The ELOVL6 gene is a strong positional and functional candidate gene involved in de novo lipogenesis and acts on the elongation of SFA and MUFA. A polymorphism in the promoter region of this gene (ELOVL6:c.-533C > T) has previously been associated with percentages of palmitic and palmitoleic FA in muscle and backfat in the BC1_LD population [38]. In addition, expression of the ELOVL6 gene was lower in the backfat of animals with the Iberian allele in comparison to those with the Landrace allele. As expected from the elongation function of this gene, a lower ELOVL6 expression was associated with a higher percentage of palmitic and palmitoleic FA in muscle and adipose tissue [38].
Regarding oleic FA (C18:1(n-9)), the main dietary FA, its content in BF was decreased in animals with the Iberian allele for the TAS regions on SSC8. It must be noted that the opposite effect was observed in the major SSC4 and SSC6 TAS regions for oleic (C18:1(n-9)) IMF content [17].
Effect of the SETD7:c.1034T > G and SETD7:c.700G > T SNPs and MAML3 microsatellites on gene expression The expression profiles of the pig MAML3 and SETD7 genes were studied in liver and BF tissues of 56 BC1_LD females by RT-qPCR. The analysis of gene expression in the F 2 generation was not possible because tissues for RNA isolation were not available. Differences in the expression of MAML3 among animals were observed, with  coefficients of variation (CV) of 35% and 42% in liver and BF, respectively. The SETD7 gene expression was less variable, with CV values of 18% and 33% in liver and BF tissue, respectively. However, no significant differences in expression of SETD7 were detected among animals classified according to the SETD7 genotypes (either SETD7:c.-1034T > G or SETD7:c.700G > T) in either tissue. Similarly, no differences in expression of MAML3 were observed among animals classified according to the MAML3_MS1 and MS2 microsatellites. In addition, no significant correlation was found between expression levels of MAML3 or SETD7 in the liver and adipose, which suggests that different and tissue-specific mechanisms control the liver and adipose tissue expression of MAML3 and SETD7.
Association study for BF FA composition with markers located in positional candidate genes Two microsatellites in the MAML3 gene (MAML3_MS1 and MAML3_MS2), one SNP in the SETD7 gene (SETD7: c.700G > T), and one SNP in the ELOVL6 gene (ELOVL6: c.-533C > T) were used to genotype 168 animals of the F 2 generation. An association analysis with these markers and the SSC8 genotypes from the 133 SNPs of our custom porcine SNP panel was performed.
For the first region (93 Mb), polymorphisms in the SETD7 and MAML3 genes were studied. For SETD7, the SETD7:c.700G > T polymorphism did not show the most significant association (Table 4). In addition, MAML3 gene microsatellites showed no significant associations for any of the traits studied. However, the SNPs showing the strongest signals (Table 4) were located within a 2 Mb interval of the SETD7 and MAML3 genes. These results suggest that other non-genotyped polymorphisms may cause the observed effects on FA composition in the 93-Mb region.
For the second region (119 Mb), a polymorphism in the ELOVL6 gene was studied. The ELOVL6:c.-533C > T polymorphism showed the highest association with percentage of palmitic and palmitoleic FA, ACL, and C18:0/ C16:0 and C18:1(n-7)/C16:1(n-7) ratios (Table 4). Hence, these results are consistent with those found in the IMF FA composition of the BC1_LD generation [38]. The clear association of the ELOVL6:c.-533C > T polymorphism with percentage of FA in IMF and BF indicates a pleiotropic effect of this gene in both tissues.
Analysis of the additive value of SNPs SETD7:c.700G > T and ELOVL6:c.-533C > T showed a higher contribution of ELOVL6:c.-533C > T SNP for all studied FA and indices. Furthermore, the additive value of the two SNPs [See Additional file 2: Table S5] showed an effect in the same direction. For instance, the Iberian alleles of both QTL increased palmitic and palmitoleic FA content and reduced the elongation ratios. These results are in accordance with the reported Iberian-Landrace breed differences in BF FA composition [39].