Comparison of Muscle Transcriptome between Pigs with Divergent Meat Quality Phenotypes Identifies Genes Related to Muscle Metabolism and Structure

Background Meat quality depends on physiological processes taking place in muscle tissue, which could involve a large pattern of genes associated with both muscle structural and metabolic features. Understanding the biological phenomena underlying muscle phenotype at slaughter is necessary to uncover meat quality development. Therefore, a muscle transcriptome analysis was undertaken to compare gene expression profiles between two highly contrasted pig breeds, Large White (LW) and Basque (B), reared in two different housing systems themselves influencing meat quality. LW is the most predominant breed used in pig industry, which exhibits standard meat quality attributes. B is an indigenous breed with low lean meat and high fat contents, high meat quality characteristics, and is genetically distant from other European pig breeds. Methodology/Principal Findings Transcriptome analysis undertaken using a custom 15 K microarray, highlighted 1233 genes differentially expressed between breeds (multiple-test adjusted P-value<0.05), out of which 635 were highly expressed in the B and 598 highly expressed in the LW pigs. No difference in gene expression was found between housing systems. Besides, expression level of 12 differentially expressed genes quantified by real-time RT-PCR validated microarray data. Functional annotation clustering emphasized four main clusters associated to transcriptome breed differences: metabolic processes, skeletal muscle structure and organization, extracellular matrix, lysosome, and proteolysis, thereby highlighting many genes involved in muscle physiology and meat quality development. Conclusions/Significance Altogether, these results will contribute to a better understanding of muscle physiology and of the biological and molecular processes underlying meat quality. Besides, this study is a first step towards the identification of molecular markers of pork quality and the subsequent development of control tools.


Introduction
Growing market demand for lean meat has directed pig breeding programs to obtain modern meat type of fattener [1]. Intense selection aiming at improving pork production efficiency through increased daily gain and carcass leanness has resulted in improved growth rate and feed conversion ratio as well as lean meat content and loin eye area, and decreased back fat thickness and carcass fat content [2]. However, some meat quality traits playing an important role in consumer acceptance of pork, like water holding capacity, colour, pH, intramuscular fat (imf) content and tenderness, were also affected [3]. Meat quality is complex and depends on the interactive effects of pig genotype, environmental conditions, pre-slaughter handling and slaughtering procedure [4]. Moreover, meat quality determination, as a result of physiological processes taking place in muscle could involve a large pattern of genes associated with both muscle structural and metabolic features. Ascertaining the transcriptome expression profiles differences between selected and non selected breeds which exhibit great differences for muscle meat quality traits, could be helpful to understand the biological processes underlying the development of meat quality.
For this purpose, the experiment was conducted to study gene expression profiles in Longissimus lumborum muscle (LM) of two contrasted pig breeds in terms of carcass fatness and meat quality, Large White (LW) and Basque (B). LW is the most predominant breed used in modern pig industry, with high lean meat productivity, low fat content and high daily gain, but with standard meat quality. By contrast, B is a local, indigenous breed with low lean meat and high fat contents, high meat quality characteristics, and which is genetically distant from other European pig breeds [5,6]. Furthermore, the present transcriptome analysis is the first one undertaken on the high meat quality B breed, despite the increasing number of publications focusing on gene expression in relation with pork quality [7][8][9][10][11].
The aim of our study was to investigate the LM transcriptome profiles of LW (n = 20) and B (n = 20) pigs in relation to muscle traits and meat quality, and thereby clarify the biological events that result in the great phenotypic differences reported in literature between these two breeds [6,12] and improve our general understanding of the determination of pork quality. These two breeds of pigs were reared either in alternative (A, indoor bedding and free access to an outdoor area; n = 10 per breed) or conventional (C, fully slatted floor; n = 10 per breed) housing systems, already demonstrated to influence some muscle and meat quality traits [13].
In order to get accurate information regarding gene expression profiles, the transcriptome analysis was undertaken using a new and specific pig muscle microarray, the 15 K Genmascqchip, in which 85% of the probes have been linked to a unique annotated sequence, and to 9169 unique genes [14]. Functional analysis of Gene Ontology (GO) biological process (BP) terms and functional annotation clustering were undertaken, to highlight main relevant biological networks and genes associated with muscle physiology and meat quality.

