Expression-based GWAS identifies variants, gene interactions and key regulators affecting intramuscular fatty acid content and composition in porcine meat

The aim of this work is to better understand the genetic mechanisms determining two complex traits affecting porcine meat quality: intramuscular fat (IMF) content and its fatty acid (FA) composition. With this purpose, expression Genome-Wide Association Study (eGWAS) of 45 lipid-related genes associated with meat quality traits in swine muscle (Longissimus dorsi) of 114 Iberian × Landrace backcross animals was performed. The eGWAS identified 241 SNPs associated with 11 genes: ACSM5, CROT, FABP3, FOS, HIF1AN, IGF2, MGLL, NCOA1, PIK3R1, PLA2G12A and PPARA. Three expression Quantitative Trait Loci (eQTLs) for IGF2, ACSM5 and MGLL were identified, showing cis-acting effects, whereas 16 eQTLs had trans regulatory effects. A polymorphism in the ACSM5 promoter region associated with its expression was identified. In addition, strong candidate genes regulating ACSM5, FOS, PPARA, PIK3R1, PLA2G12A and HIF1AN gene expression were also seen. Notably, the analysis highlighted the NR3C1 transcription factor as a strong candidate gene involved in the regulation of the 45 genes analysed. Finally, the IGF2, MGLL, MC2R, ARHGAP6, and NR3C1 genes were identified as potential regulators co-localizing within QTLs for fatness and growth traits in the IBMAP population. The results obtained increase our knowledge in the functional regulatory mechanisms involved in these complex traits.


