Global analysis of the association between pig muscle fatty acid composition and gene expression using RNA-Seq

Fatty acids (FAs) play an essential role as mediators of cell signaling and signal transduction, affecting metabolic homeostasis and determining meat quality in pigs. However, FAs are transformed by the action of several genes, such as those encoding desaturases and elongases of FAs in lipogenic tissues. The aim of the current work was to identify candidate genes, biological processes, and pathways involved in the modulation of intramuscular FA profile from longissimus dorsi muscle. FA profile by gas chromatography of methyl esters and gene expression by RNA-Seq were determined in 129 Iberian × Duroc backcrossed pigs. An association analysis between the muscle transcriptome and its FA profile was performed, followed by a concordance and functional analysis. Overall, a list of well-known (e.g., PLIN1, LEP, ELOVL6, SC5D, NCOA2, ACSL1, MDH1, LPL, LGALS12, TFRC, GOT1, and FBP1) and novel (e.g., TRARG1, TANK, ENSSSCG00000011196, and ENSSSCG00000038429) candidate genes was identified, either in association with specific or several FA traits. Likewise, several of these genes belong to biological processes and pathways linked to energy, lipid, and carbohydrate metabolism, which seem determinants in the modulation of FA compositions. This study can contribute to elucidate the complex relationship between gene expression and FA profile in pig muscle.

by genes encoding enzymes like desaturases and elongases in lipogenic tissues 2,15 . Notwithstanding, the regulatory mechanisms beyond these reactions are complex and involve the combination of TFs 2 , as well as several biological processes and pathways. In the context of animal and plant breeding, there is an increasing interest in identifying genes controlling the phenotypic variation of complex traits. That is the case of several studies focused on the genetic basis of intramuscular FA composition in pig muscle across several breeds 8,[16][17][18][19][20][21] .
The aim of this work was to study the association between the porcine longissimus dorsi (LD) muscle FA profile and its transcriptome, focusing on the identification of the most relevant candidate genes, biological processes and pathways related to intramuscular FA composition.

