Comparative Transcriptome Profile between Iberian Pig Varieties Provides New Insights into Their Distinct Fat Deposition and Fatty Acids Content

Simple Summary Iberian pigs are meat quality models due to their high fat content, high intramuscular fat, and oleic fatty acid composition. These parameters present great variability and are differentiated among the lines that make up the Iberian pig population. However, there is little information on how the genetic expression influences quality across Iberian varieties. This study aimed to compare the muscle expression profile between two varieties of Iberian pig (Torbiscal and Retinto) and their reciprocal crosses, differentiated by fatness. Our results suggest that the Retinto variety, which has the greatest fat content amongst the studied Iberian varieties, showed a higher expression of genes related to adiposity. Likewise, a higher expression of genes related to lipolysis was found in the Torbiscal variety, described as having less fat content than Retinto. Further genetic variation analysis in these Iberian varieties showed relevant associations for SNP (Single Nucleotide Polymorphism), related to these differentially expressed genes, with the meat quality traits. Thus, our findings evidence that differences in the genetic architecture and expression of Iberian varieties might explain the variability in their fat content and composition and hence, their meat quality. Abstract The high deposition of intramuscular fat and the content of oleic fatty acid are characteristic of the Iberian pig. These two parameters present great variability and are differentiated amongst the varieties that make up the Iberian pig population. Although previous studies generated evidence for causal genes and polymorphisms associated to the adipogenic potential of the Iberian pig, there is little information about how genetic expression influences this trait’s variability. The aim of this study was to analyses the expression profile between two varieties of Iberian pig (Torbiscal and Retinto) and their reciprocal crosses differentiated in their intramuscular fat (IMF) content and fatty acid (FA) composition in the Longissimus thoracis muscle using an RNA-seq approach. Our results corroborate that the Retinto variety is the fattiest amongst all studied varieties as its upregulated genes, such as FABP3 and FABP5, SLC27A1 and VEGFA among others, contribute to increasing adiposity. In its turn, Torbiscal pigs showed an upregulation of genes associated with the inhibition of fat deposition such as ADIPOQ and CPT1A. Further genetic variation analysis in these Iberian varieties showed relevant associations for SNP located within the differentially expressed genes with IMF and FA content. Thus, the differences found in the genetic architecture and the muscle transcriptome of these Iberian varieties might explain the variability in their fat content and composition and hence, their meat quality.


Introduction
The Iberian pig has always been a meat quality model due to its high fat content and its fatty acids composition [1]. In fact, Iberian products are famous and appreciated for their flavor, which depends on the volatile components produced by fatty acid's oxidation [2]. The quality of these products is influenced both by the genetic background of the Iberian breed and its nurture [1][2][3][4].
The Iberian pig population is structured in seven genetic varieties, which differ in their fatty acid profiles and intramuscular fat (IMF) content, as well as in their morphology and productivity [5,6]. Fat deposition is a heritable trait in Iberian pigs and these differences may be influenced by the genetic background of each variety [6]. For instance, the Torbiscal (TT) variety has a lower IMF content in the Longissimus thoracis muscle than Retinto (RR). Moreover, the Torbiscal variety has lower oleic and linoleic fatty acids content than Retinto [5,6]. Additionally, the reciprocal cross between both varieties (RT and TR) reveals a negative heterosis for IMF [6].
The genetic differences between Iberian pig varieties have been tackled by multiple approaches, some of which reported causal genes and polymorphisms contributing to pig breed differentiation [7,8]. However, there is scarce information on how gene expression influences the adiposity of Iberian pig varieties. Given the developments in next generation sequencing technologies, we sought to explore the muscle transcriptome differences between two Iberian pig varieties and their reciprocal crosses, which show distinct IMF content and fatty acid composition in Longissimus thoracis [6]. Using an RNA-seq approach, we identified those genes and metabolic pathways underlying the variation of their fat content and composition, and hence, their meat quality.