Growth, body composition and LM characteristics
As shown in Table 1, B and LW pigs displayed large differences regarding growth, body composition, LM composition and biophysical traits, and sensory quality of meat. At the same live weight, B pigs were older (+85 days) than the LW due to their lower growth rate, and exhibited higher backfat thickness and percentage (+75%) and lower percentage of loin (230%). Regarding LM composition, water content was slightly higher in the LW, whereas total protein content was similar between the two breeds. LM collagen content and glycolytic potential (GP) were higher, but intramuscular fat (imf) content was lower (251%) in the LW compared with the B pigs. LW pigs also exhibited higher meat drip loss and shear force, and lower scores of tenderness, juiciness and flavour than the B. No significant effect of the housing system was found on growth and body composition, as well as on LM composition and biophysical traits. Tenderness score was lower for meat from A versus C housing system (4.0 vs. 4.4, P = 0.018) but juiciness and flavour scores did not differ (P.0.10) according to housing system.

Transcriptomic analysis
Comparison of B and LW muscle transcriptome was achieved using a custom 15 K skeletal muscle pig microarray (Genmascqchip) [14] which is publicly available through Gene Expression Omnibus (GEO) [15] Platform accession no. GPL11016. Briefly, this new porcine skeletal muscle microarray is well annotated (more than 70%) and thereby allows studying a list of 9169 unique genes corresponding to 8622 human Entrez Gene ID. The WEBbased GEne SeT AnaLysis Toolkit [16] was used for the categorization of Gene Ontology (GO) terms for Biological Process (BP). The GO-slim (i.e. representing high-level GO) terms was used to focus on the most important processes. As shown in Figure 1, 13 biological processes were highlighted. The metabolic process category is the most important one (50% of the genes), whereas growth category accounts for less than 5% of the genes, and around 20% of the genes remained unclassified.
Muscle expression profiles of the two breeds (B, n = 20 and LW, n = 20) reared in the two different housing systems (C, n = 10 per breed and A, n = 10 per breed) were compared by transcriptomic analysis. LM genes expression was not modified according to the housing system, since no differentially expressed probe was found between A and C pigs. By contrast, we observed a strong breed effect on gene expression, with 12% of probes being differentially expressed between B and LW pigs (Benjamini-Hochberg (BH) adjusted P value,0.05). Genes showing a significant difference in expression between breeds were divided into 2 lists according to fold change (FC) value. Fold change value is expressed as the expression ratio of B to LW samples when genes are highly expressed in B pigs and as the expression ratio of LW to B samples when genes are highly expressed in LW pigs. The differentially expressed probes corresponded to 1233 unique annotated genes, out of which 635 were highly expressed in the B pigs (Table S1) whereas 598 genes were highly expressed in the LW pigs (Table  S2). Full details of gene name, description, identification, FC and BH adjusted P value are reported in the Tables S1 and S2. In case of redundancy (i.e. more than one probe per gene), the FC were always similar within the probe set suggesting that microarray data were highly consistent. The most differentially expressed (FC.2) and well informative (i.e. with at least one associated GO BP term) genes are shown in Table 2 for genes highly expressed in the B pigs (2,FC#2.6), and in Table 3 for genes highly expressed in the LW

Validation of microarray analysis by quantitative RT-PCR
Among the differentially expressed genes, twelve were chosen to validate the microarray differential expression results by real time quantitative PCR. Six genes highly expressed in the LW breed (ADAMTS8, SPARC, GLOD4, ANKRD1, HHATL and IGF1) and six genes highly expressed in the B breed (LIPE, ZNF24, FOS, FABP3, PPARD and FHL3) with FC extending in microarray analysis between 1.2 and 4.3, were thus analysed by RT-PCR. The results are shown in Figure 2. All these 12 genes were also found differentially expressed between the two breeds by RT-PCR methodology, and for each gene, the FC values were similar between the two methodologies used, i.e. microarray and RT-PCR.

Functional analysis of differential expression between breeds
The two lists of genes, the 635 genes highly expressed in the B pigs and the 598 genes highly expressed in the LW pigs were submitted to an enrichment analysis for GO BP using the Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatic resources [17][18][19]. Significant results (P value#0.05) are presented in Table 4. GO BP terms related to lipid metabolism and transport, carbohydrate metabolism and transcription, were enriched in the B highly expressed genes list. GO BP terms for biological adhesion, protein polymerisation, chemotaxis and cytoskeleton organization, were enriched in the LW highly expressed genes list. To reduce the redundancy and study functionally related genes into a network format, a functional annotation clustering was performed using DAVID tools [17]. We used the three GO terms, BP, cellular component (CC) and molecular function (MF), BIOCARTA (http://www.biocarta.com) and Kyoto Encyclopedia of Genes and Genomes (KEGG) [20] pathways to build on biological modules consisting of clusters of related functional terms, for both B (Table 5) and LW (Table 6) highly expressed genes lists. An enrichment score of 1.3 which is equivalent to non-log scale P-value of 0.05, has been used as threshold for cluster significance according to Huang et al. [19]. Six enrichment groups were found to be significant with an enrichment score higher than 1.3, in the B highly expressed genes list. Three groups of genes were functionally categorized as genes of cytoskeleton (cluster 1-B), vacuole and lysosome (cluster 2-B) and glucose metabolic process (cluster 4-B). The three others clusters (clusters 3-B, 5-B and 6-B) were connected with transcription process. Regarding the LW highly expressed genes list, eight clusters related to extracellular region and collagen (clusters 1-LW, 2-LW and 5-LW), polysaccharide binding (cluster 3-LW), cell motion (cluster 4-LW), contractile fiber and actin polymerization (clusters 6-LW and 7-LW) and chemotaxis (cluster 8-LW) were identified.

Discussion
Regarding genetic background, the B pig is an indigenous breed characterized as ''unique'' among 11 breeds belonging to seven European countries [5]. Despite an increasing number of publications focusing on gene expression in relation with pork quality [11], the present study is the first transcriptome analysis of this non-selected and high meat quality pig breed. Our objective was to clarify the biological events which could enlighten the muscle phenotypic differences reported in the literature between the B and LW pigs [6]. Transcriptional profiling of whole skeletal muscle tissue presents a challenge since changes in gene expression may reflect mRNA composition between various cell types existing in this tissue. However, even if we cannot ascribe expression changes to one specific cellular type we assume that myofiber is the major one and that comparison between breeds in the tissue as a whole is informative. After transcriptome analysis, 12 differentially-expressed genes between LW and B breeds have been properly validated by quantitative PCR analyses, thus demonstrating that the new GenmascqChip [14] is a powerful tool to study pig gene expression and thus get a better understanding of muscle physiology.
Even if the number of studies comparing gene expression of skeletal muscle in different pig genetic backgrounds is rather scarce [10,[21][22], the number of genes found differentially expressed between the two breeds in this study is in the same order of magnitude as found in literature. The high discrepancy in gene expression between B and LW pigs was rather balanced, with 635 and 598 genes highly expressed in B and in LW, respectively. This was associated with strong breed differences regarding growth performance, skeletal muscle characteristics and meat quality traits. In particular, a lower lean and a noticeable higher fat development were observed in B compared with LW, in agreement with previous B and LW comparisons on growth, carcass and muscle traits, and fat tissue metabolism [6,12]. In order to slaughter the pigs from the two breeds at the same time and body weight, B pigs were put on experiment two months earlier and were three months older than LW pigs at slaughter because of their slower growth rate. Moreover, pigs from both breeds received the same amount of feed at a given live weight whereas the potential growth rate and appetite are much higher in LW than in B pigs. Thus, breed, age and feeding effects are confounded.
On the contrary, no difference in LM gene expression profiling was highlighted between A and C housing systems. Accordingly, LM composition, biophysical traits and meat quality were not affected by the housing system, except tenderness. This differs from a previous comparative study [13], thereby confirming that the animal response to husbandry varies according to genotype, environmental (climatic) conditions, etc. [23]. Moreover, differences for tenderness score were much lower between housing systems than between breeds. This supports the generally higher effect of genotype, especially for highly contrasted breeds, than housing conditions on muscle and meat traits [24]. However, we can not exclude that slaughtering conditions could have masked a potential housing effect established before slaughter.
A functional analysis of differential gene expression between LW and B pigs highlighted four main relevant biological networks associated to these breed differences: 1/metabolic processes, 2/ cytoskeleton and contractile fiber, 3/extracellular matrix, and 4/ vacuole, lysosome and proteolysis. Some examples of genes belonging to each of these categories will be discussed in relation to muscle physiology and meat quality. Metabolic processes Enrichment analysis reveals that genes related to lipid metabolism process and fatty acid (FA) transport are more expressed in B than in LW breed. The higher FABP3 (muscle fatty acid binding protein) expression and imf content of the B pigs corroborate several studies indicating this gene as a candidate for the control of imf deposition in pigs [25,26]. In the B pigs, the higher expression of ACACB (acetyl-CoA carboxylase beta) considered as the rate limiting step in FA synthesis, agrees with Alfonso et al. [6] who reported a higher activity of acetyl-CoA carboxylase (ACC) in the muscle of B compared with LW. However, even if the B pigs seemed to deposit more imf than the LW, they also use more lipids as fuel substrates, and rely on fat oxidation and lipolysis to sustain their metabolic requirements. Indeed, PPARD (peroxisome proliferator-activated receptor delta), SLC25A20 (carnitine/acylcarnitine translocase, which mediates the transport of acylcarnitines into the mitochondrial matrix for their oxidation) and ETFDH (electron-transferring-flavoprotein dehydrogenase), all related to the mitochondrial oxidation of FA, were more expressed in B than LW pigs. Regarding lipolysis, PPAP2A (phosphatidic acid phosphatase type 2A) and LIPE were found in the enriched functional category from the lipid catabolic process in the B highly expressed gene list. PPAP2A would play an active role in the hydrolysis and uptake of lipids from extracellular space [27], and a higher expression of LIPE in the ''fatty'' Jinhua than in the leaner Landrace breed has been observed [28]. Altogether, this indicates that a higher FA turn-over (including transport, synthesis and catabolism) could explain the breed discrepancy for imf content. However, since investigation has been conducted on the whole muscle tissue, we cannot exclude that a contribution of a higher number of adipocytes in B pigs could have mediated gene expression variations between the two breeds. Finally, the lower SPARC (secreted protein, acidic, cysteine-rich) expression in the LM of the B pigs is consistent with their higher imf and suggests a role of this gene in controlling imf content.
Indeed, SPARC has been reported to inhibit adipogenesis and SPARC-null mice have been found to exhibit significantly more fat accumulation than wild-type mice [29].
Both functional annotation clustering (cluster 4-B) and enrichment analysis showed that glucose metabolism process is also of great importance in LM traits of B pigs. In this cluster, AGL (glycogen debranching enzyme) and PHKB (phosphorylase kinase beta) are responsible for the complete degradation of glycogen [30,31]. This might suggest that B pigs would use glycogen as a muscle metabolic substrate whereas the LW would spare more glycogen, thus explaining their higher muscle GP. However, PGM1 (phosphoglucomutase 1), GAPDH (glyceraldehyde 3 phosphate dehydrogenase) and LDHA (lactate dehydrogenase A) are up-regulated in the LW, indicating that these pigs would also rely on glucose to fulfil their energy requirements. The higher gene expression of glycolytic pathways in the LW agrees with the more glycolytic and less oxidative muscle metabolism generally observed in domestic compared to wild pigs [32].
Last, creatine kinase (CK) is an essential enzyme to maintain the ATP/ADP ratio in muscle cells and adjust energy availability for contraction. The higher expression of CKB (creatine kinase B chain; cytosolic) in the LW and of CKMT2 (sarcomeric mitochondrial creatine kinase) in the B, suggest a rapid glycolytic ATP production during contraction in LW, while in B the mitochondrial ATP production would be transferred to myofiber via CKMT2, thereby reflecting a more oxidative muscle metabolism. Moreover, in agreement with the suggested cytosolic CK as a candidate protein marker for pork drip loss [33], CKB is more expressed in the LW which exhibited higher drip loss than the B pigs. Accordingly, cytosolic CK protein content was shown to be positively associated with meat lightness, which increases with drip [34].

Cytoskeleton and contractile fiber
Three clusters, cluster 1-B (cytoskeleton) issued from B, and clusters 6-LW (contractile fiber) and 7-LW (actin filament  organization) issued from LW functional annotation clustering analyses, reveal the skeletal muscle organization and structure as important features to characterize the breed differences in gene expression profiles. The actin cytoskeleton is involved in many cellular processes [35], but the relationships between actin dynamics, cytoskeletal organization and muscle development are still unclear. Interestingly, ABRA (actin-binding Rho activating protein, also called STARS, striated muscle activator of Rho signaling) is highly expressed in the LM of B pigs. ABRA activates the serum response factor and leads to enhanced gene expression in skeletal muscle [36] and could thus contribute to the up-regulation of transcription found in the B pigs (clusters 5-B and 6-B). In the same way, FOS, the most differentially expressed gene in the B muscle, is a transcription factor known to induce myogenesis. Thus, ABRA and FOS are probably associated to the higher transcriptional activity observed in the B breed (clusters 3-B, 5-B and 6-B). Apart from transcription, ABRA is involved in skeletal muscle atrophy and hypertrophy [37]. In this cytoskeleton cluster, we also found ABLIM2 (actin binding LIM protein family, member 2) recently identified as an ABRA interacting partner [38]. We can thus hypothesize that ABRA and ABLIM2 could control the development of B muscle and maintain its cytoskeletal integrity. LMOD2 (leiomodin 2) which interacts with actin filaments to promote thin filament elongation and probably their length [39] displays a higher expression in B pigs. This might have led to longer sarcomeres and thereby contributed to improve meat tenderness in this breed, since sarcomere length is positively associated with pork tenderness [40]. Furthermore, because muscles with short sarcomeres generally exhibit high drip loss [41], the higher expression of LMOD2 could also be related to the lower drip loss of the B breed. Besides, MYOZ1 (myozenin 1) also called calsarcin 2, is a sarcomeric calcineurin binding protein specific of striated muscles [42]. Calcineurin mediates calcium signalling and plays a central role in the regeneration and regulation of hypertrophy of skeletal muscle [43]. Calcineurin activity would be inhibited by MYOZ1, as shown in MYOZ1 knock-out mice [42]. Therefore, we hypothesize that the higher expression of MYOZ1 in the B muscle relates to their lower muscle mass through a reduced calcineurin activity, compared with the LW muscle. Similarly, TRIM63 (tripartite motif containing 63, also called MURF1: muscle specific RING-finger protein-1) localized at both M-and Z-lines of the sarcomeres, has been related to muscle atrophy by gene expression profiling and knock-out studies [44]. This would suggest higher muscle atrophy in the B than in the LW, which could be related to their older age at slaughter and might contribute to their lower loin percentage. Finally, in this cytoskeleton cluster, ZYX (zyxin), a protein involved in stress fiber repair and maintenance of cytoskeleton integrity [45] was also found as highly expressed in the B pigs. However, we did not report any positive relationship between ZYX expression and meat drip loss in our study, contrarily to Ponsuksili et al. [8]. This may be explained by different experimental designs, i.e. contrasted breeds in the present work versus extreme drip loss groups within a F2 population in the study of Ponsuksili et al. [8], and indicates that the relationships between ZYX expression and pork quality remain further studies.
Apart from cytoskeleton, our results emphasise the importance of myofibrillar network and especially the contractile fiber (cluster 6-LW) in the muscle expression profile differences between B and LW. Major constituents of sarcomeres: ACTA1 (actin alpha 1), ACTA2 (actin alpha 2), MYH1 (myosin heavy chain 1, IIx), MYH3 (myosin heavy chain 3), TPM1 (tropomyosin 1) and TPM3 (tropomyosin 3) are all highly expressed in the LW, indicating that in this breed, LM is a fast skeletal muscle expressing IIx myosin. In this cluster, NEB (nebulin), which is abundantly expressed in skeletal muscle, plays a key role in thin filament length regulation, intermyofibrillar connectivity and calcium homeostasis [46]. Moreover, in Hanwoo cattle, NEB expression is associated with low marbling and high shear force [47], in accordance with the higher NEB expression in LW than B pigs, and their lower imf content and higher shear force value. The ANKRD1 (cardiac ankyrin repeat domain 1, also known as CARP) which belongs to this cluster, interacts with titin and other sarcomeric proteins to maintain sarcomeric integrity, and its expression is altered in several conditions such as exercise, muscle wasting, dystrophies and stress response [48]. In heart, ANKRD1 interacts with protein CASQ2 (calsequestrin 2) which stores Ca 2+ inside the sarcoplasmic reticulum and modulates Ca 2+ homeostasis [49]. Thus, ANKRD1 seems to be involved in both structure and calcium handling in skeletal muscle, two characteristics of great importance for meat quality. The higher expression of ANKRD1 in LW muscle associated with the great differences in meat quality between the 2 breeds, confirm the involvement of this gene in the biological processes determining pork quality, in agreement with Ponsuksili et al. [50] who suggested ANKRD1 as a candidate gene for meat quality. Interestingly, CASQ2 and ATP2A1 (sarcoplasmic reticulum Ca 2+ -ATPase 1), a protein controlling the pumping of Ca 2+ from the cytosol back to the sarcoplasmic reticulum, are both highly expressed in LW muscle. Besides, correlations between ATP2A1 mutation and imf as well as muscle water content, suggest that ATP2A1 locus could affect pork quality [51]. In conclusion, present results demonstrate the importance of muscle structure (cytoskeleton and sarcomere properties) in the differences found between breeds, thereby confirming the role of structural proteins in the determination of muscle and meat phenotypes [50] even though the relationships between gene expression and muscle traits remain further studies.

Extracellular matrix
The LM of LW pigs is characterized by enrichment clusters of the cellular component GO terms for the extracellular region (cluster 1-LW) and the extracellular matrix part (clusters 2-LW and 3-LW). These clusters (representing 81 genes) exhibit the highest enrichment scores and statistical significance in this study, thus revealing their biological importance. SPARC and SMOC2 (SPARC related modular calcium binding 2), two genes of these clusters, are members of the BM40 family which plays a key role in the cell-matrix interactions by promoting matrix assembly and cell adhesiveness. Especially, SPARC is a key matricellular protein involved in collagen I deposition and fibrillogenesis [29]. Since the LW exhibited higher muscle collagen content than the B pigs, SPARC gene expression could mediate this discrepancy. In this cluster, DCN (decorin) is involved in matrix assembly, and its targeted ablated expression strongly affects the collagen network [52]. DCN is also known to interact with TGFB1 in satellite cells proliferation and differentiation [53]. Interestingly, TGFB receptor was found in the same cluster, suggesting that this interaction could explain LW muscle development. DCN, dermatopontin (DPT) and dystonin (DST) act in the same way since DPT accelerates the assembly of collagen into fibrils [54], whereas DST deficient mice exhibited weak skeletal muscle cytoarchitecture [55]. Moreover, six genes encoding various collagen types are highly expressed in LW muscle, in agreement with their higher collagen content. Thus, all these biological processes are in accordance with the LM properties of the LW compared with the B pigs, namely their elevated collagen content and shear force value, and lower tenderness score [56].

Vacuole, lysosome and proteolysis
The vacuole and lysosome cluster (cluster 2-B) is a highly enriched CC cluster in the B pigs. Lysosomes contain many hydrolytic enzymes involved in the degradation of cytoplasmic proteins, even if the calpains and proteasome represent the main myofibibrillar proteolysis pathways [57]. Cathepsin D (CTSD, lysosomal protein) is highly expressed in the B pigs, and a mutation in CSTD gene has been associated with increased average daily gain and muscle mass, and decreased backfat deposition in both Duroc and LW pigs [58]. In this cluster, we also found ATP6V1D (ATPase, H+ transporting lysosomal), a subunit of a vacuole ATPase pumping protons from the cytoplasm to the lumen of the lysosome which might control pH homeostasis in muscle cells [59]. The potential role of NEU1 (sialidase 1) in the control of cell proliferation, collagen content and extracellular matrix remodelling in skeletal muscle [60] and in inhibition of early myogenesis [61], agrees with the lower expression of NEU1 found in LW pigs. This is also in accordance with their higher expression of collagen encoding genes and muscle development, and might thus be related to their lower meat tenderness.
The proteolysis function is not put forward by enrichment cluster analyses, but is at upmost importance in the context of meat tenderness. CAST (calpastatin), an inhibitor of calpain, one of the main proteolytic enzymatic systems involved in post-mortem muscle proteolysis and tenderization [62], was found highly expressed in LW breed. This could explain the higher shear force and lower tenderness score of the LW, as a result of a lower proteolysis level. Indeed, CAST has been suggested as a candidate gene for meat tenderness in pigs [63]. The ubiquitin-proteasome system, the main actor of nonlysosomal cytoplasmic protein degradation, appears to be also involved in the discrepancy between B and LW breeds with two out of the three ligase enzymes of the proteasome complex, TRIM63 and FBXO32 (Fbox protein 32), more expressed in the B pigs. This could also contribute to the more tender meat of the B pigs, since proteasome would be one of the main endogenous proteolytic systems contributing to post-mortem meat tenderization [64].
In conclusion, our study aimed at identifying the biological events underlying differences in muscle physiology and meat quality traits reported in literature between the contrasted B and LW breeds. From transcriptomics and functional analyses, four main biological clusters were identified. Energy metabolism and lipid deposition are associated with breed-differences in muscle gene expressions and chemical composition. Furthermore, the cytoskeleton and the contractile fibers would play a role in the determination of muscle and meat phenotypes. Last, our results suggest the extracellular matrix as an important component of the LW muscle in accordance with their elevated collagen content, which could explain their reduced tenderness.
As a whole, our results contribute to a better understanding of muscle physiology and its consequences on the development of meat quality. Besides, this study is a first step towards the identification of molecular markers of pork quality and the subsequent development of control tools.

Ethics Statement
The experiment was conducted following French guidelines for animal care and use edited by the French Ministries of High Education and Research, and of Agriculture and Fisheries (http:// ethique.ipbs.fr/sdv/charteexpeanimale.pdf). All animals were reared and slaughtered in compliance with national regulations and according to procedures approved by the French veterinary Services at INRA PEGASE facilities. Our research unit was holder of a pig experimentation agreement (Nu A35622) delivered by the Veterinary Services of the French Ministry of Agriculture. Moreover, the technical and scientific staff involved in the experiment was holder of an individual agreement for experimentation on living animals delivered by the French Veterinary Services.

Animals, husbandry and slaughtering
Forty finishing castrated males from a commercial selected LW pure breed (n = 20 issued from 10 litters produced from 9 different boars) pigs and an autochthonous B breed (n = 20 issued from 10 litters produced from 6 different boars) were reared in two different housing systems. At the average live weight of 35 kg, 2 pigs from each litter were chosen on the basis of their live weight and growth rate from birth up to 35 kg, and assigned to either a conventional (C), indoors on slatted floor (1.0 m 2 /pig), or an alternative (A) with indoor bedding (1.3 m 2 /pig) and a free access to an outdoor area (1.1 m 2 /pig) housing system at INRA-UMR PEGASE experimental farm, thus giving 4 groups of 10 pigs per breed and housing system. Pigs of both housing systems were fed a standard commercial diet with 2.5 kg/d/pig from 35 up to 110 kg, and 3.0 kg/d/pig up to slaughter at the average weight of 145 kg. Pigs were slaughtered at INRA-UMR PEGASE experimental slaughterhouse in four sessions (each including pigs from each breed and housing system), by electrical stunning and exsanguination.

Carcass, muscle and meat quality measurements
The day of slaughter, hot carcass weight and back fat thickness (mid line, between 4 th and 5 th lumbar vertebra level) were recorded. After 24 h at 4uC, the fresh carcass and wholesale cuts of right side were weighted for calculation of loin and backfat proportions. Thirty minutes after exsanguination, a sample of LM was carefully collected on all pigs (right half-carcass, last rib level) immediately frozen in liquid nitrogen and stored at 280uC until RNA extraction (see below) and determination of glycolytic potential [13]. The following day, a transversal slice of LM was taken (1 st lumbar vertebra level), trimmed of external fat and minced. Half of this sample was put under vacuum and stored at 220uC until determination of intramuscular fat content, and the remaining minced LM was freeze-dried before determination of protein and collagen contents, as previously described [65]. Water content was determined from the weight of minced muscle before and after freeze-drying, and used for calculations of protein and collagen content per gram of fresh muscle. The day after slaughter, another slice (1.5 cm depth) was taken of the LM (the 2 nd vertebra level), weighed (100610 g) and kept for 48 h at 4uC in plastic bags for determination of drip loss [66]. On all pigs, a piece of the left loin (between the 2 nd /3 rd and 6 th /7 th last ribs) was taken 24 h after slaughter, partially trimmed of external fat, kept at 4uC for 3 subsequent days, and deboned. The LM was put under vacuum, frozen and stored at 220uC before determination of shear force on cooked meat using a Warner-Bratzler cell fitted on a universal testing machine (Instron France S.A.S., Guyancourt, France) according to Honikel [66]. The shear force mean values were obtained from 10 measurements per LM sample. A piece of the right loin (between the 2 nd /3 rd and 9 th /10 th last ribs) was also taken on all pigs the day after slaughter, prepared and stored as described above until sensory analysis performed at INRA-EASM (Le Magneraud, Surgères, France). The 40 roasts were evaluated over 10 sessions, each including four roasts, one per breed and per rearing system. After thawing for 48 h at 4uC, roasts (900 g) were cooked in an oven (dry heat, 250uC for 10 min, followed by humid heat, 100uC for around 45 min up to a core temperature of 8062uC). Then, the middle part of 1-cm thick slices of roasts was presented to the 12 panellists who evaluated tenderness, juiciness and flavour on a continuous scale form 0 (absent) to 10 (very high intensity of the trait). The average of individual panellist scores from each sample was used for the statistical analysis.
Carcass, muscle and meat quality data were submitted to an analysis of variance (GLM procedure, SAS Inst. Inc., Cary, NC). The model included the fixed effects of breed, housing system, and their interaction. Least square (LS) means were calculated per breed and per housing system.

Microarray design
A custom pig skeletal muscle microarray [14] of 15198 oligonucleotides (60 mers) was used in this study. Among the 15198 probes of the GenmascqChip, 12939 probes (i.e. 85% of the oligonucleotides) have been linked to a unique annotated sequence and to 9169 unique genes (i.e., 30% of redundancy). An 8615 K oligo-microarray Agilent format was chosen, therefore one probe per microarray and eight microarrays were fitted in each slide.

RNA extraction and Microarray hybridization
Total RNA was extracted by crushing the frozen tissue in Trizol reagent (Invitrogen, Cergy-Pontoise, France) and purification using a silica-based spin-column (RNA II kit, Macherey Nagel, Lyon, France). The quality and concentration of total RNA were verified by electrophoresis using an Agilent Analyzer (Agilent Technologies France, Massy, France) and UV spectrometry (Nanodrop, Thermo Scientific, Illkirch, France). In order to compare the 40 LM samples among the experiment, each sample was compared to a reference pool composed of an equal amount of transcripts isolated from all 40 LM samples. Total RNA (350 ng) from each animal was labeled individually with Cy3 and the reference sample was labeled with Cy5, using the Quick-Amp Labeling Kit (Agilent Technologies, Santa Clara, USA) and following the manufacturer's instructions. Microarray hybridizations were carried out in Agilent's SureHyb Hybridization Chambers containing 300 ng of Cy3-labeled cRNA sample and 300 ng of Cy5-labeled reference sample per hybridization. The hybridization reactions were performed at 65uC for 17 hours using Agilent's Gene Expression Hybridization Kit. Slides were disassembled and washed in Gene Expression Wash Buffer 1 for 1 minute at room temperature and then in Gene Expression Wash Buffer 2 for 1 minute. Microarrays were scanned at 5 mm/pixel resolution using the Agilent DNA Microarray Scanner G2505B, and images were analyzed with Agilent Feature Extraction Software (Version 9.5), using the GE2-v5_95_Feb07 FE extraction protocol. These MIAME compliant microarray data have been deposited into the GEO [15] repository and are publicly available through GEO Series accession no. GSE26614.

Microarray Data Analyses and Statistics
All analyses were performed using the R software version 2.8.1 [67]. Raw spots intensities were first submitted to quality filtration based on four criteria: intensity, uniformity, saturation and outliers detection. Intensities of filtered spots were transformed into log 2 (Cy3/Cy5). Data were normalized within chips by subtraction of the sample median value across all probes from all raw values, and between chips using the ''Rquantile'' method of the Limma R package [68] to obtain experimentally consolidated gene expression values. The ''Rquantile'' method was used since the red channel was the common reference throughout the experiment. To increase the power of differential expression analysis [69], spots with the smallest expression variability across samples were filtered out using K-means algorithm (k = 4). All together, 4870 spots were finally retained for statistical analyses. Expression data were then adjusted for slide effect when significant (p,0.05) by analysis of variance, before performing differential expression analysis. Residuals were then submitted to an analysis of variance using the fixed effects of breed (B or LW), housing system (C or A) and their interaction (breed6housing system). Data were then submitted to Benjamini and Hochberg (BH) multiple testing correction procedure [70] using an adjusted P value cutoff of 0.05.

Functional analysis
Enrichment analysis for specific GO terms for BP has been carried out using the DAVID [17][18][19]. In DAVID analysis, the GO _FAT terms were selected to filter the broadest terms without overshadow the more specific ones. The lists of genes were uploaded using the ENTREZ gene ID (http://www.ncbi.nlm.nih. gov/gene). The P values for enrichment (or EASE scores) were computed by a modified Fisher's exact test, using our custom microarray (i.e. 8639 human ENTREZ gene ID) as background. The GO categories (BP, CC and MF) and KEGG and Biocarta pathways were clustered using the DAVID Functional Annotation Clustering tool [17][18][19], where the enrichment score for each cluster was computed as the negative log of the geometric mean of P values in the cluster.

Real Time PCR analysis
Complementary DNA was synthesised from 2 mg of total RNA previously used for microarray analysis, using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Foster City, CA). Primers were designed using Primer Express Software (Applied Biosystems, USA) based on Sus scrofa published nucleotides sequences (Table S3). Amplification was performed in triplicate, in 12.5 ml with 2.5 ng of reverse-transcribed RNA and both forward and reverse primers (200 nM each) in 16 PCR buffer (Fast SYBRH Green Master Mix, Applied Biosystems) containing Uracil DNA glycosylase to prevent any DNA contamination from previous PCR. A StepOnePlus TM Real Time PCR system (Applied Biosystems) was used. Thermal cycling conditions were as follows: 50uC for 2 min, 95uC for 20 s, followed by 40 cycles of denaturation at 95uC for 3 s, and annealing at 60uC for 30 s. Specificity of the amplification products was checked by dissociation curves analysis. Three genes were used as reference for normalization: HPRT1 (hypoxanthine phosphoribosyltransferase 1), B2M (beta-2 microglobulin) and 18S (18S rRNA predeveloped TaqMan kit from Applied Biosystems). Using geNorm [71] and Normfinder [72] algorithms, all three genes appeared to have a stable expression on all LM samples. For each sample, the normalized expression level (N exp ) was calculated according to the following formula: N exp = (1+E) 2DCt (sample2calibrator) /NF, where the calibrator is a pool of the 40 LM samples, E is the PCR efficiency and NF is the normalization factor calculated using geNorm algorithm. Normalized expression levels of mRNAs were then compared between B and LW samples using the Student t-test and P value#0.05 for significance.

Supporting Information
Table S1 Genes highly expressed in Basque pigs. Results were expressed as the Basque to Large White ratio of the gene expression. The p value of each gene was adjusted according to the Benjamini-Hochberg method. Difference in gene expression was considered significant if its adjusted p value was p,0.05. Redundancy represented the number of probes per gene. In this list, 73 genes had more than one probe. (XLS) Results were expressed as the Basque to Large White ratio of the gene expression. The p value of each gene was adjusted according to the Benjamini-Hochberg method. Difference in gene expression was considered significant if its adjusted p value was p,0.05. Redundancy represented the number of probes per gene. In this list, 75 genes had more than one probe. (XLS)