Results and Discussion
Gene-expression profiling. In the present study, 114 BC1_LD animals generated from the experimental IBMAP population and showing a wide range of IMF content and FA composition values were used for RNA extraction and RT-qPCR to perform eGWAS (Supplementary Table S1).
A significant sex effect (p-value ≤ 0.05) between gene-expression levels was detected in 20 out of the 45 genes analysed (44%): ACSS1, ACSS2, ATF3, CREG1, DGAT2, ETS1, FABP5, HIF1AN, IGF2, NCOA2, NCOA6, PLA2G12A, PPARA, PPARG, PPARGC1A, PRKAA1, PEX2, SCD, SP1 and SREBF1 (Fig. 1). Interestingly, some of the sex-biased genes identified are key regulators of lipid metabolism, such as PPARA, PPARG, PPARGC1A and SREBF1. Note that within the sex-biased genes, several lipogenic genes were more expressed in females (DGAT2, NCOA2, NCOA6, PPARG, PRKAA1, SCD, SP1 and SREBF1), whereas more lipolytic genes were over-expressed in males (ATF3, PPARA and PPARGC1A). Sexually dimorphic gene expression of genes involved in lipid metabolism has already been described in muscle as well as in other tissues such as liver 19,20 . Consistent with our results, PPARA and SREBF1 presented sex-biased gene expression in human skeletal muscle 21 . Furthermore, PPARGC1A and PPARA gene expression in muscle are affected under diet supplementation with 17β -estradiol (E 2 ) 21 . In addition, the SCD and PPARG mRNA expressions have been observed to be lower in male than in female human muscle cell cultures 22 . To assess the relationship between muscle gene-expression levels and phenotypes, a hierarchical cluster analysis of the correlation values between the gene expression levels of the 45 genes and the fatty acid content in muscle was performed (Fig. 2). The hierarchical cluster analysis showed that genes mainly related to lipogenic pathways (DGAT2, PPARG, SCD, MGLL, NCOA1, NCOA2, NCOA6, PRKAA1, SP1, SREBF1) clustered together, whereas genes mainly related to lipolytic pathways (ATF3, CPT1B, PPARA, PPARD, PPARGC1A) joined together (Fig. 2). Genes within the lipogenic-related cluster showed, in general, positive correlations with palmitoleic (C16:1(n-7)) and octadecenoic (C18:1(n-7)) FAs, while the lipolytic-related group showed a positive correlation with PUFAs in general, specifically with the linoleic (C18:2(n-6)) FA. Overall, these results are in agreement with a previous muscle RNA-Seq transcriptome study of animals which were extreme for IMFA composition performed by our group 16 , where a higher expression of genes related to lipogenic pathways was observed in the muscle of animals with high MUFA and SFA content in muscle.
From the associated SNPs (n = 241), 215 SNPs were successfully annotated with the Variant Effect Predictor (VEP) of Ensembl (Sscrofa 10.2 annotation release 80) of which 54% (117 SNPs) were located in intergenic regions. The remaining 46% (98) of SNPs were mapped within 76 genes: 69 (32%) in introns, 10 in the 5′ flanking region, 13 in the 3′ flanking region, three in the 3′ UTR region, two in the coding regions of genes and determining synonymous mutations, and one missense mutation (Supplementary Table S2). Twelve out of the 98 intragenic SNPs were located within a gene exerting a lipid metabolism function (ABCA3, ACSM3, LRP5, LHFPL4, BRPF1,  HRH1, PPARG, ACAD9, COPG1, MGLL, ACAD11, ARHGAP26) (Supplementary Table S2), but only six of them were expressed in muscle (ACSM3, PPARG, ACAD9, MGLL, ACAD11 and ARHGAP26), according to the Gene Expression Atlas 23 and the RNA-Seq experiment performed in individuals from the BC1_LD 16 (Supplementary  Table S2). Only one of the intragenic SNPs mapping within genes with a lipid-related function was identified in a trans eQTL: the ARHGAP26 gene located in a trans-eQTL for the PPARA gene. This gene is activated via lipid interaction, however, its role has not been well-defined 24 .
cis-eQTL identified. The Insulin-Like Growth Factor 2 (IGF2) gene is not mapped in the current Sscrofa10.2 assembly 25 . However, it has been located by linkage map in the telomeric end of the p arm of SSC2 26 . An intronic IGF2 mutation (IGF2-intron3-G3072A) was described by Van Laere et al. 27 to have a major effect on muscle growth in pigs. Although this mutation only segregated in a small family in the IBMAP F 2 population 28 , the cis-eQTL for IGF2 suggests that the intronic mutation identified by Van Laere et al. 27 may segregate in the IBMAP BC1_LD, having an important effect on the traits analysed. However, after genotyping the IGF2-intron3-G3072A mutation in the BC1_LD animals, we observed that this SNP was not the most associated (p-value = 1.07 × 10 −7 )  26 .
with the IGF2 mRNA expression levels (Fig. 3B), suggesting that other mutations in the IGF2 gene may be the responsible for the mRNA variation in the BC1_LD population. Further analysis will be conducted to validate our hypothesis. Concerning the Monoglyceride Lipase (MGLL) eGWAS results, one of the annotated cis-SNPs (ASGA0103932) was mapped within an intronic region of the MGLL gene (Supplementary Table S2). Nevertheless, this SNP was not the most significant, associated SNP (ASGA0103932; p-value = 2.34 × 10 −8 ), suggesting the presence of another polymorphism within or near this gene as the causative mutation affecting MGLL gene expression levels. The most significant cis-SNP for MGLL (ASGA0093606; p-value = 2.20 × 10 −9 ) was located less than 0.69 Mb downstream of the MGLL gene. In addition, an SNP (ISU10000701) in the MGLL cis-eQTL was annotated in the upstream region of the PPARG gene, more than 4 Mb from the MGLL gene (Supplementary Table S2). PPARG and MGLL were reported to be co-associated in a previous study of our group for growth and fatness traits, where PPARG was described as a major regulator 17 . Additionally, the literature-based analysis with Genomatix also identified an interaction between these two genes ( Supplementary Fig. S1), and a chromatin immunoprecipitation (ChiP) experiment performed in epithelial cultured cells revealed PPARG binding sites in the distal MGLL promoter 29 . Therefore, it is likely that this region, initially considered as a single cis-acting eQTL, in fact comprises two eQTLs. The correlation between PPARG and MGLL gene expression in muscle was moderate (r = 0.48; p-value = 8.28 × 10 −8 ), reinforcing the hypothesis that at least two factors may determine the different levels of MGLL mRNA. In this regard, the MGLL gene expression would be affected by a variant present in the same gene (MGLL), and also by PPARG. However, both markers (ASGA0093606 and ISU10000701) were in strong LD (R 2 = 0.99), not being able to assess whether MGLL gene expression is affected by cis-, trans-or both cis-/trans-factors on SSC13. Further investigation with other populations is required to reveal this. Remarkably, the significantly associated SNP inside the PPARG gene (ISU10000701) was one of the three main co-associated SNPs in the network identified in Puig-Oliveras et al. 17 involved in the determination of growth and fatness traits. Nonetheless, no association was found between the SNP within PPARG (ISU10000701) and the PPARG gene expression.
Finally, three SNPs (ASGA0090088, ASGA0105223 and SIRI0001454) in complete linkage disequilibrium (R 2 = 1) were the most significant ones in the cis-eQTL (p-value = 7.12 × 10 −14 ) associated with the mRNA levels of the Acyl-CoA Synthetase Medium-Chain Family Member 5 (ACSM5) ( Table 1). The ASGA0090088 marker was the closest cis-SNP, mapping at approximately 798 kb from the upstream region of the ACSM5 gene (Supplementary Table S2). Notably, after amplifying and sequencing the proximal promoter region of the ACSM5 gene in ten BC1_LD animals, three polymorphisms were identified at positions g26260919C > T (rs323520560), g26260505G > A (rs339587799) and g26260422G > A (rs331702081), according to the Ensembl ENSSSCG00000026453 sequence, which seem to be in linkage disequilibrium. The most proximal 5′ mutation (g26260422G > A) was genotyped in the BC1_LD animals, and the association analysis of this polymorphism with the ACSM5 mRNA expression values revealed that this SNP was the most associated one (p-value = 7.11 × 10 −15 ; Fig. 3A). Further analyses are necessary to determine the role of these polymorphisms in the transcriptional regulation of the ACSM5 gene.
trans-eQTL identified. Table 2 summarizes all the relevant lipid-related genes mapped in the trans-eQTL regions associated with the genes analysed.
ACSM5 gene expression was also associated in trans with two chromosomal regions on SSC3 and SSC10. The most associated SNP on SSC10 was an intronic polymorphism (ASGA0090778; p-value = 2.55 × 10 −8 ) in the Component of Oligomeric Golgi Complex 7 (COG7) gene (Supplementary Table S2). Although the COG7 gene has not been described to have a direct relationship with lipid metabolism function, it is reported that members of the COG complex are involved in intra-golgi trafficking and glycosylation of proteins and lipids 30 Table 2. Trans-eQTL annotation.
where lipids and lipid hormones bind to be exchanged between biological membranes 31 . Supporting STARD13 as a regulator of FOS gene expression, transgenic mice with pancreas STARD13 ablation showed no detectable mRNA expression of the FOS gene 32 . Two regions on SSC2 and SSC6 were associated in trans with the Peroxisome Proliferator-Activated Receptor Alpha (PPARA) gene expression. In SSC2, several genes related to lipid metabolism were identified: Rho GTPase Activating Protein 26 (ARHGAP26), Fibroblast Growth Factor 1 (Acidic) (FGF1) and Nuclear Receptor Subfamily 3, Group C, Member 1 (glucocorticoid receptor) (NR3C1) (Supplementary Table S3). ARHGAP26 is activated via lipid interaction 24 and may play a role in adipogenesis, as well as the transcription factor NR3C1, which maps at approximately 0.5 Mb of the ARHGAP26 gene. The FGF1 gene has been identified as being differentially expressed in animals phenotypically extreme for FA composition in muscle 16 . Being noteworthy, FGF1 is involved in preadipocyte differentiation and has been suggested to be acting in the PPARG system, however, the mechanisms remain unclear 33 . In accordance with this hypothesis, it is described that FGF1 can be trans-located to the nucleus, exerting a growth regulatory activity 34 . Hence, it would be important to study the FGF1-PPARA relationship. In the SSC6, PPARA trans-eQTL mapped the Palmitoyl-Protein Thioesterase 1 (PPT1) gene, Metallophosphoesterase 1 (MPPE1), Inositol(Myo)-1(or 4)-Monophosphatase 2 (IMPA2), and Cell Death-Inducing DFFA-Like Effector A (CIDEA). Interestingly, CIDEA-null mice showed a decreased PPARA, PPARG and SREBF1 gene expression and a decreased de novo FA synthesis in the liver 35 . Furthermore, the CIDEA-PPARA interaction identified in the eQTL analysis was also captured by Genomatix literature-based analysis (Supplementary Fig. S1). In addition, two melanocortin receptor genes (MC2R and MC5R) mapped in the PPARA SSC6 trans-eQTL. In adipocyte cells of MC2R knockdown mice alterations in fatty acid composition were observed: a reduction in the C16:1/C16:0 and C18:1/C18:0 ratios and an increase in the arachidonic acid content 36 . On the other hand, MC5R has been associated with obesity phenotypes such as body mass index and fat mass in humans 37 Table S3). Interestingly, the PRKAA2 gene identified within the SSC6 trans-eQTL for PIK3R1 is a catalytic subunit of the AMP-Activated Protein Kinase (AMPK) which is in charge of regulating lipid synthesis by phosphorylating lipid metabolic enzymes such as ACACA, ACACB, ACC, GYS1, HMGCR, HSL and LIPE to inactivate them. Therefore, it acts regulating key enzymes of fatty acid uptake, esterification, lipolysis and oxidation 38 . In the same way, PRKAA2 knockdown affects Akt activation. Therefore, PIK3R1 and PRKAA2 are both associated with the PI3K-Akt signalling pathway. Supporting these results, recently published studies suggest that AMPK activates Akt via regulating PI3K 39 . These results highlight PRKAA2 as a strong candidate gene to explain the variation in the mRNA levels of PIK3R1.  Table S3) were mapped. Interestingly, a polymorphism in the porcine MTTP gene has been associated with the lipid transfer activity of the MTTP protein and with the fatty acid composition of fat 40 . ADH5 was the most expressed gene in pig muscle; meanwhile, the other ADH members (ADH7 and ADH4), jointly with MTTP and HPGDS genes, showed very low expression levels in muscle 16 , which is in concordance with what was reported in human muscle 23 . Two potential regulators of PLA2G12A were identified within the trans-eQTL at the 79.9 Mb position on SSC6: Mitochondrial Trans-2-Enoyl-CoA Reductase (MECR) and Sestrin 2 (SESN2) ( Supplementary Table S3). Finally, only one lipid-related gene was identified within the second trans-eQTL at 9.7 Mb, the WW Domain Containing Oxidoreductase (WWOX) gene, which has recently been described to play a role in cholesterol homeostasis and triglyceride biosynthesis 41 .
On SSC9, an eQTL at around 117 Mb affected the expression of genes PLA2G12A and the Hypoxia Inducible Factor 1, Alpha Subunit Inhibitor (HIF1AN) ( Table 1). Their mRNA expression was highly correlated (r = 0.60; p-value = 1.98 × 10 −12 ), suggesting a common element regulating their transcriptional levels. The three most significant SNPs (H3GA0028012, ASGA0044215, ALGA0117195; p-value = 4.48 × 10 −6 ) within this eQTL were in complete linkage disequilibrium (R 2 = 1) ( Table 1). In this region PIK3CG was mapped, which encodes a Phosphatidylinositol-4,5-Bisphosphate 3-Kinase which is known to participate in different functions, including the insulin-signalling pathway and lipid metabolism 42 , and the DLD gene which encodes a Dihydrolipoamide Dehydrogenase involved in acetyl-coA biosynthesis (Supplementary Table S3). The PIK3CG gene, unlike other PI3K family members, is activated by interaction with G-protein-coupled receptors and by silencing this gene the PI3K-Akt signalling pathway is inhibited 43 . Accordingly, studies in cell lines have suggested that HIF1AN gene expression is repressed by a mechanism involving PI3K signalling 44 . Altogether, these results suggest that PRKAA2 may regulate the Class I PI3K Regulatory Subunit 1 (PIK3R1), which is able to form a heterodimer with the PIK3CG catalytic subunit, activate the Akt pathway and inhibit HIF1AN gene expression.
For IGF2 gene expression, apart from the SSC2 cis-eQTL, a trans eQTL located approximately at 162 Mb in the same chromosome was identified, and the Sirtuin 3 (SIRT3) gene mapped in this region (Supplementary Table  S3). It has been described that SIRT3 knockout mice exhibited decreased oxygen consumption and increased oxidative stress in skeletal muscle 45 . Moreover, SIRT3 knockout mice showed a down-regulation of Akt phosphorylation 45 (Supplementary Fig. S2). Members of the SIRT gene family have been described as being relevant genes controlling lipolysis and promoting fat mobilization in white adipose tissue 46 , and the SIRT1 member has been identified as one of the most central genes in a liver co-association network for intramuscular FA composition in IBMAP animal material 18 . Recently, the SIRT3 member, which can be activated by the AMPK protein, has been suggested as playing a major role in obesity-related diseases 47 .
Finally, within the trans-eQTLs identified for CROT, FABP3, MGLL, and NCOA1 gene expression, no strong candidate gene exerting a known lipid metabolism function could be detected.