Ethics Approval
All animal procedures were carried out in accordance with the Spanish Guidelines for Animal Experimentation and were approved by the Ethical Committee of Institut de Recerca i Tecnologia Agroalimentàries (IRTA-2012-0054-C02-01).

Phenotypic Records and Experimental Design
Twenty-eight castrated male Iberian pigs from the Retinto (RR) (n = 7) and Torbiscal (TT) (n = 7) varieties and their reciprocal crosses Retinto × Torbiscal (RT) (male × female) (n = 7) and Torbiscal × Retinto (TR) (n = 7) were used. The two varieties used in this study are recognized in Spain's official Iberian herd book (Spanish Association of Iberian Purebred Pig Breeders, AECERIBER). During the experiment, pigs were kept from birth to slaughter under intensive rearing conditions such as those used in pig commercial farms, in the CAP-IRTA center (Pig Control and Evaluation) of Monells (Girona, Spain). All of them were reared in similar conditions to reach 102.8 ± 6.8 kg body weight (BW) at 242 ± 12.0 days of age. Then, each variety and cross were reared indoors ad libitum. After this fattening period, pigs were slaughtered in the commercial slaughterhouse "Jamón y Salud" of Llerena (Badajoz, Spain) at 299.3 ± 12.1 days of age and 153.5 ± 10.4 kg BW. For each animal, a sample of 200 g of Longissimus thoracis muscle was collected and stored at −32 • C until analyzed. The total IMF and the fatty acid profile were quantified according to the method described by [9].

Sample Extraction and Sequencing Process
Twenty-eight tissue samples (seven per variety) were taken for RNA sequencing from the Longissimus thoracis muscle of the animals reared ad libitum. RNA was isolated by the acid phenol method using TRI-reagent (Sigma-Aldrich, Tres Cantos, Spain) and following the manufacturer's instructions [10]. Then, the samples were pair-end sequenced using an Illumina Hiseq 2500 platform. On average, the number of reads was 11.79 ± 0.79 millions of reads per sample.

Quality Control
The quality control of the reads was performed using FastQC v11.6 and MultiQC v1.5 [11][12][13]. Alignment and counting processes were carried out using STAR v2.4.0.1 package [14]. MultiQC software was used to visualize the reads before and after the quality control, and after aligning and counting the high-quality sequences.

Differential Expression Analysis
The differentially expressed gene (DEG) analysis was performed with the R package edgeR [15]. Previously, genes whose counts were less than 10 reads were removed to prevent confounders in the normalization. Therefore, 16,252 protein-coding genes were analyzed in this study. Gene counts were normalized using the TMM (trimmed mean of M-values) method, which takes the sequencing depth and the total counts per sample. This procedure prevents sample collection and sequencing errors [14]. In the next step, for each gene, a negative binomial regression was fitted: where Y ij is the normalized expression level of each gene, L i is the variety of each sample, P ij is the covariate, which corresponds to the BW at slaughter of each pig, and ε ij is the residual term. The statistical test used in this analysis was the empirical Bayes (EB), which enables contrasting the expression of each gene between varieties [14]. We implemented multiple test corrections using false discovery rate (FDR) [16]. Lastly, genes were considered differentially expressed (DEG) across groups when the fold change > 1.5 and the FDR < 0.1.
Gene ontology and metabolic pathway enrichment analyses were performed using David, GeneCards and the R package enrichR v1.0. [17,18]. Only those pathways with a p-value < 0.1 were considered.

