Dissecting the Gene Expression Networks Associated with Variations in the Major Components of the Fatty Acid Semimembranosus Muscle Profile in Large White Heavy Pigs

Simple Summary The amount and fatty acid composition of intramuscular fat are important features for the qualitative characteristics of processed and fresh meat products, but the knowledge of the key molecular drivers controlling these traits is still scant. To this aim, the present study investigated the co-expression networks of genes related to variations in the major fatty acids deposited in pig Semimembranosus muscle. Palmitic and palmitoleic acid contents were associated with a downregulation of genes involved in autophagy, mitochondrial fusion, and mitochondrial activity, suggesting that the deposition of these fatty acids may be enhanced in muscles with a reduced mitochondrial function. A higher proportion of oleic acid and a reduction in the percentages of n-6 and n-3 polyunsaturated fatty acids were related to changes in the mRNA levels of genes involved in Mitogen-Activated Protein Kinase (MAPK) signaling. The obtained results indicated gene expression networks and new candidate genes associated with the studied traits. Further studies are needed to confirm our results and identify in the discussed genes molecular markers for future selection schemes aimed at improving pork nutritional and technological quality. Furthermore, as pigs are considered reliable animal models for several human conditions, the obtained results may also be of interest for improving the knowledge of the molecular pathways associated with obesity and diabetes. Abstract To date, high-throughput technology such as RNA-sequencing has been successfully applied in livestock sciences to investigate molecular networks involved in complex traits, such as meat quality. Pork quality depends on several organoleptic, technological, and nutritional characteristics, and it is also influenced by the fatty acid (FA) composition of intramuscular fat (IMF). To explore the molecular networks associated with different IMF FA compositions, the Semimembranosus muscle (SM) from two groups of Italian Large White (ILW) heavy pigs divergent for SM IMF content was investigated using transcriptome analysis. After alignment and normalization, the obtained gene counts were used to perform the Weighted Correlation Network Analysis (WGCNA package in R environment). Palmitic and palmitoleic contents showed association with the same gene modules, comprising genes significantly enriched in autophagy, mitochondrial fusion, and mitochondrial activity. Among the key genes related to these FAs, we found TEAD4, a gene regulating mitochondrial activity that seems to be a promising candidate for further studies. On the other hand, the genes comprised in the modules associated with the IMF contents of oleic, n-6, and n-3 polyunsaturated FAs (PUFAs) were significantly enriched in Mitogen-Activated Protein Kinase (MAPK) signaling, in agreement with previous studies suggesting that several MAPK players may have a primary role in regulating lipid deposition. These results give an insight into the molecular cascade associated with different IMF FA composition in ILW heavy pigs. Further studies are needed to validate the results and confirm whether some of the identified key genes may be effective candidates for pork quality.


Introduction
Global meat demand is estimated to be 16% higher in 2025 than in the 2013-2015 period, with poultry and pork representing the most consistent meat production and demand increase in developing countries [1]. Italy is among the main producers of processed meats, contributing about one-third of the European charcuterie [2] and being the first producer in Europe for Protected Designation of Origin (PDO) certified meat products [2]. The contribution of Italy to PDO production is relevant, particularly for high-quality dry-cured hams (such as Parma and San Daniele PDO hams) [3]. These products are obtained from hind legs from heavy pigs, which are specifically selected to fit the objectives required for these products [4,5], including slaughter at a minimum age of nine months and about 150-170 kg live weight. The genetic types used for heavy pig production show specific genetic and phenotypic features; carcasses and cuts have a balanced deposition of lean mass and fat, which make these meats suitable for processing.
The amount and fatty acid (FA) composition of subcutaneous and intramuscular fat (IMF) is a significant feature for the qualitative characteristics of processed and fresh meat products. Indeed, dry-cured ham quality greatly depends on thigh covering fat and IMF deposition, as an appropriate fat layer and higher amount of IMF are required for improved organoleptic characteristics and have been proven to prevent hams from experiencing excessive processing losses [4]. The amount and composition of IMF are of great importance for the organoleptic and nutritional quality of fresh meat, affecting also consumers' perception and acceptance of pork products [6]. Even though an increase in the degree of lipid unsaturation in meat could be beneficial for human health [7,8], polyunsaturated FAs (PUFAs) are more likely to incur oxidative phenomena [9], making these FAs undesirable for the processing industry. Therefore, the pork processing industry's technological requirements and consumers' dietary demands do not completely match. These contrasting requirements create the necessity to find a balance between these divergent needs. So far, efforts have been made to improve both technological and nutritional quality traits, and an extensive body of literature focuses on the multiple variation factors affecting meat quality and FA composition of the final meat products [5,6,9,10]. Studies carried out in Italian Large White (ILW) pigs found that the FA composition of IMF can be improved through selection [11] and is associated with genetic markers [12]. In fact, muscle FA composition showed low to medium heritability estimates, ranging from 0.157 to 0.237 in ILW pigs. In particular, n-6 FA showed positive genetic correlations with carcass lean mass, while being negatively related to backfat thickness, suggesting that FA deposition in muscle is also related to carcass composition [13]. However, the exact molecular basis of muscle FA composition is still relatively unknown.
In the last ten years, high-throughput sequencing technologies have been used for transcriptome analysis, which has allowed the exploration of the gene co-expression networks in an unprecedented manner in terms of accuracy and data insight [14]. To date, the molecular networks associated with muscle FA composition have been investigated in several pig populations [15][16][17], but no information exists about gene co-expression networks related to FA composition in ILW heavy pigs. In Landrace x Iberian crossbred pigs slaughtered at about six months of age, individuals with high PUFA deposition in the Longissimus dorsi muscle were characterized by transcriptomic profiles, indicating an inhibition of glucose uptake and lipogenesis, which would produce a decrease in the muscle triglyceride storage [15]. In a sample of White Duroc × Erhualian pigs, muscle and liver FA composition was shown to be associated with genes involved in insulin resistance and in the immune system [16]. The latter result was in agreement with the literature, where adipogenesis has been indicated in humans and pigs to be implicated in inflammatory responses [18][19][20], suggesting that studying fat traits in pigs may also provide insights for a better comprehension of human disease genetic background.
In humans, several pathological conditions, such as obesity and diabetes, have been associated with changes in the FA composition of muscle and adipose tissue of the patients [21][22][23][24]. In humans, low contents of long-chain PUFAs and high amounts of palmitic acid in muscle were found to be linked with obesity [24]. On the other hand, high contents of n-3 and n-6 PUFA in adipose tissue were related to smaller adipocyte size [23]. Thus, as suggested by Zhang et al. [16], reaching a better understanding of the gene networks associated with FA composition in porcine tissues could provide information concerning lipid metabolism that may be of interest also for a better understanding of obesity-related pathological conditions in humans.
Therefore, shedding light on the genetic factors underlying muscle FA composition may be of paramount importance for the pig production industry and may also provide useful information for other species. Our previous study reported the gene expression networks associated with Semimembranosus muscle (SM) IMF deposition and the genes differentially expressed in two groups of ILW heavy pigs divergent for SM IMF content [13]. The current study starts from the results described in our previous research [13] and aims at investigating the co-expression networks of genes related to the variability noticed for the major FAs (palmitic, palmitoleic, stearic, and oleic acids) and FA classes (n-6 and n-3 PUFAs) in pig SM. These results allow for the identification of genes associated with FA composition and suggest possible genetic markers for the improvement of meat nutritional and technological quality.

