Expression analysis of candidate genes for fatty acid composition in adipose tissue and identification of regulatory regions

The aim of this work was to study the genetic basis of the backfat expression of lipid-related genes associated with meat quality traits in pigs. We performed a genome-wide association study with the backfat gene expression measured in 44 genes by qPCR and the PorcineSNP60 BeadChip genotypes in 115 Iberian x Landrace backcross animals. A total of 193 expression-associated SNPs located in 19 chromosomal regions were associated with expression levels of ACSM5, ELOVL6, FABP4, FADS2, and SLC27A4 genes. Three expression quantitative trail loci (eQTLs) corresponding to ACSM5, FABP4, and FADS2 were classified as cis-acting eQTLs, whereas the remaining 16 eQTLs have trans-regulatory effects. Remarkably, a SNP in the ACSM5 promoter region and a SNP in the 3′UTR region of FABP4 were the most associated polymorphisms with the ACSM5 and FABP4 expression levels, respectively. Moreover, relevant lipid-related genes mapped in the trans-eQTLs regions associated with the ACSM5, FABP4, FADS2, and SLC27A4 genes. Interestingly, a trans-eQTL hotspot on SSC13 regulating the gene expression of ELOVL6, ELOLV5, and SCD, three important genes implicated in the elongation and desaturation of fatty acids, was identified. These findings provide new data to further understand the functional regulatory mechanisms implicated in the variation of fatty acid composition in pigs.

Pork meat is an appreciated, all-purpose lean meat, and represents one of the main sources of animal meat for humans 1 . Fat and fatty acids (FAs), are fundamental to various aspects of meat quality and play a crucial role to meat nutritional value, both in adipose tissue (backfat, BF) and muscle (intramuscular fat 2 ). It is demonstrated that FA composition is dependent on physiological status, nutrition conditions 2,3 , and genetic factors 4 .
In the last few years, genome-wide association studies (GWAS) have been performed in attempts to uncover the genetic basis of FA composition traits. Studies of our group and others have benefited from this approach and have identified genomic regions significantly associated with intramuscular fatty acid (IMFA) composition by using different experimental and commercial populations [5][6][7][8][9] . However, in most cases, the GWAS approach only allows the detection of genetic variants that explain a modest proportion of the total heritability of the analyzed traits 10 . In addition, the way from the genomic statistical association to the identification of true causal genetic variants is plagued of difficulties. Hence, it has become evident the necessity to integrate new approaches to better understand the biological significance of GWAS findings. Recently, the association between genetic variants and gene expression levels has been described and used to identify expression quantitative trait loci (eQTLs). The eQTLs regulating the transcript abundance of the mRNAs can be identified systematically using high-throughput technologies and have recently been proposed as a good strategy to deepen the study of the genetic architecture of complex traits 11,12 . Liver, skeletal muscle and adipose tissue are three of the most important tissues involved in FA metabolism 13 . Adipose tissue, one of the main energy reserves in animals, is composed of adipocytes embedded in a matrix of connective tissue with a highly developed vascular system. The adipocytes are dynamic cells that play a relevant role in energy balance and overall body homeostasis. Their main metabolic functions are to accumulate lipids, by synthesis of triacylglycerols, and lipid mobilization, through hydrolysis of triacylglycerols 14 .
Previous studies in our group have identified genomic regions, candidate regulators and regulatory polymorphisms in the liver and muscle tissues of individuals of an Iberian x Landrace backcross population (BC1_LD 15,16 ). Furthermore, a RNA-Seq transcriptome study of the adipose tissue of two groups of pigs with phenotypically extreme IMFA composition in the BC1_LD animals identified metabolic pathways differentially modulated between groups controlling lipid and FA metabolism 17 . Taking into account the relevant role of adipose tissue in the regulation of lipid metabolism, the goals of the present work were (1) to study the expression of candidate genes related with lipid metabolism in the adipose tissue of BC1_LD animals and (2) to analyze the chromosomal regions significantly associated with the gene expression levels to characterize the regulatory mechanisms influencing gene expression phenotypes. With these results we aim to increase our knowledge of the regulation of key genes determining FA composition in pigs in adipose tissue.