Gene Variation Analysis and Association Study
In this analysis, SNP calling in coding and untranslated (UTR) regions of the top DEG (Table A3) was performed on the mapping files generated by STAR v2.4.0.1 package [14]. The BAM files were processed, and the variants were called using SAMtools v0.1.18 [19]. Only variants with a minimum root mean square (RMS) mapping quality of 30 were selected for further analysis. We also performed a visual inspection of alignments within the genomic regions of the genes most differentially expressed (Table A3) between varieties using the integrative genomics viewer (IGV) [20].
A selection of 12 SNP variants were genotyped in additional RR (44) and TT (74) pigs from an independent experiment [6,7]. The genotype was based on custom-design allelic discrimination assays (Applied Biosystems, Thermo Fisher Scientific, Waltham, MA, USA) in an OpenArray platform (Applied Biosystems). The association analysis between the SNPs and the IMF and their main fatty acid profiles were performed using plink 1.9 [21] and the statistical model included the sex as a fixed effect and age at slaughter as covariate and the two first principal components (PC) based on a previous PC analysis on this population with 50k SNPs [7].

Quality Control
Quality control analysis revealed high-quality reads with a sequencing precision of 99%. The number of reads per sample after quality control did not vary significantly; at first, there were on average 11.74 ± 0.74 million reads and after the quality control it decreased to 11.65 ± 0.73 million. Moreover, the length of the reads was greater than 96 bp which does not differ from the sequencing read size (100 bp), and there was no adapter contamination higher that 0.1% at each sample. Lastly, the GC content of these samples (40-60%) was within the accepted CG content interval in Iberian pig [22].
The results of the alignment showed that the RNA sequences were correctly aligned with 94% of the sequence length mapped, whereas only 2% of the percentage of sequences aligned to multiple loci and short sequences and, for gene counts, 75% of genes were counted at each sample, and the number of genes mapped to multiple loci and unmapped was approximately 5%.

Differential Gene Expression Analysis
From a total of 16,252 expressed genes, 3799 unique genes were differentially expressed (DEG) across the Iberian pig varieties and reciprocal crosses (RR, TT, RT, and TR) ( Table 1). The RR variety had the greatest number of DEG compared to others (Table 1). In particular, the biggest divergences were found between RR and TR varieties followed by the contrast between RR and TT. Otherwise, the comparison of genetic expression levels between the TT, RT and TR resulted in a much lower number of DEG (Table 1). These results are in line with previous phenotypic analysis (6) which indicated a line, maternal/paternal and heterosis effects between these strains (RR and TT) for IMF and fatty acids composition. Additionally, despite of the small sample size (seven individuals) of our phenotypic data, a difference in IMF was found between RR and TR (−3.11, p-value 0.03) which could explain the greatest number of DEG found between these two crosses. This result would also suggest a possible paternal imprinting of the Torbiscal variety, but it cannot be validated with the current analysis and is falling outside the scope. Nevertheless, the DEG and phenotypic results showed that the RR variety is clearly distinct from the other variety and reciprocal crosses [6]. Thus, given the distinctive nature of the RR variety, our study has focused on analyzing the expression divergence between RR and the rest of Iberian genetic types analyzed in this study.