Functional network analysis of genes mapping in eQTLs.
For trans-eQTLs, all genes located within a ± 1 Mb interval were selected for gene annotation. Conversely, for cis-eQTLs, only the candidate gene studied was considered (ACSM5, IGF2, and MGLL) for further analyses. In the 18 eQTLs, a total of 292 protein-coding genes, 13 miRNA, one miscRNA, six pseudogenes, one rRNA, eleven snoRNAs and four snRNAs were annotated (Supplementary Table S4). From the 292 protein-coding genes with Ensembl Gene ID, 256 had at least one human orthologous gene and were submitted to IPA to perform a functional categorization (Supplementary Table S4).
Focusing on the first network comprising energy production, small molecule biochemistry and drug metabolism functions, we identified the Serine/Threonine Kinase Effector (Akt) complex as central in the network (Supplementary Fig. S2). Remarkably, the Akt complex, which is involved in glucose transport and lipogenesis, was also identified in the muscle transcriptome study of 12 BC1_LD animals which were extreme for intramuscular FA composition 16 . In agreement with these results, several genes (PIK3CG, PPAP2B, PRKAA2, PTPN2, and SIRT3), identified as potential regulators within the eQTLs of the 45 lipid-related genes, are strongly related to the Akt pathway.

Identification of master regulators. A total of 298 genes, including (1) the 253 genes annotated in the
trans-eQTL intervals and, (2) the 45 genes studied of the present study were analysed with iRegulon Cytoscape plugin 48 . We observed that the EP300 gene was the most enriched transcription factor motif (enrichment score threshold for the motif, NES = 4.737). What was noteworthy was the EP300 gene was identified as a key regulator of FA composition and IMF traits in the same material using a gene co-association network 18 . On the other hand, five of the 45 genes studied were identified as strong regulators of several genes analysed with iRegulon: NR1H3 (NES = 5.348; 10 target genes), MLXIPL (NES = 4.988; 29 target genes), PPARA (NES = 4.785; 32 target genes), NFKB1 (NES = 3.855; 6 target genes) and PPARG (NES = 3.536; 10 target genes). The PPARA motif was present in the highest number of target genes (32 out of 45 genes). Interestingly, the NR3C1 gene identified within the PPARA trans-eQTL on SSC2 was also identified with iRegulon among the regulators for the 45 genes analysed and MatInspector (Genomatix software) predicted a binding site for this transcription factor in the promoter of PPARA. The NR3C1 gene may play a negative role in adipogenesis, regulating the lipolytic and anti-lipogenic gene expression. In human studies, polymorphisms in the NR3C1 gene have been suggested as contributing to obesity 49 . Afterwards, by using the Genomatix interface, it was observed that NR3C1 interacts with NFKB, CREB, NCOA1 and NCOA2 genes, and may be affected by the Corticotropin-Releasing Hormone (CRH) and Insulin (INS) (Supplementary Fig. S3). Accordingly, CRH was identified in both network analysis for fatness and growth traits 17 and the insulin-signalling pathway in a RNA-Seq study comparing animals extreme for their intramuscular FA composition 16 , both studies performed with IBMAP animal material. Thus, we hypothesize that NR3C1 may be a master regulator of lipid metabolism through the regulation of PPARA gene expression. However, further studies are necessary to corroborate this hypothesis.
QTL and eQTL co-localization. Transcriptomic data can be used to identify candidate genes underlying QTLs through the co-localization with eQTLs. Thus, our group looked for the overlapping of the 241 SNPs with the QTLs described in PigQTLDB 50 . A total of 234 SNPs (97%) co-localized in 132 QTLs for fatness traits and 157 SNPs (65%) within 10 different QTLs for FA composition (Supplementary Table S5), confirming a high co-localization of eQTLs and fat-related QTLs. The large number of QTLs described for these traits, covering a large extension of the porcine genome, provides evidence for a complex genetic pleiotropic regulation basis of these traits.
Notably, five genes identified within the eQTLs (i.e. ARHGAP6, IGF2, MC2R, MGLL, NR3C1) overlapped with QTLs described in the IBMAP population. For instance, the IGF2 gene, for which a cis-eQTL was detected, was identified within a confidence interval of one of the epistatic regions affecting muscle fibre traits 51 . The MGLL gene, showing a cis-acting eQTL, was located close to a GWAS interval affecting FA composition traits: palmitic (C16:0) and oleic (C18:1(n-9)) FAs, the polyunsaturated/saturated ratio (PUFA/SFA ratio), SFA, and unsaturated index (UI) 12 . Furthermore, MGLL maps within a QTL interval affecting growth 11 . The MC2R gene, identified in a PPARA trans-eQTL, maps within an IBMAP QTL region for IMF content and backfat thickness 52 . Moreover, ARHGAP6 and NR3C1 genes also identified in a trans-eQTL for PPARA gene expression are located within a QTL for growth traits in the IBMAP population 11 .