Results and Discussion
Selection of genes related with lipid metabolism in the adipose tissue. Using the information generated in previous studies of our group, strong candidate genes affecting FA composition of BF and intramuscular fat in the BC1_LD generation were identified by GWAS, RNA-Seq and co-association network approaches 5,[17][18][19][20][21] (Supplementary Table S1). A total of 44 candidate genes related with lipid metabolism were selected to study their expression pattern in BF. Fourteen of them (ARNT, CYP2U1, EGF, ELOVL6, FABP4, FABP5, FADS1, FADS2, FADS3, NFKB1, PLA2G12A, PLCB2, PLPP1, and USF1) are functional and positional candidate genes related with lipid metabolism which were identified in GWAS analyses for BF and IMFA composition in the BC1_LD animals 5,21 . We also included two candidate genes differentially-expressed (ELOVL6 and SCD) in a RNA-Seq analysis of the adipose tissue of two phenotypically extreme groups of animals for IMFA composition in the BC1_LD cross 17 . Lipid metabolism genes identified in gene co-association networks for FA composition (ACSM5, ANK2, ARNT, FABP4, FABP5, MGLL, and PPARG) 19 were also included; two of which, ARNT and PPARG, were also identified in gene co-association networks for fatness and growth traits 20 . In addition, in order to complete the set of genes, we included genes which have been described in the literature to play different roles in lipid metabolism such as transporters (RBP4, SCAP, SLC27A1, and SLC27A4), enzymes (AGPAT2, CPT1A, CROT, DGAT1, DGAT2, ELOVL5, LIPC, LPIN1, PEX2, and PNPLA2) and transcriptional factors, cofactors or nuclear receptors (CD36, ESRRA, MLXIPL, NR1H3, POU2F1, PPARA, PPARD, PPARGC1A, RXRG, and SREBF1). Finally, we added the ADIPOQ gene, which encodes an important adipokine of white fat tissue, exerting multiple biological processes on carbohydrate and lipid metabolism (review in Shehzad et al. 22 ). Although this gene is not annotated in the current Sscrofa10.2 assembly, it has been mapped to the SSC13q36-41 interval 23 .
Four reference genes (ACTB, B2M, HPRT, and TBP) were also selected as reference controls.
Expression genome-wide association analysis. In the present study, the adipose tissue expression of 48 genes (44 target and four reference genes) was measured by quantitative real-time PCR (qPCR) in 115 BC1_LD animals. For the PPARGC1A gene, a poor PCR efficiency was obtained and it was discarded for further analysis. Expression GWAS (eGWAS) were performed with the gene expression values of the 43 remaining target genes and the genotypes of 40,460 SNPs of the PorcineSNP60 BeadChip (Illumina) that passed the quality control (see Methods section). At whole genome level (q-value ≤ 0.05), significant association signals in five of the analyzed genes were detected (Table 1, Fig. 1, Supplementary Fig. S1): ACSM5, ELOVL6, FABP4, FADS2, and SLC27A4. The ACSM5, FABP4, FADS2, and SLC27A4 genes presented more than one associated eQTL (Table 1). Three out of 19 eQTLs were identified as cis-acting for the ACSM5, FABP4, and FADS2 gene expression ( Fig. 1), suggesting the presence of proximal polymorphisms regulating the expression of these genes. These results showed a difference in the ratio of cis-and trans-eQTLs. In general, studies performed in animals have identified regulatory trans-eQTLs in a higher ratio than those performed in humans 12,24 . Accordingly, previously eGWAS analyses performed in muscle 15 and liver 16 of BC1_LD animals identified a prevalence of porcine trans-eQTLs compared wih cis-eQTLs.
From the associated SNPs, and according to the Variant Effect Predictor (VEP) of Ensembl (Sscrofa10.2 annotation release 84), 49.2% (95 SNPs) were located in intergenic regions. The remaining 50.8% (98 SNPs) mapped within a total of 68 genes: 66 in introns, 11 in 5′ upstream and 16 in 3′ downstream gene regions, three in the 3′UTR region, and two in the coding region of a gene determining synonymous mutations (Supplementary  Table S2). A total of 86 SNPs (44.6%) were located inside the cis-acting eQTLs, whereas 107 SNPs (55.4%) were in trans-eQTLs (Supplementary Table S2).
In the following sections, the SNPs and genes mapped in the eQTL regions associated with the expression phenotypes of ACSM5, ELOVL6, FABP4, FADS2, and SLC27A4 are discussed in detail. The selection of genes in the eQTLs was performed by: (1) extracting the annotated genes from the eQTL intervals including an additional Mb at each end, (2) selecting genes with a lipid-related function and (3) expressed in adipose tissue. Table 2 summarizes all the relevant lipid-related genes mapped in the trans-eQTL regions.
The ACSM5 gene. The most significant cis-SNPs (p-value = 1.11 × 10 −17 ) for the Acyl-CoA Synthetase Medium-Chain Family Member 5 (ACSM5) gene expression were located between 23 disequilibrium (LD) (r 2 = 0.42-0.99). Recently, the polymorphism ACSM5:g.26260422 G>A (rs331702081), located at the proximal promoter region of the ACSM5 gene, has been described as the most associated with the expression of ACSM5 in Longissimus dorsi muscle of BC1_LD animals 15 . Hence, the ACSM5:g.26260422 G>A SNP was genotyped in the 115 BC1_LD animals and included in the eGWAS study, showing that is one of the SNPs located in the 23. 44-27.94 Mb block with the lowest p-value (p-value = 1.11 × 10 −17 ; Fig. 1A). Consequently, this SNP may be a strong candidate to explain the mRNA variation of the ACSM5 gene in BF and muscle. Even so, the correlation value between the ACSM5 gene expression in BF and muscle is of r = 0.60 (p-value = 3.0 × 10 −12 ), suggesting that other factors than the SSC3 cis-eQTL may be differentially regulating the expression of ACSM5 in both tissues.
The effect of ACSM5:g.26260422 G>A SNP on the binding of transcription factors (TFs) was determined by the LASAGNA-Search version 2.0 software 25 . The aryl hydrocarbon receptor nuclear translocator (ARNT) gene (TF ID = M00236 and M00539) and the signal transducer and activator of transcription 6 (STAT6) gene (TF ID = M00500) were identified to bind only when the A allele is present. ARNT gene plays an important role in the regulation of hepatic lipogenesis and gluconeogenesis 26,27 , and was identified as one of the most central genes in a liver co-expression network analysis of IMFA composition in pigs 19 . STAT6 gene has been described to interact with the peroxisome proliferator activated receptor gamma gene (PPARG), and the cooperative binding of the two genes led to an increase response of PPARG 28 . The importance of PPARG gene lies on its regulation of adipocyte differentiation and glucose homeostasis, and it was identified as a major regulator for growth and fatness related traits in a co-association network in muscle of BC1_LD individuals 20 . Different transcription binding sites for PPARG in the ACSM5 proximal promoter region were also predicted by LASAGNA (Supplementary Table S3).
Five chromosomal regions (Table 1) were also associated in trans with the ACSM5 gene expression. Likewise, two of these regions on SSC3 (100.91-101.46 Mb) and SSC10 (0.05-0.20 Mb) have been recently associated with the mRNA expression of the ACSM5 gene in the Longissimus dorsi muscle of BC1_LD animals (SSC3: 100.35 Mb; SSC10: 0.17 Mb) 15 . These results, together with the cis-eQTL identification for the ACSM5 expression, further confirm the existence of common regulatory mechanisms implicated in the expression of ACSM5 in BF and muscle. No strong candidate genes exerting a known lipid metabolism function were detected in the trans-eQTL identified on SSC16 for ACSM5. Conversely, a list of lipid-related genes identified on SSC2, SSC3 and SSC10 ACSM5-trans associated regions is showed in Table 2.  The ELOVL6 gene. In the present study, previously described ELOVL6 polymorphisms (three in the promoter region, one in exon 4 and two in the 3′UTR region) 18,29 were genotyped in the 115 BC1_LD animals and incorporated in the eGWAS. The obtained results were consistent with our previous studies on SSC8 18 , showing a chromosome level significant region associated with ELOVL6 mRNA expression in BF ( Supplementary  Fig. S2). The most significant peak was for two SNPs with the same p-value (ALGA0049135 and ALGA0049139; p-value = 4.60 × 10 −5 ; â = 0.563), located in an intron of ANK2, at 2.4 Mb to the ELOVL6 gene. In contrast the ELOVL6:c.-533C>T and ELOVL6:c.-394G>A polymorphisms, which were in full LD and also showed high association (p-value = 6.42 × 10 −04 ; â = 0.460), did not achieve statistical significance after multiple testing corrections (q-value ≤ 0.05; Supplementary Fig. S2). This observed discrepancy could be due to the different animals used between the two studies (94 samples overlapping between the two studies) and to the subtle differences observed between the two methods used for mRNA quantification (Pearson correlation coefficient, r = 0.91).
Interestingly, ANK2 is one of the most central genes in an adipose co-expression network related with IMFA composition in the BC1_LD 19 . In a recent study, knockin mice expressing ANK2 genes with human nonsynonymous mutations showed altered glucose homeostasis contributing to increased adiposity. This phenotype was caused by the reduction in ANK2 protein levels which produce an elevation of cell surface GLUT4 and increased glucose uptake in skeletal muscle and fat 30 . Thus, we cannot discard an association of this gene with the mRNA expression levels of ELOVL6. Further analyses are necessary to corroborate this hypothesis.
Regarding the trans-associated regions with the BF mRNA ELOVL6 expression, Corominas et al. 18 did not identify any region at whole genome level (q-value ≤ 0.05), in contrast with the trans-eQTL on SSC13 identified in the present study. No genes related to lipid metabolism function were identified in the SSC13 trans-eQTL. However, we have to highlight that the same region on SSC13 was associated at chromosome level with the expression of the ELOVL5 and SCD genes. ELOVLs and SCD genes are implicated in the elongation and desaturation of FAs in the endoplasmic reticulum membrane. These metabolic functions are essential to the maintenance of lipid homeostasis (reviewed in Guillou et al. 31 ). Moreover, it has been reported that the expression of these genes are primarily regulated at transcriptional level (reviewed in Guillou et al. 31 ). In this study, highly significant Pearson's correlations (p-value = 2.22 × 10 −16 ) among the mRNA expression of these genes were obtained (r ELOVL5-ELOVL6 = 0.90; r ELOVL5-SCD = 0.84; r ELOVL6-SCD = 0.86), supporting the involvement of common elements regulating their mRNA expression. In addition, we found high correlations with the genes coding for the TFs SREBF1 and PPARs (r ELOVL5- which contribute to the regulation of these genes 17,31,32 . High correlation values were also found between the mRNA expression of DGAT1 and DAGT2 genes with ELOVLs and SCD genes ( . This result agrees with the interrelated function of these genes which are implicated in lipogenesis (ELOVLs and SCD) and triglyceride (TG) synthesis (DGAT1 and DGAT2), the main function of adipose tissue 14 .
ELOVL6 and SCD genes were over-expressed in the adipose tissue transcriptome (RNA-Seq) of BC1_LD animals with higher content of intramuscular monounsaturated fatty acids (MUFA) and saturated fatty acids (SFA), when compared with animals having more polyunsaturated fatty acids (PUFA) 17 supporting the relevance of these genes in the determination of FA composition in the BC1_LD.
The FABP4 gene. The Fatty Acid Binding Protein 4 (FABP4) gene is not annotated in the current Sscrofa10.2 assembly. Nevertheless, radiation hybrid (RH) and linkage maps located FABP4 gene close to FABP5 on SSC4 33,34 , in agreement with the human comparative map. For this reason, we used the position of the annotated FABP5 gene (SSC4: 60.31 Mb) in order to define cis/trans eQTLs for FABP4 gene. Recently, an indel located in the intron 1 of FABP4 (FABP4:g.2634_2635insC) was identified as the most associated polymorphism with the FABP4 mRNA levels in BF 35 . In the present study, this indel, which was genotyped in the 115 BC1_LD animals and included in the eGWAS study, was one of the most significant associated polymorphisms (p-value = 1. 27 × 10 −6 ; â = −0.219) (Fig. 1B), but the most significant cis-SNP (ALGA0025337; p-value = 8.67 × 10 −7 ; â = 0.328) was located in an intergenic region (Supplementary Table S2).
Taking into account the highly polymorphic nature of the porcine FABP4 gene 36 and a possible regulatory role of miRNAs in the FABP4 gene expression 37 , the 3′UTR of the FABP4 gene was amplified and sequenced using 10 animals with extreme values of FABP4 mRNA expression in BF. A polymorphism in the 3′UTR (FABP4:g.6723 A>G) was identified and genotyped in the 115 BC1_LD animals. The association analysis with the FABP4 mRNA expression values revealed that the SNP in the 3′UTR region has the lowest p-value (p-value = 3.07 × 10 −8 ; â = −0.260) (Fig. 1B).
Interestingly, FABP4:g.6723 A>G SNP was inside a putative miRNA binding site for the putative human miRNA hsa-miR-3182 and it is predicted to bind only when FABP4:g.6723 G allele is present. However, we did not find the homologous Sus scrofa miR-3182 in the current assembly (Sscrofa10.2) of the pig genome. Although the 3′UTR SNP (FABP4:g.6723 A>G) is a clear candidate to explain the differences of FABP4 mRNA levels among animals, we cannot discard also a role of the indel (FABP4:g.2634_2635insC) in the FABP4 gene regulation. Indeed, the indel polymorphism was predicted to be located in a target binding site for PPARG and NR4A2 TFs 35 . In our analysis the correlation between the mRNA expression of PPARG and FABP4 was r = 0.51 (p-value = 9.27 × 10 −9 ). Further functional analyses are needed to test the role of these polymorphisms in the The eGWAS performed for the FABP4 gene also revealed four trans-eQTLs on SSC2, SSC3, SSC4, and SSC9 (Table 1); where several genes involved in lipid metabolism were mapped ( Table 2).
It is noteworthy to highlight some genes identified in the trans-eQTLs, such as the Berardinelli-Seip Congenital Lipodystrophy 2 (Seipin) (BSCL2) gene located on SSC2. Although the function of this gene is still being investigated, its relationship with adipogenesis, with the genesis of lipid droplets and the regulation of the metabolism of phospholipids and triacylglycerides has been established 38 . Liu et al. 39 performed an adipose-specific Seipin knockout mice with FABP-mediated BSCL2 deletion exhibiting decreased lipolysis in response to β-adrenergic receptors agonists in vivo. In this chromosomal region we also identified the Lectin Galactoside Binding Soluble 12 (LGALS12) gene, an intracellular galectin preferentially expressed in adipocytes which regulates lipolysis and whole-body energy metabolism 40 . Furthermore, FABP4 and LGALS12 genes were identified in the glycolysis/ gluconeogenesis pathway in an expression analysis performed in human adipose-derived stem cells 41 . Finally, the Phospholipase A2 Group XVI (PLA2G16) gene, which encodes the most abundant PLA 2 in the adipose tissue 42 , was also identified on SSC2. PLA2G16 has been described an important regulator of the rate of adipose tissue lipolysis via the production of eicosanoid mediators 42 . On SSC3, we identified the Cytochrome P450 Family 1 Subfamily B Member 1 (CYP1B1) gene. Liu et al. 43 have demonstrated that the reduction in body weight gain and white adipose tissue in CYP1B1 deficient mice exhibited coordinated decreases in FA synthesis (regulated by FABP4, among others) when compared with wild type ones.
The FADS2 gene. The most significant cis-SNP for FADS2 gene expression (ASGA0104083; p-value = 5.98 × 10 −8 ; â = 0.466) was located less than 0.23 Mb upstream of the FADS2 gene (Fig. 1C). FADS2 cluster together with FADS1 and FADS3 in a region on SSC2 44 . In liver, a high correlation between the mRNA expression levels of FADS1 and FADS2 has been observed (r FADS1-FADS2 = 0.92; p-value = 1.11 × 10 −17 ) and a common trans-associated chromosomal region was identified on SSCX 16 . Hence, these results pointed to common regulatory mechanisms controlling the mRNA expression of both desaturases in liver 16 . Conversely, in BF the correlation between the mRNA expression levels of FADS1 and FADS2 was lower (r FADS1-FADS2 = 0.63; p-value = 8.26 × 10 −14 ) and only the cis-acting element associated with the FADS2 mRNA expression was identified. Altogether these results suggest that common and specific elements are regulating the expression of desaturases in BF.
Conversely, the low correlation observed between FADS2 mRNA levels in BF and liver (r = 0.23; p-value = 1.91 × 10 −2 ) suggests a different mechanism of regulation in both tissues. In addition, the FADS2 gene was differentially-expressed between sexes depending on the tissue. Whereas in adipose tissue the FADS2 gene showed a higher expression in males, in liver the expression was higher in females 16 . Similar results have also been shown in rats, where the mRNA expression of FADS2 is higher in the liver of females than males 45 . To find polymorphisms which may modulate FADS2 expression a search of polymorphisms was performed by the analysis of whole-genome sequence data from seven founders of the IBMAP population with IGV 46 . Three polymorphisms were identified at positions g9118843C>T (rs331050552), g9118813G>A (rs321384923), and g9118721G>A (rs336076510) according to the Ensembl ENSSSCG00000013072 sequence. The most proximal 5′ mutation (g9118843C>T) was genotyped in the BC1_LD animals and included in the eGWAS analysis. However, no significant association was found between this mutation and the FADS2 mRNA expression. Further analyses are necessary to find new candidate cis-acting polymorphisms implicated in the regulation of FADS2 gene expression.
Three regions on SSC6, SSC8, and SSC10 were associated in trans with FADS2 gene expression where several lipid-related genes were identified (Table 2), highlighting the Lysine Demethylase 1A (KDM1A) gene, also known as LSD1, located on SSC6. This gene has been recently identified as a novel regulator of lipid metabolism as it is required for the expression of SREBF1 and for the efficient binding of SREBF1 to their target gene promoters 47 . Remarkably, SREBF1 regulates the expression of desaturases 48 , suggesting the KDM1A gene as a potential candidate regulator of FADS2 gene expression. Supporting this hypothesis, a transcription binding site for SREBF1 in the FADS2 proximal promoter region was predicted by LASAGNA (Supplementary Table S4).
The most associated SNP for FADS2 trans-eQTL on SSC8 was an intronic polymorphism (MARC0041089; p-value = 1.93 × 10 −5 ; â = 0.516) in the Insulin Like Growth Factor Binding Protein 7 (IGFBP7) gene (Supplementary Table S2). This gene encodes a member of the insulin-like growth factor (IGF)-binding protein and through their regulation of IGFs and insulin may influence the metabolism of adipocytes, having implications in the development of obesity and insulin resistance 49 .
Sex effect on gene expression in the adipose tissue. The control of FA homeostasis has been described to have a pronounced sexual dimorphism 51 . In agreement, sex differences in the control of gene expression have been identified in different tissues like muscle, liver, adipose tissue and brain in mice 52,53 . In the present work, a significant sex effect (p-value ≤ 0.05) on gene expression levels was detected in 20 out of the 43 genes analyzed (47%): ACSM5, AGPAT2, ANK2, ATGL, DGAT2, EGF, ELOVL5, ELOVL6, FADS2, MGLL, MLXIPL, NFKB1, PEX2, PLA2G12A, PLPP1, PPARA, PPARD, SCAP, SCD, and SREBF1 (Fig. 2). Overall, a higher number of genes were more expressed in females. Conversely, only three genes (FADS2, NFKB1, and PLPP1) showed a higher expression in males (Supplementary Table S6). Among them, some key regulators of lipid metabolism such as PPARA, SCD, and SREBF1 which may be determining the bias observed in the female over-expressed genes. In accordance, lipogenic genes such as SCD and SREBF1 were also more expressed in females in the Longissimus dorsi muscle 15 and liver 16 of pigs. However, we observed a different pattern of expression for PPARA, showing a higher expression level in females in liver and adipose tissue and a higher expression in males in muscle 15,16 . Finally, for ACSM5, AGPAT2, MLXIPL, EGF, ELOVL5, PEX2, and SCAP the same female-biased gene expression as in pig liver tissue 16 was observed. This agree with the results observed in mice, in which the analysis of genes differentially expressed between sexes in liver, adipose, muscle, and brain revealed common functionalities of steroid and lipid metabolism among the sexually dimorphic genes between liver and adipose tissues. Furthermore, the overlap of sexually dimorphic genes was higher between adipose and liver tissues (22.9%) than between adipose and muscle tissues (6.6%; 52 ). The highest rate of coincidence of female-biased gene expression among liver and adipose tissue may suggest common regulatory mechanisms in both tissues associated with sex differences.