Methods
Animals, sample collection and phenotypic data. A total of 129 animals generated by an experimental backcross named as BC1_DU (25% Iberian and 75% Duroc) were employed. All pigs were maintained under the same intensive conditions and fed ad libitum with a commercial cereal-based diet and free access to water. A more detailed description of the backcross BC1_LD generation, experimental design, animal raising, and feeding is provided in Martínez-Montes et al. 22 Animal procedures were carried out according to the Spanish Policy for Animal Protection RD1201/05, which meets the European Union Directive 86/609 about the protection of animals used in experimentation. This study was conducted in accordance with relevant guidelines and regulations of the animal care and use committee of the Institut de Recerca i Tecnologia Agroalimentàries (IRTA), which adopts "The European Code of Conduct for Research Integrity" and "Good Experimental Practices". Likewise, the experimental protocol was approved by the Ethical Committee of the IRTA. Our study is also reported in full compliance with ARRIVE guidelines (https:// arriv eguid elines. org/). Animals were slaughtered in five batches in a commercial abattoir of Mollerussa (Spain). Samples of longissimus dorsi (LD) muscle were collected, immediately snap frozen in liquid nitrogen and stored at − 80 °C until analysis. In addition, the distribution by sex was 59 females and 70 males, males were not castrated. At slaughter, pigs had an average age of 190 days (range 174-205 days), and 73.70 kg of carcass weight (range 46.10-109.20 kg).
FA composition in the C14-C22 range was determined using a gas chromatography of methyl esters protocol as described by Mach et al. 23 in intramuscular LD muscle (n = 129). In brief, 200 g of LD muscle samples of 129 BC1_DU pigs were homogenized and used to measure the FA profile in duplicate. Additional information on the LD muscle FA composition in BC1_DU population is indicated in Crespo-Piazuelo et al. 20 .The FA composition (n = 17 FAs) was expressed as percentage of total identified FAs. Total percentages of SFA, MUFA, and PUFA were obtained through the sum of the individual FAs (Table 1). FA and metabolic ratios were calculated from the ratio between individual FA percentages as it is shown in Table 1. In addition, we further calculated the following FA metabolic indices: Total RNA isolation and sequencing. Total RNA was isolated from the LD muscle (100 mg) of 129 animals using the the RiboPure™ Isolation kit for High Quality Total RNA (Ambion ® ; Austin, TX) following the manufacturer's recommendations. RNA quantification and purity was done with a NanoDrop ND-1000 spectrophotometer (NanoDrop products, Wilmington, DE, USA). RNA integrity was checked by Agilent Bioanalyzer-2100 equipment (Agilent Technologies, Inc., Santa Clara, CA, USA), using only those samples with an RNA integrity number (RIN) greater than 7 for the RNA-Seq experiment.
Library preparation and sequencing was performed at CNAG institute (Centro Nacional de Análisis Genómico, Barcelona, Spain). For each sample, one paired-end library was prepared using TruSeq Stranded mRNA kit (Illumina, Inc.; San Diego CA, USA). To discriminate among samples, libraries were labeled by barcoding and pooled to be run in Illumina HiSeq 3000/4000 instruments (Illumina, Inc.; San Diego CA, USA). In brief, in this study 2 × 75 bp reads, a mean of 45.09 million of paired-reads per sample, and an average of 90.06% (ranging from 80.51 to 96.09%) of uniquely mapped reads were generated. Bioinformatic analyses. Quality control and basic statistics of reads were performed using FastQC (v0.11.9) 24 and MultiQC (v0.7) 25 programs, respectively. Sequencing reads were mapped employing the STAR software (v2.7.9a) with default parameters 26 , and using the Sscrofa11.1 pig genome assembly as reference. Then, gene expression was quantified using RSEM (v1.2.28) 27 software with default parameters and annotation from pig Ensembl Genes 97. Filtering was performed to keep rows of those genes with at least 129 reads in total, thus, 14,870 genes were retained. Normalization of the read counts was done as indicated below.

Concordance analysis.
In order to have a preliminary quality control based on the biological role of the results of gene expression, a study of overlap between the number of associated genes with each trait was carried out. This analysis was organized by grouping the FA into five categories: (i) saturated FAs, (ii) monounsaturated FAs, (iii) polyunsaturated FAs, (iv) FA ratios, and (v) metabolic ratios and indices. Within each group, we explored the number of genes detected by each trait as well as the intersection size using the ComplexUpset package v1.3.3 32 .
Gene ontology enrichment analysis. Gene Ontology (GO) enrichment analysis with the list of associated genes (by category among the trait class) were implemented using the ClueGO v2.5.8 plugin 33 in Cytoscape v3.9.0 software 34 . All genes expressed in LD muscle were employed as background in the overrepresentation analysis across biological processes, molecular functions and pathways (KEGG). The statistical significance was assessed with a hypergeometric test using Benjamini and Hochberg method 30 for multiple testing correction (BH corrected P-value < 0.05). In addition, a minimum k-score of 0.44 was used, GO tree interval levels set from three to eight, and a minimum of three genes per cluster with at least 4% in selected genes. Results with and without the fusion feature "GO Term Fusion" were generated to evaluate the redundant parent-child terms. To this end, ClueGo output was visualized using an R script that allowed to facet by categories via the ggplot2 package v3.3.5 35 .
Correlation analysis. Pearson correlation analysis was carried out between the gene expression values and the FA phenotypic values using the R cor.test function. Here, the same normalization scale for the variables to be correlated was applied, as it has been described in the previous step of the ELMSeq algorithm. The Benjamini-Hochberg procedure 30 was used to correct the raw P-value corresponding to the correlation coefficients. The goal of this strategy was to provide information about the direction of the association reported by ELMSeq, and therefore to obtain interpretable results regarding positive or negative relationships between the FA composition and the gene expression.
Selection of candidate genes associated with the FA profile in muscle. The list of detected candidate genes were prioritized based on the following criteria: (i) candidate genes delimited into functional categories across biological processes, molecular functions and pathways (i.e., lipid metabolism, carbohydrate metabolism, amino acid metabolism, and nucleic acid metabolism terms); (ii) genes that correlate both specifically and simultaneously with FA phenotypes; and (iii) literature mining plus novelty assessment via Geneshot (including: rare, uncommon, common and very common levels).

Results
Identification of candidate associated genes with FA profile in muscle. In the first step we performed an association analysis using the ELMSeq algorithm to identify genes (n = 14,870 genes) associated with the profile of 36 FA in 129 BC1_DU pigs. We observed a variable number of associated genes across FA (1022 genes in total, Supplementary Table S1). Among the 36 traits, 22 FA were associated with a range of 1-553 genes. For each trait, we estimated the distribution of the number of associated genes across metabolic processes, lipid metabolism, TFs or co-factors, and novel genes (including long non-coding RNAs) (Supplementary Table S2).
Gene ontology enrichment analysis. In the functional analysis, no significant GO terms/pathways were found for the list of associated genes (i.e., 50 in total) with the SFA category. Instead, for the other categories a total of 75 terms (biological process, molecular functions and pathways) were significantly enriched (Fig. 2). The full results (including genes by each GO terms) with and without "GO Term Fusion" can be consulted in Supplementary Table S3. As shown in Fig. 2, terms related to lipid metabolism were found on some categories, for example, "KEGG:03320 PPAR signaling pathway", "GO:0008203 cholesterol metabolic process", "KEGG:00564 Glycerophospholipid metabolism", and "GO:0044539 long-chain fatty acid import". Moreover, GO terms involved in carbohydrate metabolism were also enriched such as "GO:0006086 acetyl-CoA biosynthetic process from pyruvate", "GO:0004738 pyruvate dehydrogenase activity, "KEGG:00020 Citrate cycle (TCA cycle)", and "KEGG:00,620 Pyruvate metabolism". Other GO terms were found related to carbohydrate and lipid homeostasis "KEGG:04310 Wnt signaling pathway" and control metabolic processes "KEGG:04150 mTOR signaling pathway", respectively. In addition, GO terms related to some biological process, such as "GO:0015980 energy derivation by oxidation of organic compounds", "GO:0006091 generation of precursor metabolites and energy", "KEGG:00190 Oxidative phosphorylation", "GO:2001256 regulation of store-operated calcium entry", and "GO:0032543 mitochondrial translation", were also enriched (Fig. 2). Finally, due to the low number of GO terms (n = 2, GO:0050880 and GO:0035150) found for the different metabolic ratios and indices, they were grouped into the single category term of "FA ratios and metabolic indices" (Fig. 2).
Remarkably, when we explored the existence of genes inversely correlated between traits, we observed: (1) genes positively correlated with MUFA but negatively correlated with PUFA (PLIN1, SFRP5, PPP1R1B, and TRARG1), and vice versa (CYCS and TFRC); (2) genes positively correlated with MUFA and metabolic indices but negatively correlated with FA ratios (LEP); (3) genes positively correlated with MUFA and FA ratios but negatively correlated with PUFA (LGALS12); (4) genes negatively correlated with SFA but positively correlated with FA ratios and/or metabolic indices (NMNAT2); and (5) genes positively correlated with SFA traits but negatively correlated with metabolic indices (FBP1).

Discussion
Studies on pork meat nutritive values and quality have received special attention over the last decade. Among other factors, meat nutritive value is determined by FA composition, playing also an important role in meat quality traits. Likewise, gene expression is assumed to be able to control the variation of such traits. In the present work, we aimed to study the association of muscle transcriptome with intramuscular FA profile from LD muscle to identify candidate genes, biological processes and pathways related with FA composition in pigs.
Here, we used the ELMSeq approach 28 , which has been extensively employed in the study of human diseases like cancer or metabolic disorders. However, to the best of our knowledge, this is the first study that reports the use of ELMSeq approach to test the association between muscle transcriptome and FA composition in pigs or other livestock species.
Here, we reported well-known candidate genes, but also promising novel genes associated to the intramuscular FA profile (Supplementary Tables S1 and S2). In order to prioritize genes, the list of candidates was processed through an analytical pipeline including: (i) concordance analysis between the number of associated genes across categories of FAs (Fig. 1), (ii) ontology-based functional classification ( Fig. 2 and Supplementary Table S3), and (iii) correlation analysis of FA phenotypes and respective gene expression of the associated genes (Fig. 3). As results, we highlighted 57 candidate genes including protein-coding, transcriptional regulators, and long-non coding RNAs (Fig. 3). In the following sections, we discuss some of these genes accordingly to their biological function, as well as the novelty of our findings (Supplementary Table S2 HMGCR and SC5D). A detailed examination of the association study revealed that fructose-bisphosphatase 1 (FBP1) showed a synergistic effect (positively associated) with C16:0, C18:0 and SFA traits, but antagonistic (negatively associated) effect with UI trait (Fig. 3). FBP1 catalyzes the hydrolysis of fructose 1,6-bisphosphate to fructose 6-phosphate and inorganic phosphate, acting as a gluconeogenesis regulatory enzyme (GeneCards ID: GC09M094603) 36 . The expression of ELOVL fatty acid elongase 6 (ELOVL6) was highly associated with MUFA (Fig. 3). ELOVL6 encodes an enzyme responsible for condensation reaction, one of the four-step process (condensation, reduction, dehydration, and one further reduction) for elongation of very-long-chain (ELOVL) FAs. Further, ELOVL6 preferentially utilizes SFA and MUFA as a substrate (e.g., C16:0 and C16:1n-7, respectively) 37 . To be noted, previous studies in Iberian × Landrace pigs (BC1_LD) provided evidence of the link between a polymorphism in the promoter region of ELOVL6 (ELOVL6:c.-533C > T) and the percentage of C16:0 and C16:1n-7 FAs in LD muscle and backfat 38 . In addition, the ELOVL6:c.-394G > A polymorphism was suggested as the causal mutation for the QTL on pig chromosome 8 (SSC8) that affects FA composition in BC1_LD pigs 39 . Despite this, our findings did not detect association between ELOVL6 and SFA, and just a suggestive association with C18:1n-9 and the MUFA/PUFA ratio (Supplementary Table S1). Among the genes with lipolytic effect, we detected malate dehydrogenase 1 (MDH1) associated to C18:1n-9, MUFA, and C18:2n-6/C18:3n-3 and ω6/ω3 ratios (Fig. 3). MDH1 encodes an enzyme which catalyzes reversible oxidation of malate to oxaloacetate in many metabolic pathways (including the citric acid cycle), utilizing the NAD/NADH-dependent system (GeneCards ID: GC02P063557) 36 . The function of MDH1 is primarily related www.nature.com/scientificreports/ to the production of aerobic energy for muscle contraction 40 , and it was reported as candidate gene of meat quality traits in Chinese pig breeds 41 . The expression of another gene, lipoprotein lipase (LPL), was negatively associated with C18:1n-7 (Fig. 3). LPL belongs to the PPAR signaling pathway (Fig. 2), and has the dual functions of triglyceride hydrolase and ligand/bridging factor for receptor-mediated lipoprotein uptake (GeneCards ID: GC08P019901) 36 . Therefore, LPL hydrolyzes circulating triglyceride containing chylomicrons and very lowdensity lipoproteins to produce free FA. These free FAs can be assimilated by different tissues, such as muscle and adipose tissue 42 .
Regarding the genes involved in lipogenesis, we detected sterol-C5-desaturase (SC5D) as positively associated with C20:2n-6, and with C18:2n-6/C18:3n-3 and ω6/ω3 ratios (Fig. 3). SC5D (also known as SC5DL) encodes an enzyme of cholesterol biosynthesis (GeneCards ID: GC11P121292) 36 . SC5DL has been found to be one of the downregulated genes related to cholesterol metabolism in various tissues (including muscle) of lambs 43 . In a similar way, 3-hydroxy-3-methylglutaryl-CoA reductase (HMGCR ) was positively associated with C18:2n-6/C18:3n-3 ratio (Fig. 3). HMGCR is a rate-limiting enzyme in cholesterol synthesis (GeneCards ID: GC05P075336) 36 . HMGCR is also related to several biological processes including pyruvate dehydrogenase Figure 2. Functional enrichment analysis of the associated genes with FA profile and grouping by traits in BC1_DU pigs. The barplot with facets (from left to right) shows the enriched GO terms (BH adjusted P-value < 0.05) when grouping by traits. The traits grouped in this analysis are displayed in Fig. 1A-E. Legend: The full names to symbols " > …" are regulation of endoplasmic reticulum unfolded protein response; oxidoreductase activity, acting on the CH-OH group of donors, NAD or NADP as acceptor and establishment or maintenance of epithelial cell apical/basal polarity. GO terms of each individual category are presented in Table S3. www.nature.com/scientificreports/ activity, coenzyme metabolic process and signal transduction. A previous report in a Duroc population provided evidence of the links between the expression of HMGCR gene in muscle with several traits such as carcass lean percentage, C18:0 and C18:2n-6 contents, but also showed a positive correlation with cholesterol-related traits, intramuscular fat (IMF), and C18:1n-9 and C16:0 FA content 44 .
Our results also reported the association of acyl-CoA synthetase long chain family member 1 (ACSL1) and glutamate oxaloacetate transaminase 1 (GOT1) genes with C18:2n-6/C18:3n-3 and ω6/ω3 ratios (Fig. 3). ACSL1 participate in the long-chain fatty acid import and signal transduction biological process (Fig. 2). The protein encoded by ACSL1 is an isozyme of the long-chain fatty-acid-coenzyme A ligase family (GeneCards ID: GC04M184755) 36 . ACSL1 is involved in the synthesis of long-chain acyl-CoA esters, FA degradation and phospholipid remodeling 45 . Likewise, this gene was found to be associated with lipid metabolism and mitochondrial oxidation of FAs in pigs 46 . In relation to FA degradation, Zhou et al. 47 , suggested that GOT1 was crucial for providing oxaloacetate at low glucose levels, likely to maintain the redox homeostasis. In our results, GOT1 was also related to organonitrogen compound biosynthetic process and signal transduction GO terms (Fig. 2).
Genes related to carbohydrate and lipid metabolism. We also identified several genes belonging to carbohydrate and lipid metabolism that were associated with various of the analyzed traits such as PLIN1, CYCS, SFRP5, LEP and PPP1R1B (Fig. 3). Among the well-known genes, we observed perilipin 1 (PLIN1) positively correlated to MUFA traits and metabolic ratios (MUFA/SFA and MUFA/PUFA), and negatively correlated with C20:2n-6 and C18:3n-3 traits (Fig. 3). PLIN1 was also enriched in the PPAR signaling pathway (Fig. 2). PLIN1 codifies for a protein belonging to the family of perilipins, which play a role in regulating intracellular lipid stor- www.nature.com/scientificreports/ age and mobilization 48 . Association between the expression of PLIN1 and porcine IMF deposition and adipocyte differentiation has already been eported 49 . In addition, a lipolytic function of PLIN1 and increased expression in subcutaneous adipose biopsies from Iberian pigs fed a carbohydrate-enriched diet have been reported 50 . In our results, the expression of cytochrome c somatic (CYCS) was negatively correlated with MUFA and positively with correlated PUFA. CYCS encodes a small heme protein that functions as a central component of the electron transport chain in mitochondria (GeneCards ID: GC07M025118) 36 . The expression of leptin (LEP) was significantly correlated with more than one trait (Fig. 3) (i.e., positively correlated with MUFA traits, C16:1n-7/C16:0 and C18:1n-9/C18:0 ratios, and metabolic ratios as MUFA/SFA and MUFA/PUFA), but it was negatively with C18:1n-7/C16:1n-7 ratio. LEP encodes a protein that plays a major role in the regulation of energy homeostasis (GeneCards ID: GC07P128241) 36 . In addition, as a pleiotropic adipocytokine it can regulate several physiologic functions. LEP constitutes a circulating hormone that orchestrates behavioral and metabolic responses to nutrient intake 51 . Furthermore, LEP interacts with other hormonal mediators, regulators of energy status and metabolism (e.g., insulin or glucagon) to regulate growth and reproduction processes 52 . Indeed, LEP can regulates complex biological effects through its receptors. By using the Iberian pig as translational model for studies on metabolic syndrome, type 2 diabetes and nutrition-associated diseases, a polymorphism of the leptin receptor (LEPR c.1987C > T) has been informed to increase insatiability and obesity (i.e., state that in human medicine is called as leptin resistance) 53 , and as such, the Iberian pigs would be resistant to leptin-induced lipolysis 50 . However, in our study the LEPR gene was not associated with any FA traits.
In the present work, we observed the transferrin receptor (TFRC) as negatively correlated with C18:1n-9 and MUFA traits, but positively correlated with C20:2n-6 and C18:2n-6/C18:3n-3 ratio (Fig. 3). TFRC (also known as TFR1) encodes a cell surface receptor necessary for cellular iron uptake by the process of receptormediated endocytosis (GeneCards ID: GC03M196027) 36 . In mammals, TFR1 imports the transferrin-bound iron from the extracellular environment into cells. Interestingly, the intracellular labile free iron is indispensable for lipid peroxidation and ferroptosis execution 54 . On the other hand, protein phosphatase 1 regulatory inhibitor subunit 1b (PPP1R1B) was positively associated with MUFA traits but negatively correlated with PUFA traits (Fig. 3). PPP1R1B (also known as DARPP-32) encode a bifunctional signal transduction molecule, while stimulation of dopaminergic and glutamatergic receptors regulates its phosphorylation and function (GeneCards ID: GC17P039626) 36 . This gene is a potent inhibitor of protein phosphatase 1 (PPP1, previously known as PP1) when phosphorylated at Thr34 by cAMP-dependent protein kinase (PKA) 55 , and this protein could be important in the control of glycogen metabolism and muscle contraction, among other activities. Thus, it suggests that phosphorylation sites are implicated in the fine regulation of DARPP-32 function and a modulation of the regulation of PPP1 via DARPP-32. In addition, PPP1R1B has been reported as differentially abundant in the LD muscle of phenotypically extreme pigs 13 , and positively correlated with the IMF content of LD muscle in Duroc × Luchuan pigs 56 . FA and glucose transport genes. Five out of the 57 candidate genes were functional classified as FA and glucose transporters. Regarding FA transporters, our results suggested a significant association of four members of the solute carrier (SLC) gene superfamily with FA ratios and C20:2n-6 ( Fig. 3). Members of SLC superfamily encode membrane-bound transporters, which play essential roles in transporting a variety of substrates across cellular membranes. These include amino acids, glucose, inorganic cations and anions, FAs and lipids, acetyl coenzyme A, and vitamins, among others 57 . Out of all the SLC members detected, the solute carrier family 27 member 1 (SLC27A1) has also been reported as a transporter of the predominant substrates of long-chain FAs 58 . Regarding glucose transporter, ENSSSCG00000017801 was identified among the novel genes reported in our study, which was significantly associated with MUFA and PUFA traits, and metabolic ratios (Fig. 3). However, ENSSSCG00000017801 is a novel gene recently annotated in pigs and thus, there is limited information in the literature about its functions. Nevertheless, the human orthologous of this novel gene is the trafficking regulator of GLUT4 1 (TRARG1), which is predicted to be involved in endosome to plasma membrane protein transport and glucose import in response to insulin stimulus (NCBI Gene ID: 286,753) 59 . Likewise, TRARG1 is of particular interest as it was previously recognized to be located in the glucose transporter type 4 (GLUT4) storage vesicles, and to positively regulate GLUT4 trafficking or translocation 60 .

Regulators including transcription factor and co-factors.
In order to reinforce the biological significance of our study, we also focused the discussion on regulators including those TFs and co-factors found among the list of candidate genes. Ladybird homeobox 1 (LBX1) was positively associated with C18:1n-9 and MUFA (Fig. 3). LBX1 is a member of the ladybird-like gene family which encodes a homeodomain transcription factor 61 . LBX1 plays a putative regulatory role during the postnatal development of the porcine skeletal muscle in Meishan pigs 62 . In the present study, we also observed association between well-documented co-factors of FA metabolism. Two of such genes were the nuclear receptor coactivator 2 (NCOA2) and lipin 1 (LPIN1), which were positively associated with ω6/ω3 ratio (Fig. 3). NCOA2 encodes a protein that acts as a transcriptional coactivator for nuclear hormone receptors, including steroid, thyroid, retinoid, and vitamin D. In fact, a key role of NCOA2 as modulator of the intramuscular FA composition in pigs has been reported 63 . On the other hand, LPIN1 encodes a regulatory enzyme that catalyzes the penultimate step in triglyceride synthesis, including the dephosphorylation of phosphatidic acid to yield diacylglycerol (GeneCards ID: GC02P011677) 36 . He et al. 64 reported a link between a polymorphism located in LPIN1 gene with the percentage of leaf fat and IMF in pigs.
LGALS12 was found to perform a critical www.nature.com/scientificreports/ role in lipid metabolism in mice, functioning as an intrinsic negative regulator of lipolysis, regulating lipolytic protein kinase A signaling by acting upstream of phosphodiesterase activity to control cAMP levels 65 . Lastly, Wu et al. 66 reported that LGALS12 knockdown could inhibit adipogenesis in porcine adipocytes by downregulating lipogenic genes and activating the PKA-Erk1/2 signaling pathway.

Conclusions
Taken together, our results identify candidate genes linked to intramuscular FA composition in muscle, including well-known and novel candidate genes involved in biological processes and pathways mainly related to energy, lipid, and carbohydrate metabolism that appear to be determinant in the modulation of intramuscular FA profile in pigs.

Data availability
All relevant data produced or evaluated in this research are disclosed in the paper as well as its supplementary information files. The RNA sequencing data used and analyzed in the current study are available from sequence read archive (SRA), NCBI BioProject under the accession number PRJNA882638 (https:// www. ncbi. nlm. nih. gov/ sra).