Sampling and Phenotypes
The present research was conducted on the same set of samples already used in our previous paper [13]. This set of samples consisted of 12 individuals, chosen from a purebred population of 949 sib-tested ILW pigs, for their extreme and divergent contents of SM IMF (low IMF vs. high IMF group). A detailed description of the population of origin and of the 12 samples was reported in Zappaterra et al. [11] and Zappaterra et al. [13], respectively. The pigs entered the testing station at about 30 kg live weight and were reared in the same environmental conditions. During the testing period, pigs were single-stabled and fed the same finishing diet at a quasi ad libitum feeding level (i.e., about 60% of pigs were able to ingest the entire supplied ration). The sib test ended when the pigs reached an average final live weight of about 150 kg at about eight months of age. Animals were then transported to a commercial slaughterhouse according to Council Rule (EC) No. 1/2005 on the protection of animals during transport and related operations. The 949 pigs were slaughtered on 27 different days; after being electrically stunned, the animals were bled in a lying position in agreement with Council Regulation (EC) No. 1099/2009 on the protection of animals at the time of the slaughtering. After slaughter, SM samples were collected from the same carcass side of the 949 ILW pigs. The SM samples were then immediately frozen in liquid nitrogen and stored at −80 • C in a deep freezer until further analysis. As previously described [25], IMF content was determined for all gathered samples by extraction with petroleum ether from 1 g of SM using an XT15 Ankom apparatus (Ankom, Macedon, NY, USA), according to Official Procedure AOCS Am 5-04 [26]. IMF was determined in % (g/100 g of muscle tissue). SM FA composition was determined as described in Zappaterra et al. [11]. The FA and FA classes were expressed as mg/g of IMF. To avoid as much as possible other confounding effects, the two groups of samples (low and high IMF groups) were obtained balancing for sex and avoiding full and half-sibs. The 12 samples were slaughtered on nine different days.

RNA Extraction, Library Preparation, and Sequencing
A detailed description of the RNA extraction, library preparation, and sequencing is reported in Zappaterra et al. [13]. After total RNA extraction and quality checks were performed, stranded total RNA libraries were prepared, converting RNA to cDNA and adding sequencing adapters. RNA libraries were analyzed with a paired-end sequencing strategy on Illumina HiSeq2500 (Illumina Inc., San Diego, CA, USA). This strategy consisted of extracting, through ultra-high-throughput sequencing, short reads from both ends of cDNA fragments. The libraries were tagged, and pairs of libraries were run on a single lane. Paired-end reads of 100 bp were generated. The obtained raw sequence data have been deposited in the Gene Expression Omnibus (GEO) expression database under the accession number GSE144780.