Fatty acid composition and gene expression in the adipose tissue.
To explore the relationship between gene expression levels and FA composition percentages in BF and muscle, Pearson's correlation analyses were performed. In general moderate correlations were identified for BF ( Supplementary Fig. S3) and muscle ( Supplementary Fig. S4).
The cis-7 hexadecenoic acid (C16:1(n-9)) also showed negative correlation (r = −0.39; p-value = 2.05 × 10 −5 ) with SCD gene expression, which plays an important role in the regulation of the expression of genes involved in lipogenesis 56 . In accordance, Fernández et al. 57 have reported that SCD haplotypes, constituted by SCDg.2108 C>T, g.2228 T>C and G.2520 G>A, are associated with SCD gene expression in backfat and the fatty acid desaturation ratio in backfat and intramuscular fat in a Iberian x Duroc backcross.
Furthermore, the octadecenoic acid (C18:1(n-7)) was negatively correlated with CROT (r = −0.36; p-value = 1.02 × 10 −4 ) and USF1 (r = −0.35; p-value = 1.29 × 10 −4 ) gene expressions. The CROT gene plays a role in the pathway of FA beta-oxidation, and provides a crucial step in the transport of medium and long-chain acyl-coA out of the mammalian peroxisome to the mitochondria 58 . Moreover, USF1 gene is involved in the maintenance of high levels of FA synthase transcription under lipogenic conditions 59 .
For the FA composition percentages in muscle, PLPP1 mRNA expression was moderately correlated with the percentage of α-linolenic acid (C18:3(n-3)) (r = 0.39; p-value = 2.18 × 10 −5 ), the double bond index (DBI) (r = 0.31; p-value = 7.69 × 10 −4 ), the unsaturation index (UI) (r = 0.35; p-value = 1.33 × 10 −4 ), and the ratio PUFA/SFA (r = 0.31; p-value = 8.76 × 10 −4 ). Negative correlations were also identified with palmitic acid (C16:0) (r = −0.39; p-value = 1.56 × 10 −5 ). The great variability between the PLPP1 gene expression and many FA composition percentages may be explained by its multiple functions, PLPP1 converts phosphatidic acid to diacylglycerol, and function in de novo synthesis of glycerolipids as well as in receptor activated signal transduction mediated by phospholipase D 60 . Furthermore, this gene has been described as a catabolic enzyme for lysophosphatidic acid 61 , and has been identified in the enriched functional category from the lipid catabolic process as a highly expressed gene in a transcriptome analysis of Basque pigs 62 . Functional network analysis of genes mapping into eQTLs. To take a deep look at the regulatory mechanisms that are influencing the gene expression phenotypes, gene annotation of the 19 eQTLs chromosomal regions was performed. For trans-eQTLs all the genes located within ±1 Mb were selected for gene annotation. Conversely, for cis-eQTLs only the studied candidate gene was considered for further analyses. In the 19 eQTLs, a total of 474 protein-coding genes, 2 processed transcripts, 11 miRNAs, 3 miscRNAs, 10 pseudogenes, 6 rRNAs, 17 snoRNAs, and 21 snRNAs were annotated (Supplementary Table S5). From the 474 protein-coding genes with Ensembl Gene ID, 393 have at least one human orthologous gene (Supplementary Table S7) and were submitted to the Ingenuity Pathway Analysis software (IPA, Ingenuity Systems Inc.) to perform a functional categorization.
Several networks related to lipid metabolism were identified: (i) lipid metabolism, small molecule biochemistry, infectious diseases (score = 42; interval = 8); (ii) lipid metabolism, small molecule biochemistry, cellular assembly and organization (score = 32; interval = 8); (iii) carbohydrate metabolism, organ morphology, reproductive system development and function (score = 30; interval = 3), and (iv) amino acid metabolism, carbohydrate metabolism, molecular transport (score = 18; interval = 9) (Supplementary Table S8). Interestingly, two of the identified networks have the Akt complex as central (ID = 10 and 26; Supplementary Table S9). It is known that the PKB/Akt plays an important role in the insulin regulation of glucose transport 63 . In addition, it has been also identified as central in the main over-represented pathways in a muscle transcriptome study between individuals phenotypically extreme for IMFA composition 64 , in a muscle eQTL analysis of 45 lipid-related genes 15 , and in a liver eQTL analysis of 44 candidate genes related with lipid metabolism 16 . Remarkably, several of the potential regulators annotated in trans-eQTL for FABP4 and FADS2 (BSCL2, FADS1, FADS2, KDM1A, LGALS12, PLA2G16, PLB1, and PPP1CB) were identified inside the Akt pathway. These findings suggest this pathway as central in the genetic determination of FA composition traits in pigs.