Conclusions
In the present study, we provide a list of potential candidate genes and genetic variants that may be regulating the transcriptional level of relevant lipid-related metabolism genes. Combined assessment of the results obtained from eGWAS and GWAS may provide complementary information of genes and variants determining the IMF and FA composition traits. Therefore, the present results identify potentially key genes and variants affecting pork-meat quality. However, more efforts should be made to validate our results, for instance, the involvement of the NR3C1 gene as a major regulator in muscle FA metabolism.
Scientific RepoRts | 6:31803 | DOI: 10.1038/srep31803 Methods Animal samples and phenotypes. The IBMAP population was obtained by crossing three Guadyerbas Iberian boars with 31 Landrace sows 6 . In the present study, 114 animals were used (65 females and 49 males) belonging to the BC1_LD generation of the IBMAP population obtained by crossing five F1 boars with 26 Landrace sows. Animals were fed ad libitum with a cereal-based commercial diet and slaughtered at 179.8 days ± 2.6 days. Animal care and procedures were performed following national and institutional guidelines for the Good Experimental Practices and approved by the Ethical Committee of the Institution (IRTA-Institut de Recerca i Tecnologia Agroalimentàries). Samples of the Longissimus dorsi muscle were collected, snap-frozen in liquid nitrogen and stored at − 80 °C until further RNA isolation. Genomic DNA was extracted from blood samples according to the phenol-chloroform method 53 .