Gene Ontology and Pathways Enrichment Analysis
Gene ontology and pathways enrichment analysis identified 16 metabolic pathways significantly associated with an adjusted p-value < 0.1, most of which are related to the lipid metabolism (Table A2). In Table 2, we show three metabolic pathways with the highest number of DEG whose function is associated with lipids absorption, fatty acids synthesis and catabolism, and adipocyte differentiation. The influence of these three metabolic pathways in pig IMF content was previously reported in other independent studies that compared the fatness and fatty acid composition of different pig populations [22][23][24]. Interestingly, 30 DEGs, from a total of 220 genes associated with these lipid metabolic pathways, promoted lipid absorption, synthesis, and fat cell differentiation. In fact, most of these 30 genes were upregulated in the RR variety, except for the ADIPOQ, CPT1A and SREBF1 genes, with some exceptions (Table A2). CPT2; ACAA1; ACADS * Number of differentially expressed genes found in the pathway; + Bonferroni correction (P.Val adj); # Acronyms of the differentially expressed genes in each pathway; $ Gene Ontology (GO).
These 30 genes can be assembled in multiple gene groups, such as ACAD, ACAA, CPT, COX and FABP. The first four gene families are involved in the fatty acid betaoxidation, either in the first reactions, carrying these molecules through the different membranes of the mitochondrion (CPT), or in the latest reactions (ACAD, ACAA and COX). Specifically, within these gene groups, we found that ACADL, ACADS, ACADVL, ACAA1, ACAA2, COX1, COX2, COX3, COX4I1, COX4I2, COX5B, COX6A1, COX6A2, COX7A1, COX7A2, COX7B, COX7C, CPT1B and CPT2 genes were upregulated in the RR variety, when compared to the other genetic types analyzed, while CPT1A was downregulated. Knockout studies using Acadl tm1Uab and Acadvl tm1Vjeque mice showed that the lack of the coding proteins of ACADL and ACADVL, respectively, caused lipidosis, hypoglycemia and an elevated content of fatty acids in serum [25,26]. Interestingly, the CPT family includes the downregulated CPT1A gene, which encodes a carrier protein that controls the entrance of long fatty acids in the mitochondrion, promoting the beta-oxidation. Still, DE analysis showed that CPT1B and CPT2 genes were upregulated in the RR variety. These genes have an analogous function in the mitochondrion although CPT1B is specifically located in the muscle. Previous studies suggested that the expression differences between analogous genes, CPT1A and CPT1B, had an impact in the glucose and fat metabolism [27,28]. Thus, these results could explain why RR variety had a higher content of IMF, PUFA and SFA as their gene expression profile promotes reduced fatty acid catabolism [6].
The FABP family is responsible for uptaking fatty acids through the cell membrane of different tissues. Previous studies showed that genes within this family (FABP3, FABP4 and FABP5) play a key role in the fatness differences between pig lines [24]. Interestingly, FABP3 and FABP5 genes were upregulated in the skeletal muscle of RR animals. These genes have a high affinity towards long fatty acids, like the oleic acid (C18:1, n9), which is one of the biochemical characteristics of the RR pigs [6,23]. In addition, a mice knock-out model, Fabp5 tm1Hota , showed that this gene regulates adipocyte differentiation and reduces triglyceride concentration in plasma [29]. This result is in line with published comparative studies in "Duroc × white pig varieties" and "Iberian × Landrace cross" which reported that the upregulation of FABP family, was associated with higher IMF and backfat content, respectively [22,24,[30][31][32]. Lastly, multiple sequence variants within the FABP3, FABP4 and FABP5 genes and near genomic regions have been associated with pig fatness, specifically with increasing FA and IMF content [23,32,33]. Thus, the upregulation of FABP3 and FABP5 genes in RR variety supports its higher fat content within its meat products and endorses that these genes are good candidates for meat quality traits, such as IMF and backfat content.
In addition to these families, the DGE analyses identified additional genes playing a relevant role in lipid metabolism. LPL, PLA2G7 and SLC27A1 genes carry lipids through tissues, cells, or organelles. Our DEG analysis reports that these three genes are upregulated in RR as compared to TT, RT, and TR. LPL participates in lipid's regulation, such as HDL biosynthesis, fat content distribution and as a lipid carrier between lipid and carbohydrate metabolisms [34]. PLA2G7 promotes the transport of oleic fatty acid, which is responsible for the exceptional quality of the Iberian products. The different expression of PLA2G7 may explain the meat quality variations amongst these Iberian pigs. Lastly, SLC27A1 is a long-chain fatty acid membrane transporter that is active in many cell types, and it is highly expressed in pig's gluteus medius, Longissimus dorsi, diaphragm and heart muscles [35]. Previous studies with knockout mice for SLC27A1 gene, Slc27 a1tm1Jkk , showed that its deficit is linked with a decrease in the intramuscular fat deposition [36]. In addition, this gene was found to be upregulated in an RNA-Seq experiment comparing Iberian pigs and other white pigs (i.e., Duroc) [23]. SHDB and VEGFA are upregulated genes in the RR variety associated with the beta-oxidation of fatty acids and with angiogenesis, migration, and cellular growth [37], respectively. Different studies using knockout mice for the SHDB gene reported that the lack of its expression activates the energy metabolism, producing obesity resistance and the development of thyroid lipomatosis [38]. VEGFA gene has been widely studied due to its influence in many cellular processes, such as nutrient diffusion and angiogenesis. Previous studies reported that its upregulation is associated with a higher fat content in pig varieties [22,31]. Thus, the high expression of SLC27A1 and VEGFA in Retinto variety might cause a major absorption and accumulation of IMF and subcutaneous fat in those pigs, respectively, which coincides with the phenotypic characteristics of this variety [6].
Our differential expression analyses also showed other downregulated genes in RR variety, such as ADIPOQ and SREBF1. ADIPOQ encodes a specific hormone of the adipose tissue, which controls energy homeostasis, glucose, and lipid metabolism, and SREBF1 is a transcription factor responsible for lipid's homeostasis. Previous experiments reported that the expression of the ADIPOQ hormone controls multiple metabolic routes in different cell and tissue locations. As an example, the upregulation of ADIPOQ has been associated with insulin-sensitizing and anti-lipotoxic phenotypes in pigs [39,40]. A knockout mice study (Adipoq tmPesch ) showed that the lipid's metabolism, especially beta oxidation, is deactivated and there was a reduction in the insulin sensitivity [41]. In pigs, published studies suggest that the differential expression of ADIPOQ and its receptors in two pig lines with different IMF content, such as "Landrace × Meishan" [42], "low × high IMF Duroc pigs" [43,44] and "Duroc × Landrace × Yorkshire" vs. Laiwu [45], is associated with differences in its IMF deposition. Specifically, these studies show that the expression level of the ADIPOQ gene is lower in those pig varieties or crosses with a higher IMF content. Then, the downregulation of the ADIPOQ gene in RR might be associated with an increase in fat deposition, whereas its upregulation in TT, RT, and TR genotypes might lead to a lower IMF content and to upregulating metabolic pathways associated with FA oxidation [46], which agrees with the upregulation of CPT1A in these varieties. Hence, the characteristic expression level of the ADIPOQ gene in these two varieties agree with their adipogenic potential, as the cells in RR are increasing its fat content and TT cells are stimulating FA metabolism.
By clustering these DEGs according to their metabolic function and cellular location (Figure 1), we show how the muscle of RR pigs upregulates genes and metabolic pathways leading to an increased fat content. Specifically, RR upregulates genes involved in lipid transport and absorption from the plasma to the cell and to the organelles, like the endoplasmic reticulum and the mitochondrion, where the fatty acids are unsaturated and elongated, or oxidized.

Gene Variation Analysis and Association Study
Sequence analysis of the 30 DEGs identified 12 SNPs in UTR and coding regions (Table A3). These SNPs were genotyped in additional RR (44) and TT (73) pigs of an independent experiment with animals from the same populations (see [6] for more details) in which Longissimus thoracis IMF and their fatty acid profiles were also recorded. The association analysis showed that rs326546232 and rs346045742, SNPs located in ADIPOQ and ACAT1 genes, respectively, had a significant association with meat quality traits (Table 3), particularly with saturated fatty acids like the stearic fatty acid, monosaturated fatty acids, like palmitoleic, and polyunsaturated fatty acids, like the linolenic acid. These SNPs were also associated with other fatty acids like oleic, the ratio between saturated and monosaturated fatty acids as well as the intramuscular fat.
Interestingly, RR pigs had distinctive allele frequencies in all the studied SNPs comparing with TT, (Table A3). The heterogeneity in these regions of the Iberian pig varieties genome might explain the diversity in the meat quality traits. Then, these results and its associations with the studied meat quality traits support our differential expression analysis previously described as it corroborates that the adipogenic potential of these Iberian genetic types is influenced by its genetic expression and architecture.

Gene Variation Analysis and Association Study
Sequence analysis of the 30 DEGs identified 12 SNPs in UTR and coding regions (Table A3). These SNPs were genotyped in additional RR (44) and TT (73) pigs of an independent experiment with animals from the same populations (see [6] for more details) in which Longissimus thoracis IMF and their fatty acid profiles were also recorded. The association analysis showed that rs326546232 and rs346045742, SNPs located in ADIPOQ and ACAT1 genes, respectively, had a significant association with meat quality traits (Table 3), particularly with saturated fatty acids like the stearic fatty acid, monosaturated fatty acids, like palmitoleic, and polyunsaturated fatty acids, like the linolenic acid. These SNPs were also associated with other fatty acids like oleic, the ratio between saturated and monosaturated fatty acids as well as the intramuscular fat.
Interestingly, RR pigs had distinctive allele frequencies in all the studied SNPs comparing with TT, (Table A3). The heterogeneity in these regions of the Iberian pig varieties genome might explain the diversity in the meat quality traits. Then, these results and its associations with the studied meat quality traits support our differential expression analysis previously described as it corroborates that the adipogenic potential of these Iberian genetic types is influenced by its genetic expression and architecture.  In summary, the profile of upregulated and downregulated genes in RR are related with an increasing fat absorption, deposit and metabolization in the muscle cells (Figure 1). In addition, SNPs located within multiple DEGs, which play a key role in the lipid's metabolism, are more polymorphic in this variety and are also associated to the fatty acid profile of intramuscular fat, contributing to the nutritional value of meat. Thus, the genetic architecture and expression profile of this variety agrees with phenotypic studies that established that RR has the greatest adipogenic potential of the Iberian varieties [6].

Conclusions
The results of this experiment support that RNA-seq is a useful technique to elucidate which genes and pathways influence phenotypic differences between populations. In this study, we corroborate that Retinto pigs have the greatest fat content comparing to Torbiscal variety and their reciprocal crosses. This finding is supported by its gene expression profile as the upregulated and downregulated genes identified in this variety are associated with an increasing transport, absorption, differentiation, and accumulation of fatty acids in Retinto pigs. Thus, variation in gene expression profiles between these two varieties may explain the differences in their adipogenic potential. Remarkably, the different regulation of ACAA, ACAD, CPT1B, LPL, FABP3, FABP5, SLC24A1, PLA2G7, SREBF1, SHDB, VEGFA and ADIPOQ expression promotes the activation of anabolic routes in RR pigs, and the higher activation of beta-oxidation and catabolic routes in TT pigs.
Further gene variation and association analysis showed that the SNPs rs326546232 and rs346045742, located within ADIPOQ and ACAT1 DE genes, respectively, are associated with all studied meat quality traits, like the saturation of fatty acids content which is related to the nutritional value of meat. In addition, these SNPs showed distinctive allele frequencies between RR and TT varieties, evidencing the polymorphic and then more diverse meat characteristics of Retinto.
Hence, variations in Iberian meat quality traits are influenced by the gene expression and genetic architecture of each breed variety. Although our study evidence that upregulated and downregulated genes and SNPs explain part of the meat quality traits variance between Iberian pig varieties, further studies are required to confirm that meat quality differences are mainly caused by these reported genes and genetic variants.

Data Availability Statement:
The results from all data analysed during this study are included or displayed in this article (and its Tables A1-A3, such as the Metabolic pathways, Fold changes of Differential expressed genes and Associated SNPs). However, other individual results may be made available from the corresponding author on request.
Appendix A Table A1. Metabolic pathways associated with lipid metabolism that include differentially expressed genes (DEG) of our study, adjusted p-value by Bonferroni correction (P.Val adj), number of DEGs in each pathway and DEG names.