Prediction analysis of transcription factor binding sites.
To identify master regulators of global gene expression in all the genes analyzed in the present study, an in silico identification of the TF binding sites in the promoter region of (1) the 415 genes with human orthologous ID annotated in the trans-eQTL intervals and, (2) the 43 candidate lipid-related genes was performed using the iRegulon Cytoscape plugin 65 . The PPARG gene was the most enriched TF motif (normalized enrichment score, NES = 5.134; 167 target genes, 30 out of the 44 analyzed lipid-related candidate genes; Supplementary Table S9). This is in accordance with the fact that this gene is an adipogenic TF considered an important regulator of lipid and carbohydrate metabolism 66 . In our analysis the mRNA expression of PPARG was highly correlated with the mRNA expression of genes implicated in lipid transport, lipogenesis and TG synthesis (Supplementary Table S10), suggesting an important role of this TF in lipid storage. Interestingly, this TF has been identified as a key regulator of FA composition in the same animal material using a co-association network analysis 19 and in a muscle eQTL analysis of 45 lipid-related genes 15 .

Conclusions
In the present study, potential candidate genes and genetic variants regulating the expression of lipid-related genes associated with meat quality traits in pigs have been identified by eGWAS analysis in the backfat adipose tissue of Iberian x Landrace backcrossed pigs. While some genetic variants identified in cis seem to be tissue specific (FADS2, FABP4), other may exert a similar regulatory control in different tissues (ACSM5). In addition to the well-known transcription factors SREBF1 and PPARs, that regulate the expression of lipogenic genes, a trans-eQTL on SSC13 identified for the ELOVL6 gene was associated with the expression of ELOVL5 and SCD genes, which are implicated in the elongation and desaturation of fatty acids, supporting the involvement of common elements regulating their mRNA expression in adipose tissue. These results increase our knowledge of the regulatory mechanisms implicated in adipose tissue lipid metabolism and its consequences in lipid-related traits. Animal Samples. The IBMAP resource population was used in this study. This population was established by crossing 3 Iberian (Guadyerbas line) boars with 31 Landrace sows 67 , and 5 F 1 males and 25 Landrace sows were crossed to obtain the BC1_LD generation (25% Iberian × 75% Landrace). The Iberian pig breed is characterized by high levels of IMF with higher MUFA and SFA content and an extremely fat body composition. In contrast, the Landrace breed has an efficient meat production, with low IMF and high PUFA content, and with a rapid growth and leaner carcass. Hence, the IBMAP experimental cross was generated to study the genetic basis of growth and meat quality traits differentiating the two breeds 68 . Here, we report results based on 115 BC1_LD pigs. All animals were maintained under intensive conditions and feeding was ad libitum with a cereal-based commercial diet. Backcross animals were slaughtered at an average age of 179.9 ± 8.9 days, and samples of BF tissue were collected, snap-frozen in liquid nitrogen and stored at −80 °C until analysis. Genomic DNA was extracted from diaphragm samples of all animals by the phenol-chloroform method 69  a minor allele frequency (MAF) less than 5% and SNPs with more than 5% of missing genotypes. The SNPs not mapped in the Sscrofa10.2 assembly were also excluded. After the quality control filter, a subset of 40,460 SNPs remained.