Mapping, Assembly of the Reads, and Production of the Gene Count Matrix
Raw reads were obtained in FASTQ format and were quality-assessed using the FastQC program (retrieved from URL: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/, accessed on 18 March 2020). As reported in Zappaterra et al. [13], the quality of the raw reads was measured according to sequence-read lengths and base-coverage, nucleotide contributions and base ambiguities, and quality scores and over-represented sequences. Trimmomatic utility [27] was used to trim out terminal low-quality bases and adaptor sequences, according to the FastQC reports. In particular, sequences with a quality below 15 in a sliding window of 4 consecutive bases were filtered out, and the remaining reads were retained if they had at least a 36-nucleotide length. The splice-aware read mapper HiSat2 [28] was used to align the clean reads against Ensembl reference genome Sus scrofa v. 11.1 (retrieved from URL: http://igenomes.illumina.com.s3-website-us-east-1.amazonaws. com/Sus_scrofa/Ensembl/Sscrofa11.1/Sus_scrofa_Ensembl_Sscrofa11.1.tar.gz, accessed on 18 March 2020). After alignment, one BAM file was obtained for each couple of fastq files. BAM files were further processed with Stringtie [29], with the aim of assembling the sequences into known transcripts. HTSeq version 0.6.1 [30] was then used to quantify the reads and obtain the file with the gene counts. The raw gene counts were then normalized with regularized-logarithm transformation (i.e., rlog) using the DESeq2 package [31] in the R environment [32]. The normalized gene counts were used to perform a Principal Component Analysis (PCA) in the R environment [32] to evaluate the dataset multivariate structure.

Weighted Gene Correlation Network Analysis (WGCNA)
The genes whose expression was not detected in most of the samples or with variance zero were filtered out from the normalized gene count matrix. To identify strong coexpression networks of genes related to SM FA composition, we employed a co-expression analysis approach using the "WGCNA" package [33] in the R environment [32]. The gene count matrix was analyzed together with the IMF contents of palmitic, palmitoleic, stearic, oleic, n-6, and n-3 FAs. These FAs were chosen for their relative abundance in IMF (palmitic, palmitoleic, stearic, and oleic are the saturated FAs (SFAs) and monounsaturated FAs (MUFAs) most represented in pork fat tissues [6,11]) and for their nutritional value (n-6 and n-3 PUFAs).
To obtain the scale-free undirected co-expression networks between the genes in the gene count matrix, an adjacency matrix was built using Pearson's correlations between each gene couple. The soft threshold power values were estimated through the function pickSoftThreshold() in the WGCNA package, and a soft threshold power (β) of 6 was used to raise the adjacency matrix. This value was chosen as the scale-free topology index (R 2 ) reached the peak (R 2 > 0.70) for the first time when β = 6, and the minimum module size was 30 genes (Supplementary Figure S1). After β was determined, the adjacency matrix was calculated using the topological overlap measure (TOM) and the corresponding dissimilarity (dissTOM = 1 − TOM). The latter was used as a distance for gene hierarchical cluster, and then DynamicTree Cut algorithm [33] was used to identify the modules of genes. The default minimum cluster merge height of 0.25 was retained. After the construction of the gene co-expression network, WGCNA restituted a list of gene modules, which were named using color labels and clustered highly interconnected genes.
The principal component of each module was defined as the module eigengene (ME); MEs represented the expression value of each module and were used to detect modules that may comprise genes having a biologically relevant role in the variations of the traits. The module-trait relationship (module membership, MM) was calculated as the Pearson's correlation between the ME and the traits of interest. For each gene, the gene MM and the relative p-values were obtained and carefully evaluated, as they indicated the importance of that gene in a module. Due to the high number of gene modules showing a significant correlation with the analyzed traits, we decided to choose the gene modules for further analysis based on the threshold reported by Pampouille et al. [34]. Based on this study, only the gene modules with an absolute value of module-trait correlation higher than 0.7 were further analyzed. After selecting the gene modules exceeding the set threshold, the genes participating in those modules with a MM p-value < 0.001 were further considered to perform the functional enrichment analysis. The intramodular connectivity and gene significance (i.e., a coefficient and relative significance value indicating the importance of that gene for the trait, based on the correlation between the gene expression profile and the trait) were used to identify key genes in the networks.

Functional Enrichment Analysis of the Genes in the Significant Modules
To explore the biological functions of the genes significantly entering (MM p-value < 0.001) the identified gene modules (module-trait correlation higher than |0.7|), we performed Gene Ontology (GO) term enrichment analysis, as well as pathway ontology analyses using the Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 online tool [35] and Cytoscape version 3.8.2 [36]. Cytoscape functional enrichment analysis was obtained using ClueGO [37] and CluePedia [38] plug-ins. Homo sapiens was used as the reference organism for both DAVID and Cytoscape functional enrichment analyses. Bonferroni adjusted p-values < 0.05 were considered significant for DAVID results, and Bonferroni step-down adjusted p-values < 0.05 were considered significant for Cytoscape functional enrichment analysis. For a better representation of the Cytoscape results, we applied the GO term fusion, in order to reduce redundant terms [37].

Mapping Results and Weighted Gene Co-Expression Network Analysis
As we detailed in our previous study [13], a total of 3,235,579,132 paired-end reads were obtained from the SM transcriptome sequencing of the 12 samples. After trimming and filtering steps, 1,155,025,434 reads were retained and mapped against the Sus scrofa reference genome 11.1. The total % of mapped reads (unique and multi-reads) was between 90% and 94%, and the uniquely mapped reads ranged between 79.1% and 85.7% across samples (Supplementary Table S1). The mapped reads were then assembled into transcripts, and the normalized gene counts were used to assess the multivariate structure of the dataset with a PCA. The results of the PCA showed no clear clusterization between samples (Supplementary Figure S2).
The gene count matrix with the rlog values was then submitted for weighted gene coexpression analysis with the WGCNA package. After filtering out 5554 genes with missing expression values in most of the samples or with zero variance, a total of 20,046 genes were used for weighted gene co-expression network construction. As a result of network construction, 113 modules (i.e., cluster of co-expressed genes) were found by WGCNA analysis. These modules were identified by different color names, and their size ranged from 3374 to 30 genes with ME p-value < 0.05. The correlations between the 113 gene modules and the SM traits are displayed with a heatmap in Figure 1. The heatmap shows the strength of the correlations between modules and traits: the darker (dark blue and dark red) the colors, the stronger the correlations. As can be noticed, palmitic and palmitoleic acids show similar patterns of colors (and thus similar correlations with the modules) when compared to IMF. Oleic acid shows opposite correlations when compared to n-3 and n-6 FA. In fact, the modules that were positively correlated with oleic acid were negatively related to the percentage of n-3 and n-6, and vice versa. On the other side, stearic acid showed patterns of correlations different from the ones noticed for the other traits and of weaker magnitude, as indicated by the lighter colors found in the stearic acid column. The complete list of the module-trait relationship values and the relative significances are reported in Supplementary File S1. The strongest and most significant moduletrait correlations are reported in Table 2. The modules grey60, skyblue1, and laven-derblush3 were already discussed in our previous study on IMF deposition [13]. These The complete list of the module-trait relationship values and the relative significances are reported in Supplementary File S1. The strongest and most significant module-trait correlations are reported in Table 2. The modules grey60, skyblue1, and lavenderblush3 were already discussed in our previous study on IMF deposition [13]. These modules showed significant module-trait correlations also with some of the FAs considered in the present study, but additional modules were identified. Palmitic acid was negatively correlated with lavenderblush3 (r = −0.80; p-value = 0.002) and with steelblue (r = −0.80; p-value = 0.002). The same two modules were also correlated with palmitoleic, but with correlation coefficients of lower magnitude (r = −0.60 and −0.65; p-values = 0.038 and 0.022, respectively). Palmitoleic acid showed only one module (yellowgreen, r = 0.71; p-value = 0.009) above the threshold of |0.70|. Oleic, n-6, and n-3 FA were correlated with the same modules (i.e., skyblue1, antiquewhite2, darkorange, darkseagreen3, darkgreen, and steelblue). In some cases, these modules had module-trait correlation values ≥ |0.70| for all three traits (i.e., antiquewhite2, darkorange, darkgreen); in other cases, the modules had module-trait correlation values ≥ |0.70| for some FAs among oleic, n-6, and n-3 PUFAs but were significant (p-value < 0.05) also for all of them (i.e., skyblue1, darkseagreen3, and steelblue). Stearic acid did not show module-trait correlations above the set threshold. As can be noticed from Table 2, the considered FAs were clustered into two groups: on one hand, palmitic and palmitoleic acids were associated with the same gene modules; on the other hand, oleic acid was associated with the same gene modules found for n-3 and n-6 PUFAs, but with opposite correlation coefficients. That is to say, genes whose expression was positively correlated with n-6 and n-3 PUFAs (such as those included in the darkorange and darkgreen modules) were negatively correlated with oleic acid deposition, and vice versa. For these reasons, two functional analyses were performed: one for the genes included in the modules mostly associated with palmitic and palmitoleic, and the other for the genes included in the modules related to oleic acid, n-3, and n-6 PUFAs.
The complete list of the gene MM and the relative p-values (indicating the importance of each gene in the identified modules) are reported in Supplementary File S2. The overall gene significance of each gene on the studied traits is reported in Supplementary File S3.

Functional Enrichment Analyses
The genes significantly entering in the modules most correlated with the FA composition were then submitted for functional analysis. Supplementary File S4 reports the functional categories identified with DAVID for the genes significantly entering in the modules associated with palmitic and palmitoleic acids (i.e., yellowgreen, lavenderblush3, and steelblue). Figure 2 shows the Cytoscape functional annotation enrichment for the same genes. The genes included in yellowgreen, lavenderblush3, and steelblue modules were significantly enriched in "GO:0016482 Cytosolic transport" (p = 0.01), "GO:0048312 Intracellular distribution of mitochondria" (p = 0.02), "GO:0001970 Positive regulation of activation of membrane attack complex" (p = 0.02), "GO:0010508 Positive regulation of autophagy" (p = 0.02), and "GO:0010594 Regulation of endothelial cell migration" (p = 0.04).
On the whole, both functional analysis tools suggested that the genes included in yellowgreen, lavenderblush3, and steelblue modules were significantly enriched in pathways controlling the organization of intracellular organelles. Supplementary File S5 reports the functional categories identified with DAVID for the genes significantly entering in the modules associated with oleic acid, n-6, and n-3 PUFAs (i.e., skyblue1, antiquewhite2, darkorange, darkseagreen3, darkgreen, and steelblue). Figure 3 shows the Cytoscape functional annotation enrichment for the same genes. These genes were significantly enriched in "GO:1905897 Regulation of response to endoplasmic reticulum stress", "GO:0004712 Protein serine/threonine/tyrosine kinase activity", "GO:1901844 Regulation of cell communication by electrical coupling involved in cardiac conduction", "GO:1903598 Positive regulation of gap junction assembly", "GO:0031098 Stress-activated protein kinase signaling cascade", "GO:0023014 Signal transduction by protein phosphorylation", "GO:0071900 Regulation of protein serine/threonine kinase activity", "GO:0043405 Regulation of MAP kinase activity", "GO:0004672 Protein kinase activity", "GO:0018209 Peptidyl-serine modification", "GO:0046777 Protein autophosphorylation", "GO:0018105 Peptidyl-serine phosphorylation", and "GO:0004674 Protein serine/threonine kinase activity" (all p < 0.001). On the whole, both functional analysis tools suggested that the genes included in skyblue1, an-tiquewhite2, darkorange, darkseagreen3, darkgreen, and steelblue modules were significantly enriched in pathways involved in serine/threonine kinase activity. Supplementary File S5 reports the functional categories identified with DAVID for the genes significantly entering in the modules associated with oleic acid, n-6, and n-3 PUFAs (i.e., skyblue1, antiquewhite2, darkorange, darkseagreen3, darkgreen, and steelblue). Figure 3 shows the Cytoscape functional annotation enrichment for the same genes. These genes were significantly enriched in "GO:1905897 Regulation of response to endoplasmic reticulum stress", "GO:0004712 Protein serine/threonine/tyrosine kinase activity", "GO:1901844 Regulation of cell communication by electrical coupling involved in cardiac conduction", "GO:1903598 Positive regulation of gap junction assembly", "GO:0031098 Stress-activated protein kinase signaling cascade", "GO:0023014 Signal transduction by protein phosphorylation", "GO:0071900 Regulation of protein serine/threonine kinase activity", "GO:0043405 Regulation of MAP kinase activity", "GO:0004672 Protein kinase activity", "GO:0018209 Peptidyl-serine modification", "GO:0046777 Protein autophosphorylation", "GO:0018105 Peptidyl-serine phosphorylation", and "GO:0004674 Protein serine/threonine kinase activity" (all p < 0.001). On the whole, both functional analysis tools suggested that the genes included in skyblue1, antiquewhite2, darkorange, dark-seagreen3, darkgreen, and steelblue modules were significantly enriched in pathways involved in serine/threonine kinase activity.

Discussion
The purpose of this study was to investigate the gene co-expression networks related to the palmitic, palmitoleic, stearic, oleic, n-6, and n-3 PUFA percentages in ILW pigs' SM. This study starts from the results of Zappaterra et al. [13], where we investigated, in the same samples used in the present research, the differentially expressed genes and the gene networks associated with the divergent contents of SM IMF.
Palmitic, palmitoleic, and oleic acids showed trends of gene module-trait correlations similar to those observed for IMF content, suggesting that these traits may share, at least in part, a common genetic background. This evidence is consistent with our previous study [11], where positive genetic correlations were identified between these three FAs, and between these FAs and SM IMF content in ILW pigs. In agreement with these observations, oleic acid content was correlated with IMF deposition also in other porcine breeds and muscles. The content of oleic acid showed, indeed, positive genetic correlations with the IMF deposition in Longissimus dorsi [39], Semimembranosus [39], and Gluteus medius muscles of Duroc pigs [39,40]. Furthermore, the gene modules mostly associated with IMF and discussed in our previous study [13] also showed significant positive correlations with some FAs considered in the present study. In Zappaterra et al. [13], the gene networks considered for further discussion were those comprised in the modules skyblue1, darkturquoise, lavenderblush3, and grey60, which contained genes regulating DNA transcription and cell differentiation, primary cilia morphogenesis, ERK/MAPK, and G protein signaling cascades. In the current study, the module skyblue1 was associated with oleic and n-6 PUFAs, lavenderblush3 was significantly related to palmitic acid proportion, and grey60 showed positive module-trait correlations with palmitic and palmitoleic acid contents (whose value, however, did not exceed the threshold of |0.70|). This overlapping between the modules associated with IMF and the FAs considered in this study seems therefore to suggest that the same gene networks involved in IMF deposition may also play a role in the muscle metabolism of palmitic, palmitoleic, oleic, n-6, and n-3 PUFAs. In addition to those modules, the present study permitted the identification of new gene networks associated with the considered FAs, suggesting that FA composition variability relies only in part on the gene networks associated with IMF deposition.
On the whole, palmitic, palmitoleic, and oleic acids showed similar module-trait correlations when compared with IMF, while n-3 and n-6 PUFAs showed opposite trends in gene module-trait correlations when compared to those observed for IMF content. This result suggests that the gene networks linked to increased IMF deposition are associated with decreased depositions of n-3 and n-6 PUFAs. The negative relation between fat deposition and FA unsaturation degree is well known in the literature. PUFAs are indeed mainly incorporated in membrane phospholipids and play a role in membrane flexibility, inflammation control, eicosanoid production, plasma triacylglycerol synthesis, and gene expression regulation [41]. For these reasons, in individuals fed the same diet, PUFA absolute content remains quite stable in animal tissues over time and does not depend on fat deposition [42]. On the other side, increased fat deposition implies greater synthesis and storage of triacylglycerols, which are mainly constituted by SFAs and MUFAs, determining an increase in SFAs and MUFAs and a decreased proportion of PUFAs of the total FAs [42]. The trends in gene module-trait correlations observed in the present study are, therefore, in agreement with the knowledge concerning the FA synthesis and metabolism in animal tissues.
Intriguingly, unlike the above-mentioned FAs, stearic acid did not show modules with correlation values above the threshold of |0.70|. In animal tissues, stearic acid may either be obtained from the diet or synthesized by the elongation of palmitic acid. Even when compared with other SFAs, such as palmitic and myristic acids, stearic acid was found to be a less efficient substrate for triglyceride synthesis [43]. Its content in porcine subcutaneous tissue was found to increase between early life and slaughtering (at 90-110 kg live weight), suggesting that stearic acid deposition increases as subcutaneous fat becomes thicker [44]. When compared to subcutaneous fat, IMF is less dependent on diet and develops in a later stage of life [45,46], suggesting that stearic acid content may be less variable in muscle compared to subcutaneous fat depots. Furthermore, this FA is also the major substrate for the enzyme stearoyl-CoA desaturase (SCD), which catalyzes the conversion of stearic acid to oleic acid [43]. Thus, the fact that none of the observed gene modules obtained for stearic acid exceeded the defined threshold may depend on the role of this FA as an intermediate compound that undergoes several subsequent conversion steps through different metabolic pathways. The lack of strong correlations between stearic acid and the gene modules in the present study may, however, also be determined by the fact that stearic amounts in porcine tissues have smaller variation ranges when compared with other FAs [47].
Another objective of the current study was to find gene networks involved in the deposition of the considered FAs and identify possible genes of interest for the improvement of muscle FA composition. The percentage of palmitic and palmitoleic acids in IMF was associated with genes involved in intracellular distribution of mitochondria and mitochondrial activity, regulation of endothelial cell migration, positive regulation of autophagy, positive regulation of membrane attack complex, and cytosolic transport. The observed relation between the mitochondrial activity and the content of palmitic acid has already been reported in the literature. Palmitic acid is indeed a major FA used for energy supply, via adenosine triphosphate (ATP) production through FA oxidation in mitochondria. However, an excessive amount of palmitic acid was proven to cause lipid accumulation and induce endoplasmic reticulum stress and mitochondrial dysfunction in hepatic and pancreatic cells [48,49]. This evidence suggests that palmitic acid functions both as a substrate and as a regulator of mitochondrial activity. The literature reports that this FA is known to play a pivotal role in the occurrence of mitochondrial dysfunction, lipotoxicity, and cell death [50]. However, full agreement on the negative effect of this FA is lacking in the literature, as a study on rat liver mitochondria showed that palmitic acid may not cause cytotoxicity or cell death, and that these effects may be triggered by high amounts of palmitoleic acid [51]. In the present study, we found palmitic and palmitoleic acid contents associated with the same gene co-expression networks, suggesting that the deposition of these two FAs was associated with similar molecular pathways. The genes involved in mitochondrial activity and distribution were less expressed in the samples with the highest contents of palmitic and palmitoleic acids. In particular, the expression of Clustered Mitochondria Homolog (CLUH) and Microtubule Associated Protein Tau (MAPT) genes, reported in the literature as genes involved in mitochondrial activity and distribution, had strong and negative gene significance values with palmitic and palmitoleic acids. This result suggests that these genes were less expressed in pigs displaying higher SM contents of palmitic and palmitoleic acids, which may indicate a lower mitochondrial activity in those samples, possibly leading to an increased deposition of palmitic and palmitoleic acids and IMF. Though our hypothesis relies only on gene expression results, and further investigations are needed to confirm it, the role of mitochondrial function in the regulation of IMF deposition in pigs is well established [52]. However, as palmitic acid can be both a substrate and a regulator of mitochondrial activity, we are not able to discern whether the altered mitochondrial function is the leading cause of palmitic accumulation or whether the increased palmitic acid accumulation has led to a lower expression of genes related to mitochondrial function. The hypothesized altered mitochondrial function in muscle samples with high palmitic acid seems also to be supported by the downregulation of Mitofusin 2 (MFN2) gene in the same samples. This gene encodes a mitochondrial membrane protein that is critical for mitochondrial fusion and mitochondrial clustering, and its downregulation has been found to cause mitochondrial dysfunction and altered Ca2+ homeostasis in neurons [53]. Interestingly, deficiencies in mitochondrial fusion have been found to disrupt the normal metabolic routes involved in FA oxidation, rerouting towards intracellular lipid droplets the FAs that should have been used in FA oxidation [54]. During starvation, normal FA trafficking and oxidation inside the cell may be mediated by two mechanisms: one requires autophagic digestion of lipid droplets; the other is by lipolytic consumption of lipid droplets [54]. The first involves autophagosomal engulfment of the lipid droplets and fusion with the lysosome, where hydrolytic enzymes digest the lipids stored in the lipid droplets, releasing free FAs into the cytoplasm that are then oxidized by mitochondria [55]. Intriguingly, in the present study, we found several genes involved in the autophagy regulation being downregulated in samples with high palmitic acid. In particular, in addition to MFN2, high palmitic and palmitoleic acid contents in SM were associated with a downregulation of the genes La Ribonucleoprotein 1, Translational Regulator (LARP1), RNA Binding Motif Protein 14 (RBM14), and Serine/Threonine Kinase 11 (STK11), which are all genes whose expression activates autophagy [56][57][58]. Also, the gene Leucine Rich Repeat Kinase 2 (LRRK2), upregulated in the samples with high IMF, palmitic, palmitoleic, and oleic acids, has been implicated in the literature as a modulator of the autophagy process [59]. Autophagy is a complex process playing a role of paramount importance in cell survival. This process is indeed required for balancing sources of energy at critical times in cell development, eliminating misfolded or aggregated proteins (aggresomes), clearing damaged organelles (such as mitochondria), and removing intracellular pathogens [60]. Recent studies highlighted the important role of autophagy in the regulation of lipid metabolism in the liver [61], and our results seem to suggest that some sort of autophagy regulation on lipid metabolism may also exist in muscle. In the samples with high palmitic and palmitoleic acid percentages, the under-expression of genes involved in the autophagy process may suggest a downregulation of the autophagy processes in those samples. This cue may be consistent with a possible decrease in autophagosomal engulfment of lipid droplets in samples with a high percentage of palmitic and palmitoleic acids. It is, however, not possible to state with certainty whether the observed downregulation of the genes related to autophagy is actually associated with lipid droplets utilization, and the suggested hypotheses would need further verification. It is worth noting that the gene with the most negative gene significance for palmitic acid found in the present study is TEA domain transcription factor 4 (TEAD4), a transcription factor preferentially expressed in skeletal muscle that regulates the expression of muscle-specific genes [62]. At present, this gene is poorly investigated in livestock species and has been only cited once in a recently published paper investigating gene expression networks involved in cattle feed efficiency [63]. Intriguingly, this gene is highly conserved in mammals and controls mitochondrial transcription: a loss of TEAD4 impairs the recruitment of mitochondrial RNA polymerase to mitochondrial DNA (mtDNA), resulting in a reduced expression of mtDNA-encoded electron transport chain, with impaired oxidative energy metabolism in mammalian cells [64]. Based on our results and the evidence in the literature, it is possible to hypothesize that TEAD4 downregulation noticed in samples with high palmitic acid deposition may be associated with a possible reduction in mitochondrial FA oxidation. Though these hypotheses need further investigation to be supported, TEAD4 has been recently indicated as a hub gene involved in feed efficiency traits in Nelore cattle [63], thus suggesting that this gene may play a role of paramount importance in the regulation of energy metabolism in livestock species.
The modules associated with the IMF contents of oleic, n-6, and n-3 were extremely linked to Mitogen-Activated Protein Kinases (MAPKs). This result agrees with the literature, where MAPK-related pathways have often been related to FA synthesis. In neoplastic cell cultures, increased MAPK levels were related to an upregulation of Sterol Regulatory Element Binding Protein 1 (SREBP-1) transcription factor levels and of Fatty Acid Synthesis (FAS) enzyme, ultimately causing an increase in FA synthesis [65]. Among MAPK members, in the present study, Mitogen-Activated Protein Kinase 15 (MAPK15, also known as ERK7 or ERK8) gene was found to be under-expressed in samples with high oleic acid and overexpressed when samples had high percentages of n-6 and n-3 PUFAs. Remarkably, this kinase was found to act as a suppressor of adipogenesis in Drosophila [66], suggesting an anti-anabolic role for this kinase. Our result thus supports the evidence in the literature, as MAPK15 was under-expressed in the samples with a higher amount of oleic acid, whose content in muscle is positively related to IMF [11]. Consistently, this gene was upregulated in the samples with high percentages of n-3 and n-6 PUFAs, whose proportion is higher in samples with low IMF deposition [11,42]. MAPK15 has also a role of primary importance in primary cilium formation in human cells, and its depletion disrupts the organization of these cellular organelles [67]. This result seems to support the interpretations reported in our previous study, where we hypothesized a role for primary cilia in adipogenesis, and several genes related to primary cilia morphogenesis were found differentially expressed [13]. Together with MAPK15, the expression of the gene Mitogen-Activated Protein Kinase Kinase 4 (MAP2K4) was also negatively related to oleic acid and positively correlated with PUFAs. This gene suppressed the levels of peroxisomal proliferator-activated receptor γ (PPARγ) [68], a transcription factor highly expressed in adipocytes [69]. PPARγ is known to be a master regulator of adipogenesis and thermogenesis and an active coordinator of lipid metabolism and insulin sensitivity [69]. Despite the fact that its expression is particularly high in subcutaneous and visceral fat, lower levels of this gene were also found in skeletal muscle, where its mRNA expression was higher in the Longissimus dorsi muscle of pigs with high IMF contents [70]. Taken together, these results support the hypothesis that MAPK plays a role in the control of fat, and consequently oleic acid, deposition in SM, suppressing adipogenesis and therefore increasing the proportion of n-6 and n-3 PUFAs in the total FAs. Among the genes negatively related to oleic acid content is also Apolipoprotein A-5 (APOA5), which is widely studied in humans for its role in regulating plasma triglycerides [71]. Of great interest is also the recent identification of the relative APOA5 protein as a major regulator of triglyceride storage in hepatocyte intracellular lipid droplets, suggesting for this gene a role in obesity and metabolic syndromes [72]. A higher expression of this protein was associated in humans with an impaired triglyceride storage [72] and a decreased uptake of FA in cardiomyocites [73]. Mutations in the porcine APOA5 gene sequence were also found associated with meat quality traits in a Chinese pig breed characterized by high IMF deposition [74]. Taken together, these pieces of evidence suggest that this gene may also play a role in muscle intracellular lipid droplets and in triglyceride storage. The increased expression of the APOA5 gene noticed in the samples with low oleic acid (and thus lower IMF contents) may indeed support the hypothesis that, in these samples, triglyceride utilization was promoted by APOA5, which activated lipid droplet mobilization. On the other hand, the gene whose expression correlated the most with oleic acid content was WD Repeat Domain 19 (WDR19, alias IFT144). Its mRNA level was positively correlated with oleic acid and negatively related to PUFA proportions in SM, suggesting that this gene may be also involved in IMF deposition. Once again, this gene was involved in primary cilia morphology [75], emphasizing also in livestock species the possible role of primary cilium-related genes in fat deposition and FA composition. In agreement with this result, Phosphatidylinositol-4-Phosphate 3-Kinase Catalytic Subunit Type 2 Alpha (PIK3C2A) gene was also found among the genes with the highest (and positive) gene significance for oleic acid content. A wide range of biological functions have been attributed to PIK3C2A, including its effects as a regulator of autophagy, glucose transport [76], and primary cilia formation and function [77]. The phosphoinositide 3-kinase-protein kinase B/AKT (PI3K-PKB/AKT) pathway is one of the most important molecular signaling pathways implicated in several cellular processes, including energy metabolism (i.e., glucose uptake), and in cell growth, proliferation, motility, and survival [78]. In our study, the samples overexpressing PIK3C2A showed also a higher expression of Solute Carrier Family 2 Member 4 (SLC2A4, alias GLUT4). SLC2A4 gene encodes a sugar transporter protein that catalyzes hexose transport across cellular membranes through an ATP-independent mechanism and is highly expressed in adipose tissue and skeletal muscle [79]. Its upregulation in the SM samples with a high content of oleic may thus be a direct reflection of the highest content of IMF in those animals, which is highly dependent on the number of muscle inter-dispersed adipocytes [80].
On the whole, the obtained results highlight the gene networks associated with FA variability in porcine SM. Transcriptomic analyses and weighted gene co-expression networks are important tools for disentangling the molecular scenario associated with phenotypic variability and suggest new genes of interest for performing further investigations. Because of the limited sample size, the results should be considered with caution and need further confirmation. Some of the identified molecular networks may indeed be biased by the number of samples, which may have determined a background noise and spuriously correlated gene pairs. Therefore, though the molecular patterns and key genes discussed in the present study seem to be consistent with the literature, further investigations are needed to validate the hypothesized molecular cascades associated with palmitic, palmitoleic, oleic, n-3, and n-6 PUFAs.

Conclusions
To the best of our knowledge, this is the first study investigating the gene expression networks associated with the palmitic, palmitoleic, stearic, oleic, n-6, and n-3 PUFA percentages in ILW pig SM. The obtained results confirmed the results of our previous study on IMF deposition, suggesting a possible involvement of genes regulating primary cilia morphology in muscle FA composition. SM samples with high percentages of palmitic and palmitoleic acids showed an under-expression of genes involved in autophagy, mitochondrial fusion, and mitochondrial activity. Among the key genes related to these FAs, we found TEAD4, a gene regulating mitochondrial activity, which seems to be a promising candidate for further studies aimed at finding the key molecular drivers involved in the control of cell energy metabolism. On the other hand, the genes comprised in the modules associated with the IMF contents of oleic, n-6, and n-3 PUFAs were significantly enriched in MAPK signaling, and several genes encoding for kinases (in particular, MAPK15 and MAP2K4) may be important players regulating lipid homeostasis in muscle tissue. Except for APOA5, we found few genes encoding for enzymes with a direct activity on lipid substrates, suggesting that lipid metabolism and FA composition in muscle are mediated by complex intracellular signaling cascades regulating energy flux and intracellular lipid droplet utilization. The obtained results improve the knowledge of the gene networks involved in muscle FA composition, suggesting the possibility that, if confirmed, the discussed genes may be candidates to be included in selection schemes for pork quality.
Supplementary Materials: The following are available online at https://www.mdpi.com/2076-2 615/11/3/628/s1, Figure S1: Scale independence (on the right) and mean connectivity parameters (on the left) for the different values of Soft Threshold power. The highest R2 for the Scale Free Topology Model is indicated with the red line, and thus the first value of the Soft Threshold power located above the red line is the parameter chosen to obtain the best representation of the weighted gene co-expression network; Figure S2: Score plot of the Principal Component Analysis (PCA) performed on the normalized gene count matrix for the 12 samples of Semimembranosus muscle (SM). Different groups are indicated with different colors, and different shapes are used for the two sexes (F = females, M = barrows); File S1: The complete list of the module-trait relationships and the relative p-values for all identified gene modules and considered fatty acids. Module entries are arranged in alphabetical order; File S2: Gene Module Memberships (MM) and the relative p-values for all genes considered for the weighted gene co-expression analysis. Gene MM values indicate the correlation of the gene expression with the module eigengene and thus the importance of the gene in that module. Gene entries are arranged in alphabetical order; File S3: the list of the overall gene significance and the relative p-values of the genes on the considered fatty acids. The gene significance indicates the correlation between the gene expression and the trait; File S4: Results of DAVID functional enrichment analysis for the genes significantly entering in the gene modules associated with the percentages of palmitic and palmitoleic acids in the Italian Large White pig Semimembranosus muscle; File S5: Results of DAVID functional enrichment analysis for the genes significantly entering in the gene modules associated with the percentages of oleic acid, n-6 and n-3 polyunsaturated fatty acids in the Italian Large White pig Semimembranosus muscle.  Institutional Review Board Statement: The samples used in the present study were obtained from the routine assessments performed during the Sib Test selection schemes of the Italian National Association of Pig Breeders (ANAS). Sampling occurred with the permission of ANAS. The animals used in the present study were reared and slaughtered in compliance with the European rules (Council Regulation (EC) No. 1/2005 and Council Regulation (EC) No. 1099/2009) on the protection of animals during transport and related operations and at the time of the slaughtering. All slaughter procedures were monitored by the veterinary team appointed by the Italian Ministry of Health and were performed within the ANAS routine assessments. For these reasons, the present research did not need approval from a research ethics committee.

Data Availability Statement:
The data presented in this study are openly available in the Gene Expression Omnibus (GEO) expression database under the accession number GSE144780.