Selection of lipid-related metabolism genes in muscle.
In previous studies of our group, strong candidate genes affecting IMFA content and FA composition of the Longissimus dorsi muscle of the IBMAP BC1_LD, 25% Iberian and 75% Landrace, were identified by using GWAS, RNA-Seq and co-association network approaches 12,[16][17][18] . In the present study, a list of 45 genes functionally related with lipid metabolism was selected, prioritizing candidate genes for FA composition (Supplementary Table S6).
We included: (1) 54 (freely available at http://www.bioguo.org/AnimalTFDB/), enzymes (ACSS1, ACSS2, CPT1B, CROT, DGAT1, DGAT2, PDHX), as well as the IGF2 gene, which has been described as the causal factor of the imprinted QTL for muscle growth and fat deposition in a Meishan × Large White intercross 27 (Supplementary  Table S6). Moreover, polymorphisms in the IGF2 gene have been significantly associated with muscle FA composition in swine 55 . Genotyping. A total of 197 animals from the BC1_LD (160 backcrossed individuals and their respective parents) were genotyped with the Porcine SNP60 Beadchip (Illumina), following the Infinium HD Assay Ultra protocol (Illumina) 56 . Raw data were visualized with GenomeStudio software (Illumina) and trimmed for high genotyping quality (call rate > 0.99). Markers with minor allele frequency (MAF) > 5% and animals with missing genotypes < 5% were retained. After the quality control filter, a subset of 40,586 SNPs and 197 animals remained.
Gene-expression profiling. Total RNA was isolated from the Longissimus dorsi muscle of 114 samples with RiboPure kit (Ambion; Austin, TX, USA). Total RNA was quantified in a NanoDrop ND-1000 spectrophotometer (NanoDrop Products; Wilmington, DE, USA). The RNA was converted to cDNA using the High-Capacity cDNA Reverse Transcription (Applied Biosystems). The cDNA samples were loaded into a Dynamic Array 48.48 chip in a BioMark system (Fluidigm; San Francisco, CA, USA) through a NanoFlex controller following a protocol previously described 18 . For this experiment, the expressed levels of 48 genes were analysed, and a total of 45 target genes were normalized using the two most stable housekeeping genes (ACTB and TBP). Primers used for the analyses are detailed in Supplementary Table S8. Data were collected using the Fluidigm Real-Time PCR analysis software 3.0.2 (Fluidigm) and analysed using the DAG expression software 1.0.4.11 57 , applying the relative standard curve method (See Applied Biosystems user bulletin #2). Analyses were performed using the normalized gene-expression levels of each sample and assay. The animals showing abnormal gene expression levels (outliers) were removed and data obtained were normalized if necessary by performing log 2 transformation of the Normalized Quantity (NQ) value. The sex effect was also tested by using a linear model with R programme 58 .
Gene-expression association analysis. An eGWAS was also performed using the genotypes of BC1_LD animals and the expression values from muscle. A mixed model was employed implemented on Qxpak 5.0 59 : ijkl i j k l k i jkl in which y ijkl was the k th individual record, sex (two levels) and batch (five levels) were fixed effects, λ k was a − 1, 0, + 1 indicator variable depending on the k th individual genotype for the l th SNP, a l represents the additive effect associated with the l th SNP, u k is the infinitesimal genetic effect treated as random and distributed as N(0, Aσ u 2 ), where A is the numerator of the pedigree-based relationship matrix and e ijkl the residual. The same mixed model was applied to perform the association analyses with the ACSM5 (rs331702081) and IGF2 (IGF2-intron3-G3072A) polymorphisms and the ACSM5 and IGF2 mRNA expressions, respectively. To correct for multiple testing, the false discovery rate (FDR) was calculated with the q-value library of the R package setting the threshold at q-value ≤ 0.05 58,60 . Given that cis-acting eQTLs tend to have larger magnitudes of effect on gene expression than do trans-acting eQTLs 61 , and several SNPs are affected by strong LD, a more restrictive threshold to avoid false positives was applied at p-value ≤ 0.001 for those eGWAS showing a cis-eQTL (see Q-Q plot in Supplementary Fig. S4).
The SNPs identified were classified as cis when they were located within 1 Mb from the gene analysed and as trans when they were located elsewhere in the genome. The number of significant SNPs belonging to the same interval was considered among associated SNPs less than 10 Mb apart.
Gene annotation and functional analysis. The significantly associated SNPs were mapped in the Sscrofa10.2 assembly and were annotated with the Ensembl Genes 80 Database using VEP software 62 . The genomic eQTL intervals considering ± 1 Mb around the candidate chromosomal regions were annotated using BioMart software 63 .
The Core Analysis function included in the Ingenuity Pathway Analysis software (IPA 4.0, Ingenuity Systems Inc., http://www.ingenuity.com; Accessed 8 July 2015) and the Genomatix Pathway System (GePS) (Release 2.8.0) in the Genomatix software suite (Genomatix Software GmbH; Munich, Germany; available at: http://www. genomatix.de; accessed 8 October 2015) were used to perform the functional analysis of genes mapped in the 18 eQTLs regions. Specifically, the IPA software was used for data interpretation in the context of biological processes, pathways, networks and upstream regulators. All information generated in this software is derived from the Ingenuity Pathway Knowledge Base (IPKB), which is based on functions and interactions of genes published in the literature. Genomatix was used to retrieve additional information of gene functions, interactions and upstream regulators based on the literature. Furthermore, information from Mouse Genome Database 64 and Genecards 65 was used to identify gene functions affecting the phenotypes analysed. For the lipid-related genes, the Gene Expression Atlas 23 was used to determine whether they were expressed in muscle or not.
Finally, an in-silico identification of transcription-factor binding sites in the promoter region of the 256 annotated genes was performed. For this analysis, iRegulon 48 was used, which relies on the analysis of the motif enrichment for a transcription factor in the gene set using databases of nearly 10,000 TF motifs and 1,000 ChIP-seq data sets or "tracks".

Correlation of gene expression and phenotypes.
Correlations were performed among gene expression of the 45 genes to explore the relationship between genes. Furthermore, pairwise correlations among gene expression and FA composition percentages in muscle 12 were carried out to explore the relationships between gene expression and phenotypes. All values were normalized applying the log 2 of raw data if necessary. Afterwards, gene expression was corrected by sex (two levels) and batch (five levels) effects, whereas FA composition percentages were corrected by sex, batch and carcass weight. The remaining residuals of the phenotypes and gene-expression values corrected for the corresponding effects were used to obtain the pairwise correlations. The hierarchical clustering option of PermutMatrix software 66 was used to visualize the results of both traits and genes.