Ethics
In addition, ten polymorphisms were genotyped in the BC1_LD animals: two SNPs located in the proximal promoter region of the ACSM5 (g.26260422 G>A, rs331702081 15 ) and the FADS2 (rs331050552) genes, and one indel and one SNP located in the intron 1 (FABP4:g.2634_2635insC 72 ) and in the 3′UTR region of the FABP4 gene (FABP4:g.6723 A>G), respectively. For the ELOVL6 gene, three SNPs located in the promoter region (ELOVL6:c. -533C>T, ELOVL6:c.-480C>T, and ELOVL6:c.-394G>A 18,29 ), one in exon 4 (ELOVL6:c.416 C>T 18 ) and two in the 3′UTR region (ELOVL6:c.1408 A>G and ELOVL6:c.1922C>T 29 ) were genotyped. ACSM5 and FADS2 SNPs were genotyped using Taqman OpenArray ™ genotyping plates designed in a QuantStudio ™ 12 K flex Real-Time PCR System (ThermoFisher Scientific) using the SVGM platform (http://svgm.es/ca/Platforms/). The pyrosequencing protocol described by Mercadé et al. 72 and the High Resolution Melting methodology (HRM, ThermoFisher Scientific) were used for genotyping the indel and the SNP of FABP4 gene, respectively. The SNPs belonging to ELOVL6 gene were genotyped using the KASP SNP genotyping service from LGC (http://www. lgcgroup.com/genotyping/). Characterization of the 3′UTR of porcine FABP4 gene. The 3′UTR of the FABP4 gene was amplified and sequenced in 10 animals with extreme values of FABP4 mRNA expression in BF.
The 3′UTR variants of FABP4 gene were characterized by 3′-RACE PCR using UAP reverse primer and, FABP4-3NC-1-Fw and FABP4-3NC-2-Fw forward primers (Supplementary Table S11). The specific bands were excised from agarose gel and purified using NucleoSpin ® Gel and PCR Clean-up (Macherey-Nagel GmbH & Co. KG) and sequenced in forward and reverse directions.
All primers were designed using PRIMER3 software 73 based on the Y16039 sequence 33   Gene expression profiling. Total RNA was isolated from the BF of the 115 BC1_LD samples with RiboPure ™ RNA Purification 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 kit (Applied Biosystems) in 20 µl of reactions, following the manufacturer's instructions. The cDNA samples were loaded into a Dynamic Array 48.48 chip in a BioMark system (Fluidigm; San Francisco, CA, USA) through an integrated fluidic circuit controller following a protocol previously described 19 .
For this experiment, the expressed levels of 48 genes were analyzed: 44 target genes and four reference genes (ACTB, B2M, HPRT1, and TBP). The ACTB and TBP were the two most stable reference genes and were used to normalize the expression levels of the target genes. Primers used for the analysis were designed using PrimerExpress 2.0 software (Applied Biosystems) and are detailed in Supplementary Table S12. Data was collected using the Fluidigm Real-Time PCR analysis software 3.0.2 (Fluidigm) and analyzed using the DAG expression software 1.0.5.5 75 applying the relative standard curve method. Samples targeted in this study were analyzed in duplicate. The normalized quantity (NQ) values of each sample and assay were used to compare our data. Data obtained were normalized by performing log 2 transformation of the NQ value. The sex effect was also tested by using a linear model with the lm function of R program 76 . Gene expression association analysis. In order to detect SNPs associated with gene expression levels, eGWAS was performed using as phenotype the expression values of 43 genes in adipose tissue. A mixed model (1) was employed in Qxpak 5.0 77 : ijk i j l l k ijkl in which Y ijk was the k th individual expression record of each one of the 43 analyzed genes, sex (two levels) and batch (corresponding to the slaughterhouse group; five levels) were fixed effects, λ l was a −1 (aa), 0 (Aa), +1 (AA) indicator variable depending on the k th individual's 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 with random distribution N(0, Aσ u 2 ) where A is the numerator of the pedigree-based relationship matrix and the e ijkl the residual.
The R package q-value 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 77,78 , meaning that 5% of all significant p-values will be false positives.
The identified SNPs were classified as cis when they were at 5 Mb upstream or downstream of the gene position and as trans when they were located elsewhere in the genome. Significant associated SNPs located less than 10 Mb apart were considered as belonging to the same interval.
Gene annotation and functional classification. The significantly associated SNPs (FDR ≤ 0.05) identified were mapped in the Sscrofa10.2 assembly and were annotated with the Ensembl Genes 84 Database using VEP software 79 .
The genomic eQTL intervals were annotated using Biomart software [http://www.biomart.org] adding a genomic region of one Mb from the start and end of each trans-eQTLs. Only the studied candidate gene was annotated for cis-eQTLs.
The Core Analysis function included in IPA was used to perform the functional analysis of genes mapped in the 19 eQTL regions. This software was used for data interpretation in the context of biological processes, pathways and networks. All information generated is derived from the Ingenuity Pathway Knowledge Base (IPKB), which is based on functions and interactions of genes published in the literature. RNA-Seq data of BF from BC1_ LD individuals 17 and Gene Expression Atlas 80 were used to determine which of the lipid-related genes annotated in the genomic eQTL intervals were expressed in adipose tissue. Finally, a prediction analysis of TF binding sites was performed in the promoter region of the annotated genes. The iRegulon Cytoscape plugin 65 was used to analyze the TFs and their related target genes. iRegulon relies on the analysis of the regulatory sequences around each gene, and use databases of nearly 10,000 TF motifs and 1,000 ChIP-seq data sets or "tracks". In the present study, a high NES threshold value (NES > 5) was considered for the selection of potential relationships (the default NES cut-off in iRegulon is 3). A high NES score indicates (1) a motif that recovers a large proportion of the input genes within the top of its ranking and (2) higher reliability.

Gene expression and correlation analysis with phenotypes. A correlation analysis was performed
to explore the association among the expressions of the 43 genes. Moreover, to analyze the relationships between gene expression and phenotypes, Pearson's correlations among gene expression and FA composition percentages in BF 6 and Longissimus dorsi muscle 5 were performed. If required, data was normalized applying the log 2 transformation. Then, gene expression and the composition of FA in BF and muscle were corrected both by gender (two levels) and batch effects (five levels), and the composition of FA traits was also adjusted for carcass weight, using the glm R package 81 . 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 82 was used to visualize the results of both traits and genes.